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