Detailed mathematical approach

Problem Summary

Problem 825 asks for the partial sum

$T(N)=\sum_{n=2}^{N} s_n,\qquad N=10^{14},$

where the rational sequence \(s_n\) comes from the chasing game and is generated by two linear recurrences. A direct summation over \(10^{14}-1\) terms is impossible, so the solution extracts the exact dominant asymptotic term, sums that term in closed form, and keeps only a short rapidly convergent correction.

Mathematical Approach

Write

$s_n=\frac{a_n}{b_n},$

with initial values

$a_2=7,\quad a_3=35,\quad a_4=121,$

$b_2=11,\quad b_3=73,\quad b_4=395,\quad b_5=1933,$

and recurrences

$a_n=3a_{n-1}+3a_{n-2}-a_{n-3}\qquad (n\ge 5),$

$b_n=8b_{n-1}-18b_{n-2}+8b_{n-3}-b_{n-4}\qquad (n\ge 6).$

The implementations are built around the structure hidden in these recurrences.

Step 1: Factor the characteristic polynomials

The numerator recurrence has characteristic polynomial

$r^3-3r^2-3r+1=(r+1)(r^2-4r+1).$

The denominator recurrence has characteristic polynomial

$r^4-8r^3+18r^2-8r+1=(r^2-4r+1)^2.$

So the key roots are

$\lambda=2+\sqrt3,\qquad \mu=2-\sqrt3=\lambda^{-1}.$

The repeated factor in the denominator is the reason an extra linear factor in \(n\) appears there, and that is exactly what turns the ratio \(a_n/b_n\) into a shifted harmonic term.

Step 2: Solve the recurrences explicitly

Using the initial values, the closed forms become

$A=\frac{3-\sqrt3}{2},\qquad B=\frac{3+\sqrt3}{2},$

$a_n=A\lambda^n+B\mu^n-2(-1)^n,$

$b_n=\left(An-\frac12\right)\lambda^n+\left(Bn-\frac12\right)\mu^n.$

These formulas satisfy both the recurrences and the listed starting values. They also make the asymptotic behavior transparent: \(\lambda>1\) dominates, while \(\mu<1\) contributes only exponentially tiny corrections.

Step 3: Extract the shifted harmonic main term

Divide numerator and denominator by \(\lambda^n\):

$s_n=\frac{A+B\mu^{2n}-2(-1)^n\lambda^{-n}}{An-\frac12+\left(Bn-\frac12\right)\mu^{2n}}.$

Since \(\mu^{2n}=\lambda^{-2n}\), all non-dominant pieces decay exponentially. Therefore

$s_n=\frac{A+O(\lambda^{-n})}{An-\frac12+O(n\lambda^{-2n})}.$

Because

$\frac{1}{2A}=\frac{3+\sqrt3}{6},$

it is natural to define

$c=\frac{3+\sqrt3}{6}.$

Then

$An-\frac12=A(n-c),$

and hence

$s_n=\frac{1}{n-c}+O(\lambda^{-n}).$

The crucial point is that the correction is exponentially small, not merely polynomially small.

Step 4: Split the huge sum into a closed form and a short correction

Now decompose

$T(N)=\sum_{n=2}^{N}\frac{1}{n-c}+\sum_{n=2}^{N}\left(s_n-\frac{1}{n-c}\right).$

The second sum converges very quickly because the summand is exponentially small. If we truncate it at a modest cutoff \(K\), we get

$E_K=\sum_{n=2}^{K}\left(s_n-\frac{1}{n-c}\right),$

and the omitted tail is negligible for the requested \(8\)-decimal output. The implementations choose

$K=60.$

Step 5: Sum the main term with the digamma identity

The digamma function satisfies

$\psi(x+1)-\psi(x)=\frac1x.$

Applying this telescoping identity to the shifted harmonic part gives

$\sum_{n=2}^{N}\frac{1}{n-c}=\psi(N+1-c)-\psi(2-c).$

So the desired value is

$T(N)=\psi(N+1-c)-\psi(2-c)+\sum_{n=2}^{N}\left(s_n-\frac{1}{n-c}\right),$

and numerically it is enough to use

$T(N)\approx \psi(N+1-c)-\psi(2-c)+E_{60}.$

Worked Example: the checkpoint \(T(10)\)

The first few terms are

$s_2=\frac{7}{11}\approx 0.63636364,\qquad s_3=\frac{35}{73}\approx 0.47945205,\qquad s_4=\frac{121}{395}\approx 0.30632911.$

With

$c=\frac{3+\sqrt3}{6}\approx 0.78867513,$

the shifted harmonic part up to \(10\) is

$\sum_{n=2}^{10}\frac{1}{n-c}\approx 2.54851473.$

The correction over the same range is

$\sum_{n=2}^{10}\left(s_n-\frac{1}{n-c}\right)\approx -0.16616191.$

Therefore

$T(10)\approx 2.54851473-0.16616191=2.38235282,$

which matches the numerical checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations build the rational terms only up to \(n=60\). They generate the numerator and denominator sequences from the recurrences using exact integer arithmetic first, and only convert to floating point when forming each ratio \(s_n\) inside the correction sum. This keeps the short precomputation stable and accurate.

After the correction \(E_{60}\) is accumulated, the remaining task is to evaluate two digamma values: one at \(N+1-c\) and one at \(2-c\). The implementation does this numerically in two stages. First, it repeatedly uses the recurrence \(\psi(x)=\psi(x+1)-1/x\) until the argument is safely large. Then it applies the asymptotic expansion

$\psi(x)=\log x-\frac{1}{2x}-\frac{1}{12x^2}+\frac{1}{120x^4}-\frac{1}{252x^6}+\frac{1}{240x^8}-\frac{1}{132x^{10}}+\frac{691}{32760x^{12}}+\cdots.$

Finally the program combines

$\psi(N+1-c)-\psi(2-c)+E_{60}$

and prints the result to \(8\) decimal places.

Complexity Analysis

If the cutoff is denoted by \(K\), generating all needed recurrence values and summing the correction costs \(O(K)\) time and \(O(K)\) memory. The two digamma evaluations use only constant-time arithmetic and a fixed number of asymptotic-series terms. Since the implementation fixes \(K=60\), the whole computation is effectively constant-time even though \(N=10^{14}\) is enormous.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=825
2. Digamma function: Wikipedia — Digamma function
3. Linear recurrence with constant coefficients: Wikipedia — Linear recurrence with constant coefficients
4. Asymptotic expansion: Wikipedia — Asymptotic expansion
5. Bernoulli numbers in special-function expansions: Wikipedia — Bernoulli number