Project Euler

Project Euler Problem 457: A Polynomial Modulo the Square of a Prime

Let f(n) = n^2 - 3n - 1. Let p be a prime. Let R(p) be the smallest positive integer n such that f(n) mod p^2 = 0, or R(p) = 0 if no such n exists. Let SR(L) = sum of R(p) for all primes p <= L. Fi...

Source sync Apr 19, 2026

Problem #0457

Level Level 16

Solved By 955

Languages C++, Python

Answer 2647787126797397063

Length 462 words

modular_arithmeticnumber_theoryalgebra

Let \(f(n) = n^2 - 3n - 1\).

Let \(p\) be a prime.

Let \(R(p)\) be the smallest positive integer \(n\) such that \(f(n) \bmod p^2 = 0\) if such an integer \(n\) exists, otherwise \(R(p) = 0\).

Let \(SR(L)\) be \(\sum R(p)\) for all primes not exceeding \(L\).

Find \(SR(10^7)\).

Project Euler Problem 457: A Polynomial Modulo the Square of a Prime

Mathematical Analysis

Step 1: Roots modulo p

We need n^2 - 3n - 1 = 0 (mod p^2).

First solve mod p: n^2 - 3n - 1 = 0 (mod p).

Using the quadratic formula mod p: n = (3 +/- sqrt(9 + 4)) / 2 = (3 +/- sqrt(13)) / 2 (mod p)

This requires:

p != 2 (handle separately): We need 2 to be invertible mod p.
13 must be a quadratic residue mod p (i.e., Legendre symbol (13/p) = 1), OR p = 13 (handle separately).

For p = 2: f(n) = n^2 - 3n - 1. Check n mod 4: f(0) = -1 = 3 (mod 4), f(1) = -3 = 1 (mod 4), f(2) = -3 = 1 (mod 4), f(3) = -1 = 3 (mod 4). No solution => R(2) = 0.

For p = 13: f(n) = n^2 - 3n - 1 = 0 (mod 13). Discriminant = 13 = 0 (mod 13). n = 3/2 = 3 * 7 = 21 = 8 (mod 13). Double root at n=8 mod 13. Lifting to mod 169: need Hensel’s lemma with f’(8) = 28 - 3 = 13 = 0 (mod 13). Degenerate case: check f(8) mod 169 = 64 - 24 - 1 = 39 = 313 != 0 (mod 169). No lift => R(13) = 0.

For odd primes p != 13 where (13/p) = 1: Compute sqrt(13) mod p, then n = (3 +/- sqrt(13)) * inv(2) mod p. This gives two roots n1, n2 mod p. Lift each to mod p^2 via Hensel’s lemma.

Step 2: Hensel’s Lemma (lifting from mod p to mod p^2)

Given root r mod p of f(n), and f’(r) = 2r - 3 != 0 (mod p): Lifted root = r - f(r) * inv(f’(r)) (mod p^2)

Since f’(r) = 2r - 3, and if p != 2, the derivative is nonzero at a simple root.

Step 3: Algorithm

Sieve primes up to 10^7.
For each prime p, check if (13/p) = 1 using Euler criterion or quadratic reciprocity.
If yes, find sqrt(13) mod p (Tonelli-Shanks).
Compute two roots mod p, lift to mod p^2.
Take the smallest positive root mod p^2 as R(p).
Sum all R(p).

Complexity

Prime sieve: O(N log log N) for N = 10^7
Per prime: O(log p) for Tonelli-Shanks
Total: roughly O(N) with small constants

Answer

\boxed{2647787126797397063}

C++ project_euler/problem_457/solution.cpp

/*
 * Project Euler Problem 457: A Polynomial Modulo the Square of a Prime
 *
 * f(n) = n^2 - 3n - 1
 * R(p) = smallest positive n with f(n) = 0 (mod p^2), or 0 if none.
 * SR(L) = sum of R(p) for primes p <= L.
 * Find SR(10^7).
 *
 * Approach:
 *   - Sieve primes up to L.
 *   - For each prime p (skip p=2, p=13), check if 13 is a QR mod p.
 *   - If yes, find sqrt(13) mod p via Tonelli-Shanks.
 *   - Compute roots of n^2 - 3n - 1 = 0 mod p using quadratic formula.
 *   - Lift roots to mod p^2 via Hensel's lemma.
 *   - Take smallest positive root as R(p) and accumulate.
 *
 * Compile: g++ -O2 -o solution solution.cpp
 */

