The first rule of a computer algebra system is that the human must go humaly as far as possible, the program a humanly meaningful path and result. In our case...
Let N be the number 39⋅1096, and let k be the number 16⋅109. You want to calculate:
A=N⋅(N−1)⋅(N−2)⋅⋯⋅(N−k+1)N⋅N⋅N⋅⋯⋅N=1⋅(1−1N)⋅(1−2N)⋅⋯⋅(1−k−1N) .
We can try to estimate and approximate.
For instance, the product is between 1 and (1−kn)k, and for the last power use Bernoulli in order to not use the computer. So a lower bound is 1−k2N.
But we can do better:
−log(1−x)=x+12x2+13x3+14x4+…
so that it is for values of x in our range, |x|≤k/N, between x and x+x2. (We may want to refine later, if this does not work.) Then:
e−x−x2≤(1−x)≤e−x .
So our product is between the following numbers:
A0:=exp−∑1≤j<kjN+j2N2=exp(−1Nk(k−1)−1N2k(k−1)(2k−1)) ,A1:=exp−∑1≤j<kjN=exp(−1Nk(k−1)) .
It is time to use the computer:
k, N = 16. * 10^9, 39. * 10^96
A1 = exp( -k*(k-1)/N )
A0 = A1 * exp( -k*(k-1)*(2*k-1) / N^2 )
print(f"A0 = {A0}")
print(f"A1 = {A1}")
And we get of course:
A0 = 1.00000000000000
A1 = 1.00000000000000
Because A1 is in fact the exponential of the tiny number
sage: -k*(k-1)/N
-6.56410256369231e-78
and in order to also have A0 we need the contribution of the tiny-tiny number...
sage: -k*(k-1)*(2*k-1)/N^2
-5.38593030850230e-165
See https://en.wikipedia.org/wiki/Stirlin...
@maxalekseyev, stirling formula is good for factorial, but sage could not calculate the exponential expression in the denominator. If you can calculate them,can you please share your code with me ?
You don't need to calculate it explicitly to get an estimation. It sounds more like a homework pen-and-paper problem than a computational one.
it is pen and paper problem with mandatory computation