Problem 697: Randomly Decaying Sequence

Project Euler Problem 697 Solution

EulerSolve provides an optimized solution for Project Euler Problem 697, Randomly Decaying Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

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

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

Problem 697 source code

C++

#include <boost/math/special_functions/gamma.hpp>

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>

namespace {

long double log10_c_for_n(const long double n) {
    const long double ln_c = boost::math::gamma_p_inv(n, 0.75L);
    return ln_c / std::log(10.0L);
}

}  // namespace

int main() {
    const long double ln_c_100 = boost::math::gamma_p_inv(100.0L, 0.75L);
    const long double log10_c_100 = ln_c_100 / std::log(10.0L);
    assert(std::fabsl(log10_c_100 - 46.27L) < 0.01L);
    assert(std::fabsl(boost::math::gamma_q(100.0L, ln_c_100) - 0.25L) < 1e-15L);

    const long double n = 10'000'000.0L;
    const long double ln_c = boost::math::gamma_p_inv(n, 0.75L);
    assert(std::fabsl(boost::math::gamma_q(n, ln_c) - 0.25L) < 1e-12L);

    const long double answer = ln_c / std::log(10.0L);
    std::cout << std::fixed << std::setprecision(2) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    # We want to find c such that GammaCDF(c, n) = 0.75 for n = 10^7.
    # The Wilson-Hilferty transformation provides an excellent approximation
    # for large n: c ≈ n * (1 - 1/(9n) + z * sqrt(1/(9n)))^3
    # where z is the 75th percentile of the standard normal distribution.
    
    n = 10000000.0
    # z = Phi^{-1}(0.75)
    z = 0.6744897501960817
    
    term = 1.0 - 1.0 / (9.0 * n) + z / math.sqrt(9.0 * n)
    c = n * (term ** 3)
    
    ans = c / math.log(10.0)
    
    return f"{ans:.2f}"

if __name__ == '__main__':
    print(solve())

Java

import java.util.Locale;

public class Euler697 {
    public static String solve() {
        double n = 10000000.0;
        double z = 0.6744897501960817;

        double term = 1.0 - 1.0 / (9.0 * n) + z / Math.sqrt(9.0 * n);
        double c = n * Math.pow(term, 3);

        double ans = c / Math.log(10.0);

        return String.format(Locale.US, "%.2f", ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}