#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <cstdint>
#include <chrono>

using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;

const int LIMIT = 10000000;

// Sieve of Eratosthenes
vector<int> sieve_primes(int lim) {
    vector<bool> is_prime(lim + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; (ll)i * i <= lim; i++) {
        if (is_prime[i]) {
            for (int j = i * i; j <= lim; j += i)
                is_prime[j] = false;
        }
    }
    vector<int> primes;
    for (int i = 2; i <= lim; i++)
        if (is_prime[i]) primes.push_back(i);
    return primes;
}

// Modular exponentiation
ll powmod(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    if (base < 0) base += mod;
    while (exp > 0) {
        if (exp & 1) result = (lll)result * base % mod;
        base = (lll)base * base % mod;
        exp >>= 1;
    }
    return result;
}

// Modular inverse using Fermat's little theorem (p prime)
ll modinv(ll a, ll p) {
    return powmod(a, p - 2, p);
}

// Modular inverse mod p^2 using extended approach
ll modinv_p2(ll a, ll p2) {
    // Use powmod with Euler's totient: phi(p^2) = p^2 - p
    // a^(-1) = a^(p^2 - p - 1) mod p^2
    // But simpler: use __int128 based extended gcd or just powmod
    // Since p2 might not be prime, use extended gcd
    ll g, x, y;
    // Extended GCD
    ll a0 = a % p2;
    if (a0 < 0) a0 += p2;
    ll old_r = a0, new_r = p2;
    ll old_s = 1, new_s = 0;
    while (new_r != 0) {
        ll q = old_r / new_r;
        ll tmp;
        tmp = old_r - q * new_r; old_r = new_r; new_r = tmp;
        tmp = old_s - q * new_s; old_s = new_s; new_s = tmp;
    }
    ll inv = old_s % p2;
    if (inv < 0) inv += p2;
    return inv;
}

// Tonelli-Shanks: find sqrt(n) mod p. Returns -1 if not a QR.
ll tonelli_shanks(ll n, ll p) {
    n %= p;
    if (n < 0) n += p;
    if (n == 0) return 0;
    if (powmod(n, (p - 1) / 2, p) != 1) return -1;

    if (p % 4 == 3) {
        return powmod(n, (p + 1) / 4, p);
    }

    ll q = p - 1, s = 0;
    while (q % 2 == 0) { q /= 2; s++; }

    ll z = 2;
    while (powmod(z, (p - 1) / 2, p) != p - 1) z++;

    ll m = s;
    ll c = powmod(z, q, p);
    ll t = powmod(n, q, p);
    ll r = powmod(n, (q + 1) / 2, p);

    while (true) {
        if (t == 1) return r;
        ll i = 1;
        ll temp = (lll)t * t % p;
        while (temp != 1) {
            temp = (lll)temp * temp % p;
            i++;
        }
        ll b = c;
        for (ll j = 0; j < m - i - 1; j++)
            b = (lll)b * b % p;
        m = i;
        c = (lll)b * b % p;
        t = (lll)t * c % p;
        r = (lll)r * b % p;
    }
}

