Detailed mathematical approach

Problem Summary

The quantity to compute is

\[ \sum_{n=1}^{L} F_n G(n) \pmod{10^9+9}, \qquad L=10^{14}, \]

where \(F_n\) is the Fibonacci sequence and \(G(n)\) counts the residue classes \(x \pmod n\) satisfying

\[ x^2 \equiv x+1 \pmod n. \]

The published checkpoint is

\[ \sum_{n=1}^{10^3} F_n G(n)\equiv 190950976 \pmod{10^9+9}. \]

A direct sweep up to \(10^{14}\) is hopeless, so the solution rewrites \(G(n)\) in arithmetic terms, turns the Fibonacci weight into exponential weights via Binet's formula, and then evaluates the resulting sums with Möbius inversion and a sliding-window count for a binary quadratic form.

Mathematical Approach

Prime-power behavior of \(G(n)\)

The congruence is the same as

\[ x^2-x-1 \equiv 0 \pmod n, \]

whose discriminant is \(5\). That immediately explains the local structure.

For \(2^e\) there are no solutions, so \(G(2^e)=0\). For \(5\) there is exactly one solution, namely \(x\equiv 3 \pmod 5\), but it does not lift to \(25\), so

\[ G(5)=1, \qquad G(5^e)=0 \ \text{ for } e\ge 2. \]

For an odd prime \(p\neq 5\), the polynomial has roots modulo \(p\) exactly when \(5\) is a quadratic residue modulo \(p\). By quadratic reciprocity this happens precisely for

\[ p\equiv 1,4 \pmod 5. \]

In that case there are two distinct roots modulo \(p\), the derivative \(2x-1\) is nonzero at each root, and Hensel lifting gives two roots modulo every \(p^e\). If \(p\equiv 2,3 \pmod 5\), there are no roots at any power of \(p\). Hence

\[ G(p^e)= \begin{cases} 2, & p\equiv 1,4 \pmod 5,\\ 0, & p\equiv 2,3 \pmod 5, \end{cases} \qquad (p\neq 5,\ p \text{ odd}). \]

By the Chinese remainder theorem, \(G\) is multiplicative. Therefore, if

\[ n=5^\varepsilon \prod_{i=1}^{k} p_i^{a_i}\prod_{j=1}^{m} q_j^{b_j}, \]

where every \(p_i\equiv 1,4 \pmod 5\) and every \(q_j\) is either \(2\) or congruent to \(2\) or \(3\) modulo \(5\), then \(G(n)=0\) unless \(m=0\) and \(\varepsilon\in\{0,1\}\). In the surviving case,

\[ G(n)=2^k. \]

From roots to the norm form \(Q(a,b)=a^2-ab-b^2\)

The next step moves to the quadratic ring generated by a root of \(t^2-t-1\). Let \(\varphi\) and \(\psi=1-\varphi\) be the two roots of \(t^2-t-1=0\), so \(\varphi^2=\varphi+1\). In the ring \(\mathbb Z[\varphi]\), the norm of \(a-b\varphi\) is

\[ N(a-b\varphi)=(a-b\varphi)(a-b\psi)=a^2-ab-b^2. \]

This binary quadratic form is exactly the arithmetic object used by the solver:

\[ Q(a,b)=a^2-ab-b^2. \]

The prime factors that contribute to \(G(n)\) are precisely the primes that split in \(\mathbb Z[\varphi]\). Choosing one root of \(x^2-x-1\) at each split prime power is equivalent to choosing one prime factor above each rational prime. Multiplying those local choices produces an algebraic integer of norm \(n\).

Two issues remain: units and non-primitive representations. Multiplying by a unit does not change the norm, so one imposes a reduction region to select a single representative. The implementations use the classical reduced region

\[ a\ge 2b>0. \]

Also, if \(\gcd(a,b)>1\), then the representation is not primitive and corresponds to extra square factors. After enforcing both conditions, one obtains the key identity

