Detailed mathematical approach

Problem Summary

Let \(X_0=1\) and let the sequence decay by independent random factors \(U_1,U_2,\dots\), each uniformly distributed on \((0,1)\), so that

$$X_n=\prod_{k=1}^{n} U_k.$$

The target constant \(c\) is defined by the probability condition

$$\Pr\!\left(X_n\le \frac{1}{c}\right)=0.25\qquad\text{for }n=10^7,$$

and the required output is \(\log_{10} c\). The key observation is that products of independent uniform variables become sums after taking logarithms, which turns the problem into a gamma-distribution quantile computation.

Mathematical Approach

The implementations reduce the random-product statement to a single equation involving the regularized incomplete gamma function. Once that equation is inverted, the answer follows immediately.

Step 1: Convert the product into a sum

Take minus logarithms and define

$$S_n=-\ln X_n=-\sum_{k=1}^{n}\ln U_k=\sum_{k=1}^{n}(-\ln U_k).$$

Therefore the event in the problem becomes

$$X_n\le \frac{1}{c}\iff -\ln X_n\ge \ln c\iff S_n\ge \ln c.$$

So the original probability condition is equivalent to

$$\Pr(S_n\ge \ln c)=0.25.$$

Step 2: Identify the distribution of each logarithmic increment

If \(U\sim \operatorname{Uniform}(0,1)\), then \(Y=-\ln U\) has density

$$f_Y(y)=e^{-y}\qquad (y\ge 0),$$

which is the exponential distribution with rate \(1\). Because the factors are independent, the variables \(Y_1,\dots,Y_n\) are independent as well, and hence

$$S_n=Y_1+\cdots+Y_n\sim \operatorname{Gamma}(n,1),$$

meaning shape \(n\) and scale \(1\).

Step 3: Rewrite the probability with regularized gamma functions

For a gamma variable with shape \(n\) and scale \(1\), the lower and upper regularized incomplete gamma functions are

$$P(n,x)=\frac{\gamma(n,x)}{\Gamma(n)},\qquad Q(n,x)=\frac{\Gamma(n,x)}{\Gamma(n)}=1-P(n,x).$$

They satisfy

$$\Pr(S_n\le x)=P(n,x),\qquad \Pr(S_n\ge x)=Q(n,x).$$

Therefore the defining condition for \(c\) becomes

$$Q(n,\ln c)=0.25,$$

or equivalently

$$P(n,\ln c)=0.75.$$

When \(n\) is a positive integer, the upper tail can also be written as

$$Q(n,x)=e^{-x}\sum_{k=0}^{n-1}\frac{x^k}{k!},$$

but for \(n=10^7\) this identity is conceptually useful rather than computationally attractive.

Step 4: Recognize the answer as a gamma quantile

Set

$$x=\ln c.$$

Then \(x\) is exactly the 75th percentile of the \(\operatorname{Gamma}(n,1)\) distribution:

$$x=P^{-1}(n,0.75).$$

Once \(x\) is known, the requested quantity is not \(c\) itself but its base-10 logarithm, so

$$\log_{10} c=\frac{\ln c}{\ln 10}=\frac{x}{\ln 10}.$$

This is numerically important because \(c\) is astronomically large, while \(\log_{10} c\) is easy to store and print.

Step 5: Large-\(n\) approximation used by two implementations

The C++ implementation computes the gamma quantile directly. The Python and Java implementations instead use the Wilson-Hilferty approximation for a large gamma quantile:

$$x\approx n\left(1-\frac{1}{9n}+\frac{z}{3\sqrt{n}}\right)^3,$$

where

$$z=\Phi^{-1}(0.75)\approx 0.6744897501960817.$$

Because \(n=10^7\) is enormous, this approximation is already very sharp for a result that is ultimately printed to two decimal places.

Worked Example: the checkpoint \(n=100\)

Using the same approximation with \(n=100\), we get

$$x\approx 100\left(1-\frac{1}{900}+\frac{0.6744897502}{30}\right)^3\approx 106.52.$$

Hence

$$\log_{10} c\approx \frac{106.52}{\ln 10}\approx 46.27.$$

This matches the numerical checkpoint used by the exact computation to two decimal places. It also explains why the asymptotic formula is trusted for the much larger target value \(n=10^7\).

How the Code Works

All three implementations first translate the probability question into the gamma-quantile problem \(P(n,x)=0.75\) with \(x=\ln c\), and then return \(x/\ln 10\).

The C++ implementation evaluates that quantile directly with a special-function routine for the inverse lower regularized incomplete gamma function. It also performs two sanity checks: one at \(n=100\), where the expected printed value is \(46.27\) to two decimals, and one at \(n=10^7\), where substituting the computed \(x\) back into the upper regularized gamma function must recover \(0.25\) to high precision.

The Python and Java implementations take the asymptotic route. They hard-code the 75th percentile of the standard normal distribution, plug it into the Wilson-Hilferty cubic approximation for the gamma quantile, and finally divide by \(\ln 10\). That replaces a special-function inversion by a short sequence of elementary floating-point operations.

Complexity Analysis

For the single target input \(n=10^7\), all implementations use \(O(1)\) memory. The Python and Java versions run in \(O(1)\) time with only a handful of arithmetic operations. The C++ version is also \(O(1)\) at the algorithmic level for a fixed \(n\), but its constant factor is larger because it performs a numerical inversion of a special function and then verifies the result with another special-function evaluation.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=697
2. Incomplete gamma function: Wikipedia - Incomplete gamma function
3. Gamma distribution: Wikipedia - Gamma distribution
4. Exponential distribution: Wikipedia - Exponential distribution
5. Wilson-Hilferty transformation: Wikipedia - Wilson-Hilferty transformation