int main() {
    auto t0 = chrono::high_resolution_clock::now();

    vector<int> primes = sieve_primes(LIMIT);
    printf("Sieved %zu primes up to %d\n", primes.size(), LIMIT);

    ull total = 0;
    int count_contributing = 0;

    for (int idx = 0; idx < (int)primes.size(); idx++) {
        ll p = primes[idx];
        if (p == 2 || p == 13) continue;

        ll sq = tonelli_shanks(13, p);
        if (sq < 0) continue;

        ll p2 = p * p;
        ll inv2 = modinv(2, p);

        ll r1 = (3 + sq) % p * inv2 % p;
        ll r2 = ((3 - sq % p + p) % p) * inv2 % p;

        ll best = p2;

        for (int ri = 0; ri < 2; ri++) {
            ll r = (ri == 0) ? r1 : r2;

            // f(r) mod p2
            ll fr = ((lll)r * r % p2 - 3 * r % p2 + p2 - 1) % p2;
            if (fr < 0) fr += p2;
            fr = (((lll)r * r - 3 * r - 1) % p2 + p2) % p2;

            ll fpr = ((2 * r - 3) % p + p) % p;

            if (fpr == 0) {
                // Degenerate: try all lifts
                for (ll k = 0; k < p; k++) {
                    ll rr = r + k * p;
                    ll val = (((lll)rr * rr - 3 * rr - 1) % p2 + p2) % p2;
                    if (val == 0) {
                        if (rr == 0) rr = p2;
                        if (rr < best) best = rr;
                    }
                }
                continue;
            }

            // Hensel lift
            ll fpr_p2 = ((2 * r - 3) % p2 + p2) % p2;
            ll inv_fpr = modinv_p2(fpr_p2, p2);
            ll lifted = (r - (lll)fr * inv_fpr % p2 + p2) % p2;
            // More careful computation
            lifted = ((ll)((lll)r - (lll)fr % p2 * inv_fpr % p2) % p2 + p2) % p2;
            // Even more careful:
            lll val = (lll)fr * inv_fpr;
            ll sub = (ll)(val % p2);
            lifted = ((r - sub) % p2 + p2) % p2;
            if (lifted == 0) lifted = p2;
            if (lifted < best) best = lifted;
        }

        if (best < p2) {
            total += best;
            count_contributing++;
        }
    }

    auto t1 = chrono::high_resolution_clock::now();
    double elapsed = chrono::duration<double>(t1 - t0).count();

    printf("Contributing primes: %d\n", count_contributing);
    printf("SR(%d) = %llu\n", LIMIT, total);
    printf("Time: %.2f s\n", elapsed);

    return 0;
}

Python project_euler/problem_457/solution.py

#!/usr/bin/env python3
"""
Project Euler Problem 457: A Polynomial Modulo the Square of a Prime

f(n) = n^2 - 3n - 1
R(p) = smallest positive n with f(n) = 0 (mod p^2), or 0 if none exists.
SR(L) = sum of R(p) for primes p <= L.

Find SR(10^7).
"""

import math
from time import time

