Detailed mathematical approach

Problem Summary

For each odd prime \(p\), let \(R(p)\) be the least positive integer \(n\) satisfying

$n^2-3n-1\equiv 0\pmod{p^2}.$

The goal is to compute

$S(L)=\sum_{p\le L} R(p)$

for \(L=10^7\). Useful checkpoints are \(S(10)=5\), \(S(100)=1752\), and \(S(1000)=6728355\).

Mathematical Approach

Write

$f(n)=n^2-3n-1.$

The implementation does not search over all \(n\). Instead it first determines when roots can exist modulo \(p\), computes those roots explicitly, and then lifts them once to modulo \(p^2\).

Step 1: Rewrite the Congruence

Multiplying by \(4\) and completing the square gives

$4f(n)=4n^2-12n-4=(2n-3)^2-13.$

Because \(p\) is odd, \(2\) is invertible modulo \(p\). Therefore

$f(n)\equiv 0\pmod{p}\iff (2n-3)^2\equiv 13\pmod{p}.$

So the congruence has roots modulo \(p\) exactly when \(13\) is a quadratic residue modulo \(p\).

Step 2: Which Primes Can Contribute?

For odd primes \(p\neq 13\), quadratic reciprocity gives

$\left(\frac{13}{p}\right)=\left(\frac{p}{13}\right),$

because \(13\equiv 1\pmod{4}\). The nonzero quadratic residues modulo \(13\) are

$1,\ 3,\ 4,\ 9,\ 10,\ 12.$

Hence roots can exist only for primes in those six residue classes modulo \(13\). This explains the residue-class filter used by the implementation before any expensive modular square-root work is attempted.

The prime \(p=13\) must be handled separately. Modulo \(13\), the congruence becomes

$ (2n-3)^2\equiv 0\pmod{13}, $

so there is the double root \(n\equiv 8\pmod{13}\). But

$f(8)=8^2-3\cdot 8-1=39,$

which is divisible by \(13\) but not by \(13^2=169\). Therefore the root modulo \(13\) does not lift to a root modulo \(13^2\), and \(p=13\) contributes nothing.

Step 3: Explicit Roots Modulo \(p\)

Assume now that \(p\neq 13\) and that \(13\) is a quadratic residue modulo \(p\). If

$s^2\equiv 13\pmod{p},$

then from \(2n-3\equiv \pm s\pmod{p}\) we obtain the two roots

$r_1\equiv \frac{3+s}{2}\pmod{p},\qquad r_2\equiv \frac{3-s}{2}\pmod{p}.$

A modular square root \(s\) is found with the Tonelli-Shanks algorithm. Since \(2^{-1}\equiv (p+1)/2\pmod{p}\), the division by \(2\) is just another modular multiplication.

Step 4: Lift Each Root to Modulo \(p^2\)

Let \(r\) be one of the roots modulo \(p\). Any lift to modulo \(p^2\) has the form

$n=r+tp$

for some \(t\in\{0,1,\dots,p-1\}\). Expand \(f(r+tp)\):

$f(r+tp)=f(r)+tp\,f'(r)+t^2p^2,$

where

$f'(x)=2x-3.$

Reducing modulo \(p^2\) removes the last term, so the condition \(f(r+tp)\equiv 0\pmod{p^2}\) becomes

$f(r)+tp\,f'(r)\equiv 0\pmod{p^2}.$

Because \(r\) is already a root modulo \(p\), the value \(f(r)\) is divisible by \(p\). Dividing by \(p\) yields the linear congruence

$t\,f'(r)\equiv -\frac{f(r)}{p}\pmod{p},$

and therefore

$t\equiv -\frac{f(r)/p}{f'(r)}\pmod{p}.$

For \(p\neq 13\), we have \(f'(r)=2r-3\equiv \pm s\not\equiv 0\pmod{p}\), so the inverse exists. Thus each root modulo \(p\) lifts to exactly one root modulo \(p^2\).

Step 5: Determine \(R(p)\)

The two roots \(r_1\) and \(r_2\) modulo \(p\) produce two lifted roots \(n_1\) and \(n_2\) modulo \(p^2\). The definition of \(R(p)\) asks for the least positive solution, so

$R(p)=\min(n_1,n_2).$

Summing this quantity over all contributing primes gives the desired value \(S(L)\).

Worked Example: \(p=3\)

Here \(13\equiv 1\pmod{3}\), so we may take \(s\equiv 1\). The two roots modulo \(3\) are

$r_1\equiv \frac{3+1}{2}\equiv 2\pmod{3},\qquad r_2\equiv \frac{3-1}{2}\equiv 1\pmod{3}.$

For \(r=2\),

$f(2)=4-6-1=-3,\qquad f'(2)=1.$

So

$t\equiv -\frac{-3/3}{1}\equiv 1\pmod{3},$

and the lifted root is

$n=2+1\cdot 3=5.$

For \(r=1\),

$f(1)=1-3-1=-3,\qquad f'(1)=-1\equiv 2\pmod{3}.$

Since \(2^{-1}\equiv 2\pmod{3}\),

$t\equiv -\frac{-3/3}{2}\equiv 2\pmod{3},$

which gives

$n=1+2\cdot 3=7.$

Therefore \(R(3)=5\), and because \(3\) is the only contributing prime up to \(10\), we recover the checkpoint \(S(10)=5\).

How the Code Works

The C++, Python, and Java implementations all use the same structure. They begin with a prime sieve up to \(L\). After excluding \(2\) and \(13\), they keep only primes whose residue modulo \(13\) is one of \(1,3,4,9,10,12\). For each remaining prime, the implementation computes a square root of \(13\) modulo \(p\) via Tonelli-Shanks, constructs the two roots modulo \(p\), applies the Hensel correction above to obtain the two roots modulo \(p^2\), and adds the smaller lifted value to the running sum. Fast modular exponentiation is reused for Euler-criterion tests, modular inverses, and the Tonelli-Shanks substeps.

Complexity Analysis

Generating all primes up to \(L\) with a sieve of Eratosthenes costs \(O(L\log\log L)\) time and \(O(L)\) memory. The residue-class filter removes roughly half of the odd primes immediately. Each surviving prime then needs only a constant amount of modular arithmetic plus one Tonelli-Shanks square-root computation, whose cost is polylogarithmic in \(p\). In practice the sieve dominates the memory usage, and the entire method is easily fast enough for \(L=10^7\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=457
2. Quadratic reciprocity: Wikipedia - Quadratic reciprocity
3. Quadratic residue and Legendre symbol: Wikipedia - Quadratic residue
4. Hensel's lemma: Wikipedia - Hensel's lemma
5. Tonelli-Shanks algorithm: Wikipedia - Tonelli-Shanks algorithm