Problem 197: Investigating the Behaviour of a Recursively Defined Sequence

Project Euler Problem 197 Solution

EulerSolve provides an optimized solution for Project Euler Problem 197, Investigating the Behaviour of a Recursively Defined Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The sequence is defined by $$u_{0}=-1,\qquad u_{n+1}=f(u_n),\qquad f(x)=\lfloor 2^{30.403243784-x^2}\rfloor\cdot10^{-9}.$$ We are asked for $u_n+u_{n+1}$ when $n$ is very large. The important point is that this is not a symbolic closed-form problem; it is a discrete dynamical-system problem with a heavily quantized update rule. Mathematical Approach The entire solution comes from understanding what the floor and the factor $10^{-9}$ do to the orbit. They turn a nonlinear recurrence into a bounded map on a finite grid, and for the given initial value that orbit rapidly falls onto an attracting 2-cycle. A bounded nonlinear map on a decimal lattice Because $x^2\ge 0$, the exponent $30.403243784-x^2$ is never larger than $30.403243784$. Therefore $$0\le f(x)\le \lfloor 2^{30.403243784}\rfloor\cdot10^{-9}=1.420000000.$$ So after the first iteration, every term of the sequence lies in the compact interval $[0,1.42]$. Even more importantly, the floor forces every value onto the grid $$\{k\cdot10^{-9}:0\le k\le 1{,}420{,}000{,}000\}.$$ That means the recurrence is not wandering through a continuum of real numbers. It is moving through a finite set of decimal values spaced exactly $10^{-9}$ apart. Finite state space implies eventual periodicity Once a deterministic map sends a finite set into itself, every orbit must eventually repeat....

Detailed mathematical approach

Problem Summary

Let $c=30.403243784$ and define, for every real number $x$,

$$f(x)=\bigl\lfloor 2^{\,c-x^{2}}\bigr\rfloor\cdot10^{-9},$$

where $\lfloor y\rfloor$ is the floor of $y$, the largest integer not exceeding $y$. The sequence starts with $u_0=-1$ and continues with $u_{n+1}=f(u_n)$. We must find $u_n+u_{n+1}$ for $n=10^{12}$, given with $9$ digits after the decimal point. The statement gives no intermediate values. ^[1]

The whole solution in one sentence: the terms with odd index can only fall and those with even index can only rise, and since all terms are nine-decimal numbers both runs must stop, after which the sequence swaps the same two numbers forever; this happens after $516$ steps, so iterating a thousand steps and adding the last two terms gives the sum for $n=10^{12}$.

Mathematical Approach

Sections 1 to 4 prove that the sequence becomes periodic with period $2$ long before $n=10^{12}$, which gives the algorithm. Sections 5 to 7 explain, through the smooth curve behind $f$, where it settles, why neighbouring cycles give the same sum and why settling takes a few hundred steps; the exact moment comes from computing the rounded sequence.

1. A bell curve rounded down

Put $A=2^{c}\cdot10^{-9}$ and $F(x)=A\cdot2^{-x^2}$. Then $2^{c-x^2}=10^{9}F(x)$, and the definition becomes

$$f(x)=10^{-9}\,\bigl\lfloor10^{9}F(x)\bigr\rfloor.$$

So $f(x)$ is $F(x)$ cut off after its ninth decimal. The function $F$ is a bell curve: it is even, largest at $x=0$ and decreasing for $x\ge0$. The floor satisfies $y-1\lt\lfloor y\rfloor\le y$, ^[2] hence

$$F(x)-10^{-9}\lt f(x)\le F(x).$$

A high-precision evaluation gives $2^{c}=1420000000.279\ldots$ and $2^{c-1}=710000000.139\ldots$, so $A=1.4200000002793\ldots$, $f(0)=1.42$ and $f(\pm1)=0.71$: the constant $c$ makes these two values round. Every value of $f$ is therefore a multiple of $10^{-9}$ between $0$ and $1.42$, and so is every term from $u_1=f(-1)=0.71$ on. The sequence lives on a finite grid of $1\,420\,000\,001$ numbers.