def sieve_primes(limit):
    """Sieve of Eratosthenes returning list of primes up to limit."""
    is_prime = bytearray(b'\x01') * (limit + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, int(limit**0.5) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = bytearray(len(is_prime[i*i::i]))
    return [i for i in range(2, limit + 1) if is_prime[i]]

def tonelli_shanks(n, p):
    """Find sqrt(n) mod p using Tonelli-Shanks algorithm. Returns None if no square root."""
    if n % p == 0:
        return 0
    if pow(n, (p - 1) // 2, p) != 1:
        return None  # n is not a QR mod p

    if p % 4 == 3:
        r = pow(n, (p + 1) // 4, p)
        return r

    # Factor out powers of 2 from p-1
    q = p - 1
    s = 0
    while q % 2 == 0:
        q //= 2
        s += 1

    # Find a non-residue z
    z = 2
    while pow(z, (p - 1) // 2, p) != p - 1:
        z += 1

    m = s
    c = pow(z, q, p)
    t = pow(n, q, p)
    r = pow(n, (q + 1) // 2, p)

    while True:
        if t == 1:
            return r
        # Find least i such that t^(2^i) = 1
        i = 1
        temp = (t * t) % p
        while temp != 1:
            temp = (temp * temp) % p
            i += 1
        b = pow(c, 1 << (m - i - 1), p)
        m = i
        c = (b * b) % p
        t = (t * c) % p
        r = (r * b) % p

def solve(L=10**7):
    """Compute SR(L)."""
    t0 = time()
    primes = sieve_primes(L)
    print(f"Sieved {len(primes)} primes up to {L} in {time()-t0:.2f}s")

    total = 0
    count_contributing = 0

    for p in primes:
        if p == 2:
            # Check all n mod 4: f(n) = n^2 - 3n - 1 mod 4
            # f(0)=3, f(1)=1, f(2)=1, f(3)=3 mod 4 => no solution
            continue
        if p == 13:
            # Double root, f'(8) = 0 mod 13, f(8) = 39 != 0 mod 169 => no lift
            continue

        # Check if 13 is a QR mod p
        sq = tonelli_shanks(13, p)
        if sq is None:
            continue

        p2 = p * p
        inv2 = pow(2, p - 2, p)  # 2^(-1) mod p (Fermat's little theorem)

        # Two roots mod p
        r1 = ((3 + sq) * inv2) % p
        r2 = ((3 - sq) * inv2) % p

        best = p2  # Will find minimum positive root mod p^2

        for r in (r1, r2):
            # f(r) mod p should be 0, compute f(r) exactly for lifting
            fr = (r * r - 3 * r - 1) % p2
            fpr = (2 * r - 3) % p  # f'(r) mod p

            if fpr == 0:
                # Degenerate: check if f(r) = 0 mod p^2 already
                # If not, try all lifts r + k*p for k=0..p-1
                for k in range(p):
                    rr = r + k * p
                    if (rr * rr - 3 * rr - 1) % p2 == 0:
                        if rr == 0:
                            rr = p2
                        best = min(best, rr)
                continue

            # Hensel lift: r_lifted = r - f(r) * inv(f'(r)) mod p^2
            inv_fpr = pow(fpr, p - 2, p)
            # Compute lift carefully mod p^2
            # r_lifted = r - f(r) * inv(f'(r)) mod p^2
            # But we need inv(f'(r)) mod p^2 only to first order
            # Actually: lift = r - f(r)/f'(r) mod p^2
            # f(r)/f'(r) needs: f(r) is divisible by p, so f(r)/p is an integer
            # lift = (r - (fr * inv_fpr_p2)) % p2 where inv_fpr_p2 = inv of (2r-3) mod p^2
            # Simpler: use modular inverse mod p^2
            fpr_mod_p2 = (2 * r - 3) % p2
            # inv of fpr_mod_p2 mod p^2, but fpr might not be invertible mod p^2
            # Since fpr != 0 mod p, gcd(fpr_mod_p2, p^2) = 1, so invertible
            inv_fpr_p2 = pow(fpr_mod_p2, -1, p2)
            lifted = (r - fr * inv_fpr_p2) % p2
            if lifted == 0:
                lifted = p2
            best = min(best, lifted)

        if best < p2:
            total += best
            count_contributing += 1

    elapsed = time() - t0
    print(f"Contributing primes: {count_contributing}")
    print(f"SR({L}) = {total}")
    print(f"Time: {elapsed:.2f}s")
    return total

def generate_visualization(L=10**5):
    """Generate visualization of R(p) distribution for primes up to L."""
    primes = sieve_primes(L)

    rp_values = []
    rp_primes = []

    for p in primes:
        if p == 2 or p == 13:
            continue

        sq = tonelli_shanks(13, p)
        if sq is None:
            continue

        p2 = p * p
        inv2 = pow(2, p - 2, p)
        r1 = ((3 + sq) * inv2) % p
        r2 = ((3 - sq) * inv2) % p

        best = p2
        for r in (r1, r2):
            fr = (r * r - 3 * r - 1) % p2
            fpr_mod_p2 = (2 * r - 3) % p2
            fpr = (2 * r - 3) % p
            if fpr == 0:
                for k in range(p):
                    rr = r + k * p
                    if (rr * rr - 3 * rr - 1) % p2 == 0:
                        if rr == 0:
                            rr = p2
                        best = min(best, rr)
                continue
            inv_fpr_p2 = pow(fpr_mod_p2, -1, p2)
            lifted = (r - fr * inv_fpr_p2) % p2
            if lifted == 0:
                lifted = p2
            best = min(best, lifted)

        if best < p2:
            rp_values.append(best)
            rp_primes.append(p)

Project Euler Problem 457: A Polynomial Modulo the Square of a Prime

Problem Statement

Project Euler Problem 457: A Polynomial Modulo the Square of a Prime

Mathematical Analysis

Step 1: Roots modulo p

Step 2: Hensel’s Lemma (lifting from mod p to mod p^2)

Step 3: Algorithm

Complexity

Answer

Code