Calculating very complex probability equation

edit

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\cdot 10^{96}$, and let $k$ be the number $16\cdot 10^9$. You want to calculate: $$ A = \frac {N\cdot(N-1)\cdot (N-2)\cdot\dots\cdot(N-k+1)} {N\cdot N\cdot N\cdot\dots\cdot N} =1 \cdot \left(1-\frac 1N\right)\cdot \left(1-\frac 2N\right)\cdot\dots\cdot \left(1-\frac {k-1}N\right)\ . $$ We can try to estimate and approximate. For instance, the product is between $1$ and $\left(1-\frac kn\right)^k$, and for the last power use Bernoulli in order to not use the computer. So a lower bound is $1-\frac {k^2}N$.

But we can do better: $$ -\log(1-x) = x+\frac 12x^2+\frac 13x^3+\frac 14x^4+\dots $$ so that it is for values of $x$ in our range, $|x|\le k/N$, between $x$ and $x+x^2$. (We may want to refine later, if this does not work.) Then: $$ e^{-x-x^2}\le (1-x)\le e^{-x}\ . $$ So our product is between the following numbers: $$ \begin{aligned} A_0 &:=\exp-\sum_{1\le j< k}\frac jN +\frac {j^2}{N^2} =\exp\left(-\frac 1N k(k-1)-\frac 1{N^2}k(k-1)(2k-1)\right)\ ,\\ A_1 &:=\exp-\sum_{1\le j< k}\frac jN =\exp\left(-\frac 1N k(k-1)\right) \ . \end{aligned} $$ 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 $A_1$ is in fact the exponential of the tiny number

sage: -k*(k-1)/N
-6.56410256369231e-78

and in order to also have $A_0$ we need the contribution of the tiny-tiny number...

sage: -k*(k-1)*(2*k-1)/N^2
-5.38593030850230e-165