2. The first terms

The first steps can be checked with a calculator. Since $0.71^2=0.5041$ and $2^{c-0.5041}=1001242148.918\ldots$, we get $u_2=1.001242148$. Then $u_2^2=1.002485838\ldots$ and $2^{c-u_2^2}=708777686.401\ldots$, so $u_3=0.708777686$. Continuing in the same way:

$u_1=0.710000000$, $u_2=1.001242148$;
$u_3=0.708777686$, $u_4=1.002446415$;
$u_5=0.707593212$, $u_6=1.003612800$;
$u_7=0.706446531$, $u_8=1.004741400$;
$u_9=0.705337499$, $u_{10}=1.005832407$.

The terms with odd index are small and fall; those with even index are large and rise, each by about $0.001$ per pair of steps. Sections 3 and 4 show that this pattern can never reverse and must end.

3. Two steps at a time

For $0\le x\le y$ we have $c-x^2\ge c-y^2$, hence $2^{c-x^2}\ge2^{c-y^2}$, and since the floor is non-decreasing, ^[2] $f(x)\ge f(y)$. On $[0,\infty)$ the map $f$ reverses order. All its values are non-negative, so the double step $g=f\circ f$ reverses the order twice, that is, it preserves it:

$$0\le x\le y\ \Longrightarrow\ g(x)\le g(y).$$

Claim. $u_1\ge u_3\ge u_5\ge\cdots$ and $u_2\le u_4\le u_6\le\cdots$.

Proof. All terms from $u_1$ on are non-negative, and $u_{n+2}=g(u_n)$. Section 2 gives $u_3\le u_1$. If $u_{2j+3}\le u_{2j+1}$, applying $g$ gives $u_{2j+5}\le u_{2j+3}$, so by induction the odd-indexed terms never increase. Applying $f$, which reverses order, to $u_{2j+3}\le u_{2j+1}$ gives $u_{2j+4}\ge u_{2j+2}$. $\square$

In particular every odd-indexed term is at most $u_1=0.71$, and every even-indexed term is at least $u_2=1.001242148$.

4. The sequence settles long before $10^{12}$

Write $u_{2j+1}=k_j\cdot10^{-9}$. By Section 3 the integers $k_j$ satisfy $710\,000\,000=k_0$${}\ge k_1\ge k_2\ge\cdots\ge0$. As soon as two neighbours agree, $k_{J+1}=k_J$, the sequence is periodic: $u_{2J+5}=g(u_{2J+3})$${}=g(u_{2J+1})=u_{2J+3}$, and repeating the argument shows that all later odd-indexed terms are equal; the even-indexed terms after them, $u_{2j+2}=f(u_{2j+1})$, are then equal as well. Until that happens, every step lowers $k_j$ by at least $1$, which is possible at most $710\,000\,000$ times. So some $J\le710\,000\,000$ gives

$$u_{n+2}=u_n\qquad(n\ge2J+1),$$

and $2J+1\le1\,420\,000\,001\lt10^{12}$. Call the two values of the final cycle $p\le0.71$ and $q=f(p)\ge1.001242148$. They differ, so the period is exactly $2$ and not $1$. For $n\ge2J+1$ the terms $u_n$ and $u_{n+1}$ are $p$ and $q$ in some order, hence

$$u_{10^{12}}+u_{10^{12}+1}=p+q.$$

Finiteness alone already forces a repetition among the first $M+1$ terms from $u_1$ on, where $M=1\,420\,000\,001$ is the number of grid points, ^[9] but it allows any period up to $M$, and then $u_n+u_{n+1}$ could depend on $n$; the order argument shows that the period is $2$. The actual runs are far shorter than both bounds: the odd-indexed terms fall $258$ times, from $u_1$ to $u_{517}$, the even-indexed terms rise $257$ times, from $u_2$ to $u_{516}$, and $u_{n+2}=u_n$ holds from $n=516$ on, but not for $n=515$.

