# Finding the probability of one or more very unlikely successes in an enormous number of Bernoulli trials

I'm trying to find an upper bound on the complexity of a self-replicating matter pattern in order for it to be at all probable for it to spontaneously occur in the universe. The idea being, if some day in the future it can be proven that a self-replicating blob of matter cannot be composed of fewer than this many atoms, the dawn of life is by extension proven to be the great filter in the Fermi Paradox.

To set this upper bound I'm assuming that every atom in the universe (1080 atoms) has had one try to find the correct configuration once per nanosecond since the dawn of the universe (4.3 ✕ 1026 nanoseconds).

The number of matter permutations that you have to sort through is the factorial of the minComplexity, which I want to solve for. For example, I want to find what complexity yields a probability of at least one success in the entire history of the universe ≥ 0.5.

So the probability of getting it right for a single Bernoulli trial is 1/(minComplexity!)

The direct way to tackle this is to find the probability of getting all the trials wrong and inverting the result, which results in this Mathematica code:

1 - (1 - 1/minComplexity!)^(4.3 10^26 10^80)


However, I can't evaluate this for any values of minComplexity (at least, none that I've tried so far), instead getting General: Overflow or General: Underflow errors as it tries to raise a fraction to an extremely high power (driving it very close to zero). Forcing the precision with SetPrecision to ∞ doesn't make the problem go away. Solve and NSolve are also unable to handle the exponentiation as an equality or inequality (even if I replace the factorial with the gamma function).

Hoping that Mathematica's BinomialDistribution function would have more robust handling of numeric extremes in p and n, I tried

BinomialDistribution[Round[4.3 10^26 10^80], 1/minComplexity!]


and still encounter General: Underflow errors when I try to get the PDF out of it to evaluate with various values for minComplexity.

I suspect if I knew a value of minComplexity that causes the end probability to be a number far from both 0 and 1 it would work, but finding that value is the very exercise in question.

I believe in both cases, the problem boils down to raising a fraction in the interval (0,1) to an extremely large power. What is the best way to work around this in Mathematica?

P.S. The current leader for least-complex self-replicating blob of matter has 531K base pairs of DNA, not counting proteins and other components. I don't need Mathematica to know that 531K! is impossibly huge even with this generous upper bound, so if this 531K is proven to be as small as it can go we've found the great filter.

• When numbers are too large it is sometimes feasible to consider working with the Log of your expression. This johndcook.com/blog/2010/08/16/how-to-compute-log-factorial might be helpful.
– Bill
Dec 28, 2016 at 13:25
• What do you think of assuming independence of events among those 10^80 atoms? At some point should you not consider some model that does not assume independence? (I'm not suggesting that there is an obvious model of non-independence or that it would be easy. I'm just noting that assumption rarely seems questioned.)
– JimB
Dec 28, 2016 at 15:35
• With minComplexity = 73.6110 you get approximately 0.5 for 1 - (1 - 1/minComplexity!)^(43 * 10^25 * 10^80) Dec 28, 2016 at 15:55
• @JimBaldwin For an upper bound I think it's reasonable. If there is some easily-formed intermediate structure that is not self-replicating but reduces the complexity requirement for getting to a self-replicating structure, there are still guaranteed to be fewer instances of said structure than there are atoms in the universe, so considering all the atoms in the universe remains a valid upper bound. Dec 28, 2016 at 20:59

(4.3 10^26 10^80)*Log10[1 - (1 - 1/513000!)]
The result is: -1.1638*10^113`