\[ G(n)=\#\{(a,b):\ a\ge 2b>0,\ \gcd(a,b)=1,\ Q(a,b)=n\}. \]

Fibonacci weights become exponential weights

The modulus \(10^9+9\) is prime, \(5\) has a square root modulo this prime, and the two modular roots \(\varphi,\psi\) of \(t^2-t-1\) satisfy

\[ \varphi\psi=-1. \]

So Binet's formula is valid in the finite field:

\[ F_n=\frac{\varphi^n-\psi^n}{\sqrt 5}. \]

That converts the original sum into two weighted counts of quadratic-form values. Define

\[ P_w(L)= \sum_{\substack{a\ge 2b>0\\ \gcd(a,b)=1\\ Q(a,b)\le L}} w^{Q(a,b)}. \]

Then the desired answer is

\[ \sum_{n\le L} F_n G(n) = \frac{P_\varphi(L)-P_\psi(L)}{\sqrt 5} \pmod{10^9+9}. \]

Removing the coprimality condition with Möbius inversion

The form \(Q\) is homogeneous of degree two:

\[ Q(ga,gb)=g^2Q(a,b). \]

Let

\[ A_w(L)= \sum_{\substack{a\ge 2b>0\\ Q(a,b)\le L}} w^{Q(a,b)} \]

be the same sum without \(\gcd(a,b)=1\). Möbius inversion then gives

\[ P_w(L)= \sum_{g\le \sqrt L} \mu(g)\, A_{w^{g^2}}\!\left(\left\lfloor \frac{L}{g^2}\right\rfloor\right). \]

So the primitive count is recovered by inclusion-exclusion over the common divisor \(g\).

Diagonalizing the form and splitting by parity

The decisive algebraic identity is

\[ 4Q(a,b)=(2a-b)^2-5b^2. \]

Set

\[ u=2a-b,\qquad v=b. \]

Because \(a\ge 2b>0\), we have \(u\ge 3v>0\). Also \(u\equiv v \pmod 2\), so only two parity branches occur.

If \(u=2m\) and \(v=2t\), then

\[ Q(a,b)=m^2-5t^2, \qquad m\ge 3t, \qquad t\ge 1. \]

If \(u=2m+1\) and \(v=2t+1\), then

\[ Q(a,b)=m(m+1)-5t(t+1)-1, \qquad m\ge 3t+1, \qquad t\ge 0. \]

Therefore \(A_w(L)\) is the sum of two one-dimensional window sums:

\[ \sum_{t\ge 1} \sum_{m=3t}^{\lfloor \sqrt{L+5t^2}\rfloor} w^{m^2-5t^2}, \]

\[ \sum_{t\ge 0} \sum_{m=3t+1}^{\left\lfloor(\sqrt{4L+20t^2+20t+5}-1)/2\right\rfloor} w^{m(m+1)-5t(t+1)-1}. \]

Worked example: \(n=11\)

Since \(11\equiv 1 \pmod 5\), there are two solutions of \(x^2\equiv x+1 \pmod{11}\), namely \(x\equiv 4\) and \(x\equiv 8\). So \(G(11)=2\).

The reduced primitive representations of \(11\) by \(Q(a,b)\) are

\[ Q(4,1)=16-4-1=11, \qquad Q(5,2)=25-10-4=11. \]

The first pair lands in the odd branch: \((u,v)=(2\cdot 4-1,1)=(7,1)\), so \(u=2m+1\), \(v=2t+1\) with \((m,t)=(3,0)\), and indeed

\[ m(m+1)-5t(t+1)-1=3\cdot 4-0-1=11. \]

The second pair lands in the even branch: \((u,v)=(8,2)\), so \((m,t)=(4,1)\), and

\[ m^2-5t^2=16-5=11. \]

This small case shows why the pair count agrees with \(G(n)\) and why both parity branches are needed in the final summation.

How the Code Works

Arithmetic precomputation and verification

The C++, Python, and Java implementations begin by working modulo \(10^9+9\). They verify that the chosen modular constants really satisfy \(\sqrt 5^2=5\), \(\varphi^2=\varphi+1\), \(\psi^2=\psi+1\), and \(\varphi\psi=-1\). They also compare three views of the problem on small inputs: direct brute force for the congruence, the prime-factor formula for \(G(n)\), and the reduced-pair count for the quadratic form. Finally, they check the published sample sum for \(L=10^3\).

Evaluating one non-primitive weighted sum

For a fixed weight \(w\), the implementation computes \(A_w(L)\) by scanning the two parity branches separately. It does not recompute \(w^{m^2}\) or \(w^{m(m+1)}\) from scratch. Instead it updates these powers multiplicatively, using the identities

\[ w^{(m+1)^2}=w^{m^2}w^{2m+1}, \qquad w^{(m+1)(m+2)}=w^{m(m+1)}w^{2m+2}. \]

As \(t\) grows, the allowed interval of \(m\) moves to the right. The implementation keeps a running window sum: new terms are appended when the upper bound increases, and expired terms are removed when the lower bound advances. Each admissible power enters once and leaves once.

Möbius sweep and final combination

After that, the implementation iterates over \(g\le \sqrt L\) with \(\mu(g)\neq 0\), evaluates the previous routine at the scaled limit \(\lfloor L/g^2\rfloor\), and combines the results with the Möbius sign. It precomputes the sequences \(\varphi^{g^2}\) and \(\varphi^{-g^2}\), which are enough for both branches because

\[ \psi=-\varphi^{-1}. \]

So \(\psi^{g^2}\) differs from \(\varphi^{-g^2}\) only by the parity-dependent sign of \(g^2\). The final step subtracts the \(\psi\)-sum from the \(\varphi\)-sum and multiplies by \(1/\sqrt 5\). The C++ and Java implementations can split the Möbius range across threads, while the Python implementation supports the same mathematics either serially or across processes.

Complexity Analysis

The precomputed Möbius array and the tables of \(\varphi^{g^2}\) and \(\varphi^{-g^2}\) have size \(O(\sqrt L)\), so memory usage is \(O(\sqrt L)\).

For one fixed \(g\), the scaled limit is \(L/g^2\), and the two sliding-window branches together cost \(O(\sqrt{L/g^2})=O(\sqrt L/g)\). Summing over all \(g\le \sqrt L\) yields

\[ O\!\left(\sum_{g\le \sqrt L}\frac{\sqrt L}{g}\right) =O(\sqrt L\log L). \]

That is the reason the method can handle \(L=10^{14}\): it never iterates over all \(n\le L\), and it never enumerates all primitive pairs individually.

Footnotes and References

1. Project Euler problem page: Project Euler 989
2. Fibonacci numbers and Binet's formula: Wikipedia - Fibonacci number
3. Hensel lifting: Wikipedia - Hensel's lemma
4. Quadratic reciprocity and residue classes mod \(5\): Wikipedia - Quadratic reciprocity
5. Binary quadratic forms: Wikipedia - Binary quadratic form
6. Möbius inversion: Wikipedia - Möbius inversion formula