First digit of a factorial without calculating it

Question

I made this small function which use Stirling formula (https://en.wikipedia.org/wiki/Stirling%27s_approximation) and the fractional part of log10 to calculate the first digit of a factorial. However I am not sure of its correctness, I know it will likely fail for low number but if we assume n>1000 then I want to know which case it would fail. Here the gist: https://gist.github.com/hube12/ee3330b6d47a5c87a505a9f81e353b7a

def first_digit(n):
    from math import exp,log,pi,modf,sqrt
    w=log((n/exp(1)),10)
    f,i=modf(w*n)
    l=10**f
    p=sqrt(2*pi*n)
    return str(p*l)[0]

https://wikimedia.org/api/rest_v1/media/math/render/png/7fe20ccef4b13b2fc2b79b752fb595da6d855de2

Let's take the time to see what I did in math formula:

Stirling (assuming large n) : $ n! \approx \sqrt{2\times\pi \times n} \times (\frac{n}{e})^n $

Now I take the log base 10 of the fraction so $ w= \log_{10}{\frac{n}{e}} $

This allow me then to take the decimal part of $ w \times n $

So I used $ \log_{10}{(\frac{n}{e})^n} = n \times \log_{10}{(\frac{n}{e})} $

Now that I have the decimal part I can raised 10 to that part, $ l= 10^f$ This allow me to get a rational number in [0,10) which represent the highest digit of that power.

Now I only need to multiply that number l by the square root part and get the first digit of that product $ \sqrt{2\times\pi \times n} \times FirstDigit((\frac{n}{e})^n) $

The first digit of that product is the first digit of my factorial, well that's what I am asking if it is correct for every factorial (assuming n>1000 and infinite precision/ correct log)

Cheers, Neil

Answer 1

The formula you need is: $$\log(N!)=N \log(N) - N \log(e) + \log(2 \pi N)/2 + \log(1+1/(12N)) $$

The below gives me 15 correct digits for $10^9!$:

def factorial(N):
    N=Decimal(N)
    p=N*N.log10() - N/Decimal(10).ln() + (2*N*Decimal(math.pi)).log10()/2
    d=int(p)
    p-=d
    return (10**p * (1+1/(12*N)),d)

Your formula fails to work since you are using math.e which has insufficient precision. Also, you won't do better than 9 digits unless you add in the second term in the series (1/12N),

Additionally, you don't need str throughout your code, it just makes it messy.