The algorithm follows: apply $f$ to $u_0=-1$ any number $N\ge516$ of times and output $u_N+u_{N+1}$. A program could also stop at the first $n$ with $u_{n+2}=u_n$, which Sections 3 and 4 guarantee to exist; general cycle-finding methods such as Floyd's tortoise and hare ^[9] are unnecessary, since the period is at most $2$.

5. The smooth curve behind $f$

Without the floor we would iterate the smooth map $F$. A fixed point of $F$ is a value $x$ with $F(x)=x$; an iteration that reaches it stays there. Since $F\gt0$, every fixed point of $F$ is positive, and on $[0,\infty)$ the difference $F(x)-x$ is strictly decreasing. It is about $+0.0106$ at $x=0.85$ and about $-0.0096$ at $x=0.86$, so $F$ has exactly one fixed point $x^*$, and $0.85\lt x^*\lt0.86$. From $F'(x)=-2\ln2\cdot x\,F(x)$ we get

$$\begin{aligned}F'(x^*)&=-2\ln2\cdot(x^*)^2\\&\lt-2\ln2\cdot0.7225\lt-1,\end{aligned}$$

because $2\ln2\cdot0.7225=1.00159\ldots$; numerically $F'(x^*)=-1.014\ldots$. A fixed point with $|F'|\gt1$ is unstable: ^[3] a point at distance $\varepsilon$ from $x^*$ lands on the other side at distance about $1.014\,\varepsilon$. The slope has only just passed $-1$, as in a period-doubling bifurcation: the logistic map $x\mapsto rx(1-x)$ has fixed-point slope $-1$ at $r=3$, and for $r$ a little above $3$ its orbits settle into an alternation between two values. ^[4]

For $F$ the alternation can be proved directly. The function $tF(t)$ has the derivative $F(t)(1-2\ln2\cdot t^2)$, so it increases up to $t=0.849\ldots$ and decreases beyond. Hence $|F'(t)|=2\ln2\cdot tF(t)$ is at most $|F'(0.72)|\lt0.99$ on $J=[0.6,0.72]$, and at most $|F'(0.99)|\lt0.988$ on $[0.99,1.11]$, which contains $F(J)$. The double step $G=F\circ F$ is increasing, with $G(0.6)-0.6\approx+0.0078$ and $G(0.72)-0.72\approx-0.0015$, so $G$ maps $J$ into itself, and on $J$ the chain rule gives $0\lt G'(x)$${}=F'(F(x))F'(x)\lt0.978$. By the mean value theorem $G$ shrinks distances in $J$ by a factor below $0.978$, so it has exactly one fixed point $a$ in $J$, and every orbit of $G$ that starts in $J$ converges to $a$. The smooth sequence $v_0=-1$, $v_{n+1}=F(v_n)$ has $v_1=F(1)\approx0.71$, so its odd-indexed terms converge to $a$ and its even-indexed terms to $b=F(a)$, with $F(b)=a$. This is a 2-cycle with $a\lt x^*\lt b$, as $a\le0.72$ and $b\ge F(0.72)\gt0.99$.

Both points have the same slope, $F'(a)=-2\ln2\cdot ab$${}=F'(b)\approx-0.972$, so the multiplier of the cycle is

$$\begin{aligned}\mu&=G'(a)=F'(b)\,F'(a)\\&=(2\ln2\cdot ab)^2\approx0.945.\end{aligned}$$

Since $\mu\lt1$, the cycle attracts: ^[5] near $a$, each further pair of steps multiplies the distance to $a$ by about $\mu$. The rounded sequence follows this until rounding takes over: the ratio of successive distances of the odd-indexed terms to $a$ is about $0.958$ from $u_1$ to $u_3$ and $0.945$ from $u_{301}$ to $u_{303}$, fluctuates between about $0.90$ and $0.95$ in the last steps, and becomes $1$ once the terms stop at $p\ne a$.

6. Many cycles with one sum

We measure distances in units of $h=10^{-9}$, the spacing of the grid. For $x\in J$, both $F(x)$ and $f(x)$ lie in $[0.99,1.11]$. The inner truncation in $g(x)=f(f(x))$ lowers $F(x)$ by less than $h$, which by the slope bound of Section 5 raises the next value by less than $0.988h$, and the outer truncation lowers the result by less than $h$. Hence

