# Calculating very complex probability equation

I want to calculate the probability of following statement. It will be very close to $1$ ,but I need to see how close it is. I could not calculate it using standart python, so I need very sensitive tool. I think that sage can handle it. However, I could not take proper solution. Maybe, i am missing something. $$\frac{(39 \times 10^{96})!}{(39 \times 10^{96}-16 \times10^9)! \times (39 \times 10^{96})^{16 \times10^9}}$$

edit retag close merge delete

( 2023-11-08 15:13:20 +0100 )edit

@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 ?

( 2023-11-08 16:20:43 +0100 )edit

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.

( 2023-11-08 18:15:55 +0100 )edit

it is pen and paper problem with mandatory computation

( 2023-11-08 20:03:45 +0100 )edit

Sort by ยป oldest newest most voted

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

more

Using numerical computations with 100 decimals precision in GP/Pari one can check that the product A is equal

0.9999999999999999999999999999999999999999999999999999999999999999999999999999967179487158974358974359

( 2023-11-21 16:34:51 +0100 )edit