Problem 825: Chasing Game
View on Project EulerProject Euler Problem 825 Solution
EulerSolve provides an optimized solution for Project Euler Problem 825, Chasing Game, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: https://projecteuler.net/problem=825
- Digamma function: Wikipedia — Digamma function
- Linear recurrence with constant coefficients: Wikipedia — Linear recurrence with constant coefficients
- Asymptotic expansion: Wikipedia — Asymptotic expansion
- Bernoulli numbers in special-function expansions: Wikipedia — Bernoulli number
Problem 825 source code
C++
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>
using int64 = long long;
using i128 = __int128_t;
using ld = long double;
namespace {
ld digamma_ld(ld x) {
if (!(x > 0.0L) || !std::isfinite(static_cast<double>(x))) {
return std::numeric_limits<ld>::quiet_NaN();
}
ld res = 0.0L;
while (x < 10.0L) {
res -= 1.0L / x;
x += 1.0L;
}
const ld inv = 1.0L / x;
const ld inv2 = inv * inv;
ld series = inv2 * (-1.0L / 12.0L
+ inv2 * (1.0L / 120.0L
+ inv2 * (-1.0L / 252.0L
+ inv2 * (1.0L / 240.0L
+ inv2 * (-1.0L / 132.0L
+ inv2 * (691.0L / 32760.0L))))));
res += std::logl(x) - 0.5L * inv + series;
return res;
}
inline ld i128_to_ld(i128 v) {
return static_cast<ld>(v);
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
const unsigned long long N = 100000000000000ULL;
const ld sqrt3 = std::sqrtl(3.0L);
const ld c = (3.0L + sqrt3) / 6.0L;
const int K = 60;
std::vector<i128> num(K + 1, 0);
std::vector<i128> den(K + 1, 0);
num[2] = 7;
num[3] = 35;
num[4] = 121;
den[2] = 11;
den[3] = 73;
den[4] = 395;
den[5] = 1933;
for (int n = 5; n <= K; ++n) {
num[n] = 3 * num[n - 1] + 3 * num[n - 2] - num[n - 3];
}
for (int n = 6; n <= K; ++n) {
den[n] = 8 * den[n - 1] - 18 * den[n - 2] + 8 * den[n - 3] - den[n - 4];
}
ld E = 0.0L;
for (int n = 2; n <= K; ++n) {
const ld s = i128_to_ld(num[n]) / i128_to_ld(den[n]);
E += s - 1.0L / (static_cast<ld>(n) - c);
}
#ifdef VALIDATE
{
const ld s2 = i128_to_ld(num[2]) / i128_to_ld(den[2]);
const ld expected_s2 = 7.0L / 11.0L;
if (std::fabsl(s2 - expected_s2) > 1e-18L) {
throw std::runtime_error("Validation failed: S(2) mismatch");
}
ld t10 = 0.0L;
for (int n = 2; n <= 10; ++n) {
t10 += i128_to_ld(num[n]) / i128_to_ld(den[n]);
}
const ld expected_t10 = 2.38235282L;
if (std::fabsl(t10 - expected_t10) > 5e-8L) {
throw std::runtime_error("Validation failed: T(10) mismatch");
}
const int ncheck = 50;
const ld s = i128_to_ld(num[ncheck]) / i128_to_ld(den[ncheck]);
const ld delta = 1.0L / s - static_cast<ld>(ncheck);
if (std::fabsl(delta + c) > 1e-12L) {
throw std::runtime_error("Validation failed: asymptotic shift mismatch");
}
}
#endif
const ld x_big = static_cast<ld>(N) + 1.0L - c;
const ld x_small = 2.0L - c;
const ld ans = digamma_ld(x_big) - digamma_ld(x_small) + E;
std::cout.setf(std::ios::fixed);
std::cout << std::setprecision(8) << ans << "\n";
return 0;
}
Python
import math
def digamma(x):
if x <= 0.0 or not math.isfinite(x):
return float('nan')
res = 0.0
while x < 10.0:
res -= 1.0 / x
x += 1.0
inv = 1.0 / x
inv2 = inv * inv
series = inv2 * (-1.0 / 12.0
+ inv2 * (1.0 / 120.0
+ inv2 * (-1.0 / 252.0
+ inv2 * (1.0 / 240.0
+ inv2 * (-1.0 / 132.0
+ inv2 * (691.0 / 32760.0))))))
res += math.log(x) - 0.5 * inv + series
return res
def solve():
N = 100000000000000
sqrt3 = math.sqrt(3.0)
c = (3.0 + sqrt3) / 6.0
K = 60
num = [0] * (K + 1)
den = [0] * (K + 1)
num[2] = 7
num[3] = 35
num[4] = 121
den[2] = 11
den[3] = 73
den[4] = 395
den[5] = 1933
for n in range(5, K + 1):
num[n] = 3 * num[n - 1] + 3 * num[n - 2] - num[n - 3]
for n in range(6, K + 1):
den[n] = 8 * den[n - 1] - 18 * den[n - 2] + 8 * den[n - 3] - den[n - 4]
E = 0.0
for n in range(2, K + 1):
s = float(num[n]) / float(den[n])
E += s - 1.0 / (float(n) - c)
x_big = float(N) + 1.0 - c
x_small = 2.0 - c
ans = digamma(x_big) - digamma(x_small) + E
return "{:.8f}".format(ans)
if __name__ == "__main__":
print(solve())
Java
import java.math.BigInteger;
public class Euler825 {
static double digamma(double x) {
if (x <= 0.0 || Double.isInfinite(x) || Double.isNaN(x)) {
return Double.NaN;
}
double res = 0.0;
while (x < 10.0) {
res -= 1.0 / x;
x += 1.0;
}
double inv = 1.0 / x;
double inv2 = inv * inv;
double series = inv2 * (-1.0 / 12.0
+ inv2 * (1.0 / 120.0
+ inv2 * (-1.0 / 252.0
+ inv2 * (1.0 / 240.0
+ inv2 * (-1.0 / 132.0
+ inv2 * (691.0 / 32760.0))))));
res += Math.log(x) - 0.5 * inv + series;
return res;
}
public static String solve() {
long N = 100000000000000L;
double sqrt3 = Math.sqrt(3.0);
double c = (3.0 + sqrt3) / 6.0;
int K = 60;
BigInteger[] num = new BigInteger[K + 1];
BigInteger[] den = new BigInteger[K + 1];
for (int i = 0; i <= K; i++) {
num[i] = BigInteger.ZERO;
den[i] = BigInteger.ZERO;
}
num[2] = BigInteger.valueOf(7);
num[3] = BigInteger.valueOf(35);
num[4] = BigInteger.valueOf(121);
den[2] = BigInteger.valueOf(11);
den[3] = BigInteger.valueOf(73);
den[4] = BigInteger.valueOf(395);
den[5] = BigInteger.valueOf(1933);
for (int n = 5; n <= K; ++n) {
num[n] = num[n - 1].multiply(BigInteger.valueOf(3))
.add(num[n - 2].multiply(BigInteger.valueOf(3)))
.subtract(num[n - 3]);
}
for (int n = 6; n <= K; ++n) {
den[n] = den[n - 1].multiply(BigInteger.valueOf(8))
.subtract(den[n - 2].multiply(BigInteger.valueOf(18)))
.add(den[n - 3].multiply(BigInteger.valueOf(8)))
.subtract(den[n - 4]);
}
double E = 0.0;
for (int n = 2; n <= K; ++n) {
// Convert to double. Alternatively we can use BigDecimal
double nDouble = num[n].doubleValue();
double dDouble = den[n].doubleValue();
double s = nDouble / dDouble;
E += s - 1.0 / (n - c);
}
double xBig = (double) N + 1.0 - c;
double xSmall = 2.0 - c;
double ans = digamma(xBig) - digamma(xSmall) + E;
return String.format(java.util.Locale.US, "%.8f", ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}