$$\begin{aligned}-h\lt g(x)-F(F(x))\lt h\\(x\in J);\end{aligned}$$

a scan around $a$ finds values between $-0.97h$ and $+0.92h$. The pull towards the cycle, however, is weak:

$$F(F(x))-x\approx-(1-\mu)(x-a),$$

only $0.055$ units per unit of distance. Within about $1/(1-\mu)\approx18$ units of $a$ the rounding can cancel the pull, and $g(x)=x$ becomes possible. For nine decimals all $34$ grid points from about $17$ units below $a$ to about $16$ units above it are fixed points of $g$, and $f$ maps them in reverse order onto $34$ consecutive grid points around $b$. The rounded map therefore has $34$ different 2-cycles side by side, and they all have the same sum.

Claim. If neighbouring grid points $x$ and $x'=x+h$ in $J$ are both fixed points of $g$, then $x+f(x)=x'+f(x')$.

Proof. On $J$ we have $|F'|\lt0.99$, so by the mean value theorem $0\lt F(x)-F(x')\lt h$. Thus $10^{9}F$ drops by less than $1$, its floor drops by $0$ or $1$, and $f(x')$ is $f(x)$ or $f(x)-h$. In the first case $x'=g(x')=f(f(x'))$${}=f(f(x))=g(x)=x$, which is false. Hence $f(x')=f(x)-h$, and the two sums agree. $\square$

Applied along the band, the claim shows that all $34$ cycles share one sum, so it does not matter which of them the sequence enters. It enters the top one: if $p'$ is the largest fixed point of $g$ not exceeding $u_1$, then $g(x)\ge g(p')=p'$ for $x\ge p'$, so the odd-indexed terms never drop below $p'$, and as they stop at a fixed point, $p=p'$. This point lies about $16$ units above $a$, and $q=f(p)$ about $16.5$ units below $b$. The sum $p+q$ lies $0.54$ units below $a+b$, so rounding an accurate value of $a+b$ to nine decimals would give a wrong last digit.

7. Fewer decimals: a worked example and the length of the run

The arguments of Sections 3, 4 and 6 work just as well with $d$ decimals instead of nine, for the map $f_d(x)=10^{-d}\lfloor10^{d}F(x)\rfloor$; the problem is $d=9$. For $d=2$ the whole run can be done by hand. First $100F(1)=71.0000000139\ldots$, so $u_1=0.71$. Next $2^{-0.5041}=0.7051\ldots$ gives $100F(0.71)=100.124\ldots$, so $u_2=1.00$, and $100F(1)$ once more gives $u_3=0.71=u_1$. With two decimals $0.71$ is already a fixed point of $f_2\circ f_2$, and the sequence alternates between $0.71$ and $1.00$ from $u_1$ on.

With more decimals the run grows. The first index from which $u_{n+2}=u_n$ holds is $1$ for $d=1$ and $d=2$, and then

$28,\ 115,\ 214,\ 290,\ 383$ for $d=3,\dots,7$;
$454,\ 516,\ 607,\ 694,\ 802$ for $d=8,\dots,12$.

Section 5 explains this growth. For every $d$ the orbit stops within about $18$ grid units of the cycle, so each extra decimal requires the distance to $a$ to shrink by another factor of $10$, which at a factor $\mu$ per pair of steps costs

$$2\cdot\frac{\ln10}{\ln(1/\mu)}\approx81.5$$

steps per decimal; from $d=4$ to $d=12$ the observed average is $86$. The values scatter because the band of Section 6 does not always fill its range: for $d=3,\dots,12$ it contains between $3$ and $34$ cycles. For $d=9$, shrinking the distance from $0.71-a$ to about $16$ units, a factor of about $1.8\cdot10^{6}$, takes about $510$ steps at the rate $\mu$, close to the actual $516$.

How the Code Works

