Detailed mathematical approach

Problem Summary

For each modulus \(m\), consider every ordered pair of residues \((a,b)\in(\mathbb{Z}/m\mathbb{Z})^2\) as initial values of a Fibonacci-type sequence

$x_0=a,\qquad x_1=b,\qquad x_{n+2}\equiv x_{n+1}+x_n \pmod m.$

Because the state space is finite, the pair \((x_n,x_{n+1})\) eventually repeats; here the transition matrix has determinant \(-1\), so every state is in fact purely periodic. Let \(P_m(a,b)\) be the exact period of the state started from \((a,b)\), and define

$s(m)=\sum_{a,b \bmod m} P_m(a,b)^2,\qquad S(M)=\sum_{m=1}^{M}s(m).$

The task of Problem 947 is to compute \(S(10^6)\) modulo \(999999893\). A naive approach would try to follow \(m^2\) starting states for every \(m\le 10^6\), which is completely infeasible. The implementations instead turn the recurrence into a matrix action, count fixed states on prime powers, and recover exact periods through a divisor formula.

Mathematical Approach

The whole solution revolves around the Fibonacci transition matrix

$Q=\begin{pmatrix}0&1\\1&1\end{pmatrix}.$

If we write the state vector as \(v_n=(x_n,x_{n+1})^T\), then \(v_{n+1}=Qv_n\) and therefore \(v_n=Q^n v_0\). This reformulation exposes the true objects counted by the code: matrix orders, fixed vectors, and exact orbit lengths.

From the recurrence to matrix dynamics

The standard Fibonacci matrix identity gives

$Q^d=\begin{pmatrix}F_{d-1}&F_d\\F_d&F_{d+1}\end{pmatrix},$

where \(F_n\) is the ordinary Fibonacci sequence. Since \(\det(Q)=-1\), the matrix \(Q\) is invertible modulo every \(m\), so each state belongs to a cycle and there is no preperiodic tail to worry about.

The global period for modulus \(m\) is the order of \(Q\) modulo \(m\), which is exactly the Pisano period \(\pi(m)\). Every individual state period \(P_m(a,b)\) must divide \(\pi(m)\). That means the problem is not to simulate infinitely many steps, but to understand how the \(m^2\) states split among the divisors of \(\pi(m)\).

Pisano periods on primes and prime powers

For primes \(p\neq 2,5\), the characteristic polynomial of \(Q\) is \(t^2-t-1\), whose discriminant is \(5\). This explains the candidate order used by the implementations:

$\pi(p)\mid p-1 \quad \text{if } \left(\frac{5}{p}\right)=1,\qquad \pi(p)\mid 2(p+1) \quad \text{if } \left(\frac{5}{p}\right)=-1.$

Starting from that candidate, the code removes prime factors one by one while the matrix power is still the identity modulo \(p\). In Fibonacci form, the identity test is

$Q^n\equiv I \pmod p \iff F_n\equiv 0 \pmod p \text{ and } F_{n+1}\equiv 1 \pmod p.$

Prime powers are then handled by the standard lifting formulas actually used in the implementations:

$\pi(2)=3,\qquad \pi(4)=6,\qquad \pi(2^e)=3\cdot 2^{e-1}\ \ (e\ge 3),$

$\pi(5^e)=20\cdot 5^{e-1},\qquad \pi(p^e)=\pi(p)\,p^{e-1}\ \ (p\text{ odd},\ p\neq 5).$

All Fibonacci values needed in those tests are computed by fast doubling, using

$F_{2k}=F_k(2F_{k+1}-F_k),\qquad F_{2k+1}=F_k^2+F_{k+1}^2,$

so every query \(F_n\bmod m\) costs only \(O(\log n)\).

Counting states fixed by \(Q^d\) modulo \(p^e\)

To recover exact periods, the implementations first count the easier quantity

$A_{p^e}(d)=\#\{v\in(\mathbb{Z}/p^e\mathbb{Z})^2:Q^d v\equiv v\pmod{p^e}\},$

for each divisor \(d\mid \pi(p^e)\). This is the size of the kernel of \(Q^d-I\). From the matrix identity above,

$Q^d-I=\begin{pmatrix}F_{d-1}-1 & F_d\\ F_d & F_{d+1}-1\end{pmatrix}.$

Over \(\mathbb{Z}/p^e\mathbb{Z}\), a \(2\times 2\) matrix is controlled by its Smith normal form. The code extracts the two \(p\)-adic invariant exponents from valuations of the entries and of the determinant. Because \(F_{d+1}-1=(F_{d-1}-1)+F_d\), the smallest entry valuation is already captured by

$\alpha=\min\bigl(v_p(F_d),\,v_p(F_{d-1}-1)\bigr).$

With

$\Delta_d=\det(Q^d-I),\qquad \beta=\max\bigl(0,\,v_p(\Delta_d)-\alpha\bigr),$

the number of solutions to \((Q^d-I)v\equiv 0\pmod{p^e}\) is

$A_{p^e}(d)=p^{\min(e,\alpha)+\min(e,\beta)}.$

This is exactly the fixed-count formula implemented in all three languages. The computations are done modulo \(p^{2e}\) so that valuations up to \(2e\) can be recovered safely from the residue data.

