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