All three programs compute $u_N+u_{N+1}$ for a fixed $N$ by direct iteration in double precision. In the C++ program ^[10] the function f evaluates the formula of the statement, std::floor(std::pow(2.0, 30.403243784 - x * x)) * 1e-9, and solve(iterations) starts from u = -1.0, applies f iterations times, computes v = f(u) and returns u + v. The default $N=100\,000$ comes from the Options struct; parse_arguments accepts --iterations=N and --skip-checkpoints and exits with code $1$ on anything else or on a count below $1$, although the digit loop of parse_int_after_prefix has no overflow check for values beyond the int range. Unless skipped, run_checkpoints first tests that f(-1.0) is $0.71$ and that solve(10000) agrees with a stored reference value, both within $10^{-12}$, and exits with code $2$ on failure. main prints the sum with std::fixed and std::setprecision(9).

The Python program defines the same f with math.floor and the ** operator, runs $1000$ steps and prints with the format .9f; the Java program inlines the step with Math.floor and Math.pow, also runs $1000$ steps and prints with printf and %.9f. No program looks for the cycle: each relies on its step count exceeding the $516$ steps of Section 4, by a factor of about two for $1000$ and by far for $100\,000$. The C++ program with --iterations=516 already gives the final output, with --iterations=515 a different one.

Why double precision is safe here

The factor 1e-9 is not a binary fraction, ^[7] so u holds a double that is only close to a multiple of $10^{-9}$, and the exponent $c-x^2$, a number near $30$, carries rounding errors of the order of one unit in its last place, $2^{-48}\approx3.6\cdot10^{-15}$. ^[6] As the derivative of $2^{y}$ is $\ln2\cdot2^{y}$, a result near $10^{9}$ can be off by about $10^{9}\cdot\ln2\cdot3.6\cdot10^{-15}\approx2.5\cdot10^{-6}$. The power function adds its own error, which Java's Math.pow keeps within one unit in the last place, here $2^{-23}\approx1.2\cdot10^{-7}$; ^[8] for C++ and Python the margin was measured. A floor comes out wrong only if the exact value of $2^{c-x^2}$ lies that close to an integer. Along the actual orbit the Python evaluation differs from a 60-digit one by at most $1.8\cdot10^{-6}$, while the exact values never come closer than $0.0013$ to an integer, a safety factor of about $740$. So every floor yields the correct integer $k_n=10^{9}u_n$, and the C++ and Java recurrences select the same integers. Finally, the double sum u + v lies within $10^{-15}$ of the nine-decimal number $p+q$, far below the $5\cdot10^{-10}$ that could change its rounding, so the printed digits are exact.

Complexity and Verification

For $N$ steps the programs need $O(N)$ time and $O(1)$ memory; each step costs one power, one floor and a few arithmetic operations. The Python and Java programs evaluate $f$ $1001$ times, the C++ program $N+1=100\,001$ times in solve plus $10\,002$ times in its checkpoints, which takes well under a second. The stopping test of Section 4 would need only $518$ evaluations of $f$ and one comparison per step.

For this article the orbit was recomputed with exact integers $k_n=10^{9}u_n$ and 60-digit evaluations of $2^{c-x^2}$. This confirmed the values in Sections 1 and 2 and the periodicity from $n=516$ on. In each of the first $120\,000$ steps the floors in the recurrence code of all three programs select the exact integers, and the outputs of the three programs, and of the C++ program with $N=516$ and $N=10^{6}$, equal the exact value of $p+q$, while $N=515$ gives a different sum. Further high-precision computations checked the bounds of Section 5, the band of $34$ fixed points of $g$ and its image, the common sum of its cycles (and of the band cycles for every $d$ from $3$ to $12$), $\mu$ by numerical differentiation of $G$, and the run lengths of Section 7. The C++ checkpoints test $f(-1)$ and one long run, but not periodicity.

Footnotes and References

The references below supply the background facts used above. The order argument, the bound on the length of the runs, the band of cycles with a common sum and all numerical values are derived or computed in this article.