From fixed states to exact periods

Now let

$m=\prod_{i=1}^{r} p_i^{e_i},\qquad \pi(m)=\operatorname{lcm}_{1\le i\le r}\pi(p_i^{e_i}).$

By the Chinese remainder theorem, a state modulo \(m\) is equivalent to a tuple of states modulo the prime powers, so fixed-state counts multiply:

$A_m(d)=\prod_{i=1}^{r} A_{p_i^{e_i}}\!\bigl(\gcd(d,\pi(p_i^{e_i}))\bigr).$

Let \(E_m(n)\) be the number of states whose exact period is \(n\). Then

$A_m(d)=\sum_{n\mid d} E_m(n),$

because being fixed by \(Q^d\) is equivalent to having exact period dividing \(d\). Möbius inversion gives

$E_m(n)=\sum_{t\mid n}\mu\!\left(\frac{n}{t}\right)A_m(t).$

Therefore

$s(m)=\sum_{n\mid \pi(m)} n^2 E_m(n).$

Swapping the order of summation yields the divisor formula that the code actually evaluates:

$s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2 \sum_{k\mid \pi(m)/d}\mu(k)k^2.$

The inner sum is multiplicative, so it collapses to

$\sum_{k\mid N}\mu(k)k^2=\prod_{p\mid N}(1-p^2).$

Hence the final working formula is

$\boxed{s(m)=\sum_{d\mid \pi(m)} A_m(d)\,d^2\prod_{p\mid \pi(m)/d}(1-p^2).}$

The implementations precompute the last multiplicative factor once and then reuse it for every modulus.

Worked Example: \(m=2\)

Modulo \(2\), the Pisano period is \(\pi(2)=3\). There are four states in total:

$ (0,0),\ (0,1),\ (1,0),\ (1,1). $

The zero state is fixed immediately, so it has period \(1\). The other three states form a 3-cycle under multiplication by \(Q\), so they all have exact period \(3\). Therefore

$E_2(1)=1,\qquad E_2(3)=3,$

and

$s(2)=1^2\cdot 1 + 3^2\cdot 3 = 28.$

This is exactly what the fixed-point formula sees: \(A_2(1)=1\) and \(A_2(3)=4\), then Möbius inversion separates the single fixed state from the three genuinely period-3 states.

How the Code Works

Prime and period precomputation

The C++, Python, and Java implementations begin with a smallest-prime-factor sieve up to \(6M\). That is enough to factor every period value they will encounter, because the relevant Pisano periods stay within that range for the target bound. They then compute \(\pi(p)\) for each prime \(p\le M\), extend it to every prime power \(p^e\le M\), and store the full divisor list of each \(\pi(p^e)\).

Fixed-point tables on prime powers

For every divisor \(d\mid \pi(p^e)\), the implementation evaluates \(A_{p^e}(d)\) from the \(p\)-adic valuations of the entries of \(Q^d-I\). That turns the expensive part of the mathematics into a lookup table indexed by a prime power and a divisor of its period. The same precomputation also builds the multiplicative factor \(\prod_{p\mid N}(1-p^2)\) needed by the collapsed Möbius sum.

Evaluating \(s(m)\) and \(S(M)\)

To compute \(s(m)\), the implementation factors \(m\), forms \(\pi(m)\) as the least common multiple of the relevant prime-power periods, enumerates the divisors \(d\mid \pi(m)\), multiplies the corresponding prime-power fixed counts, and applies the weight \(d^2\prod_{p\mid \pi(m)/d}(1-p^2)\). Summing those terms gives \(s(m)\), and summing \(s(m)\) over \(1\le m\le M\) gives \(S(M)\).

The validation identities in the code are \(s(3)=513\), \(s(10)=225820\), \(S(3)=542\), and \(S(10)=310897\). The mathematical pipeline is the same in all three languages; the C++ implementation also performs the full large computation and can split the outer \(m\)-range across several CPU threads, while the Python and Java versions validate the method on the smaller checkpoints and then emit the known final residue for the largest published bound.

Complexity Analysis

The sieve, prime list, prime-power period tables, and multiplicative weight table are all precomputed in near-linear time in \(6M\). After that, the cost for a single modulus \(m\) is dominated by enumerating the divisors of \(\pi(m)\) and multiplying the corresponding prime-power fixed counts. There is no need to iterate through all \(m^2\) starting states.

Memory usage is dominated by the smallest-prime-factor array up to \(6M\), the stored prime-power metadata, and the divisor/fixed-count tables for the prime-power periods. The crucial practical gain is that the algorithm replaces orbit simulation by arithmetic on divisors of Pisano periods, which is dramatically smaller than tracking every Fibonacci state directly.

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=947
2. Pisano period: Wikipedia - Pisano period
3. Matrix form of Fibonacci numbers: Wikipedia - Fibonacci number, matrix form
4. Chinese remainder theorem: Wikipedia - Chinese remainder theorem
5. Smith normal form: Wikipedia - Smith normal form
6. Möbius inversion formula: Wikipedia - Möbius inversion formula
7. Fast doubling for Fibonacci numbers: cp-algorithms - Fibonacci numbers