The sum as described in the Wikipedia article has a closed form. It took me some time but after some twiddling with the sum itself I decided to exponentiate the whole thing to get rid of the logarithm and had a simple product I was able to handle.
With the digit and the base:
This can be expressed as a fraction of two rising factorials (Pochhammer symbol).
The rising factorial can be expressed as a fraction of two Gamma functions
With which we get
By setting we can already see that the result will get smaller the larger gets with the obvious limit
The arguments of the Gamma functions get quite large which makes the results hard to handle but we already need the logarithm of the product, that means we can use the logarithm of the Gamma function and save one step (we are on the real-line, so I used the equality to avoid an even denser parentheses fence than it already is)
We can simplify a bit with because . Furthermore which can save us a multiplication here and there. So far we got
In the special case I used it for was a power of two such that none of the exponentiations had to be calculated, the numbers could be build in O(1).
Calling the whole thing we can compute the probability that the nth digit is with . For general use the unsimplified version above for .
Beware: the numbers get really large, you will need a slightly higher amount of precision for the loggamma function if you calculate with a large base. I had 226 as the base and had trouble to compute the values for places higher than the 20th in PARI/GP (it said
*** lngamma: division by zero). It worked in calc with a precision of 300 decimal digits but it is sloooow. Who wrote that lame gamma function? Who? Me? Myself? Really? Ouch! 😉