Project Euler 197 — A Recursively Defined Sequence. The official statement defines f(x) = ⌊2^(30.403243784 − x²)⌋ × 10⁻⁹ and the sequence u₀ = −1, uₙ₊₁ = f(uₙ), and asks for uₙ + uₙ₊₁ with n = 10¹², given with nine digits after the decimal point.
Wikipedia — Floor and ceiling functions. The article gives the inequalities x − 1 < ⌊x⌋ ≤ x and states that the floor function is monotonically non-decreasing, the two facts used in Sections 1 and 3.
Wikipedia — Stability theory. Its section on fixed points of maps states that a fixed point a of a continuously differentiable map of the real line is stable if |f′(a)| < 1 and unstable if |f′(a)| > 1, as used for x* in Section 5.
Wikipedia — Logistic map. At r = 3 the derivative at the fixed point of the logistic map reaches −1 and a period-doubling bifurcation occurs; for r between 3 and 1 + √6 the orbit approaches a permanent oscillation between two values.
Wikipedia — Periodic point. A hyperbolic periodic point of period n of a differentiable map is attracting when the derivative of the n-th iterate has absolute value below 1 and repelling when it is above 1; Section 5 applies this with n = 2 to the cycle {a, b}.
Wikipedia — Double-precision floating-point format. A double has a 53-bit significand, good for 15 to 17 significant decimal digits, and for numbers in the normal range rounding to the nearest double has a relative error of at most 2⁻⁵³.
Python tutorial — Floating-Point Arithmetic: Issues and Limitations. Most decimal fractions cannot be represented exactly as binary fractions, and Python floats are IEEE 754 binary64 values with 53 bits of precision, so 1e-9 and the terms of the sequence are stored only approximately.
Java SE 21 — java.lang.Math. Math.pow must be within 1 ulp of the exact result, Math.floor returns the largest double not above its argument that is a mathematical integer, and the Math methods need not give bit-for-bit identical results on all platforms.
Wikipedia — Cycle detection. Iterating a function from a finite set to itself must eventually repeat a value, after which the sequence is periodic; Floyd's tortoise-and-hare algorithm and Brent's algorithm find the cycle while storing only a constant number of values.
C++ — Euler197.cpp. The linked immutable revision contains f, solve, run_checkpoints and the command-line options described above. It is the implementation record for this article.

Problem 197 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>

namespace {

struct Options {
    int iterations = 100000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--iterations=", options.iterations)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.iterations >= 1;
}

double f(const double x) {
    return std::floor(std::pow(2.0, 30.403243784 - x * x)) * 1e-9;
}

double solve(const int iterations) {
    double u = -1.0;
    for (int i = 0; i < iterations; ++i) {
        u = f(u);
    }
    const double v = f(u);
    return u + v;
}

bool run_checkpoints() {
    if (std::abs(f(-1.0) - 0.7100000000) > 1e-12) {
        std::cerr << "Checkpoint failed for first sequence value" << '\n';
        return false;
    }
    if (std::abs(solve(10000) - 1.710637717) > 1e-12) {
        std::cerr << "Checkpoint failed for converged two-cycle sum" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(9) << solve(options.iterations) << '\n';
    return 0;
}

Python

# Problem 197: Investigating the behaviour of a recursively defined sequence.
# f(x) = floor(2^(30.403243784-x^2)) * 10^-9
# u_n+1 = f(u_n), u_0 = -1. Find u_n + u_(n+1) for large n.

import math

def solve():
    def f(x):
        return math.floor(2**(30.403243784 - x*x)) * 1e-9
    u = -1.0
    for _ in range(1000):
        u = f(u)
    v = f(u)
    print(f"{u+v:.9f}")

solve()

Java

public class Euler197 {
    public static void main(String[] args) {
        double u = -1.0;
        for (int i = 0; i < 1000; i++)
            u = Math.floor(Math.pow(2, 30.403243784 - u * u)) * 1e-9;
        double v = Math.floor(Math.pow(2, 30.403243784 - u * u)) * 1e-9;
        System.out.printf("%.9f%n", u + v);
    }
}