Problem 446: Retractions B

Project Euler Problem 446 Solution

EulerSolve provides an optimized solution for Project Euler Problem 446, Retractions B, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each modulus $m \gt 1$, consider linear maps $$f(x)\equiv ax+b \pmod{m},\qquad 0 \lt a \lt m,\quad 0\le b \lt m.$$ The map is a retraction when it is idempotent on every residue class: $$f(f(x))\equiv f(x)\pmod{m}\qquad\text{for all }x.$$ Let $R(m)$ be the number of such maps. Problem 446 asks for $$F(N)=\sum_{n=1}^{N}R(n^4+4),$$ with the checkpoint $F(1024)=77532377300600$. The challenge is to evaluate $F(10^7)\bmod(10^9+7)$ without factoring each number $n^4+4$ independently. Mathematical Approach Step 1: Convert the retraction condition into arithmetic conditions Composing the map once more gives $$f(f(x))\equiv a(ax+b)+b\equiv a^2x+ab+b \pmod{m}.$$ For this to equal $ax+b$ for every $x$, both the coefficient of $x$ and the constant correction must vanish modulo $m$. Therefore $$m\mid a(a-1),\qquad m\mid ab.$$ Write $$m=\prod_{i=1}^{r}p_i^{e_i}.$$ Because $\gcd(a,a-1)=1$, each prime power $p_i^{e_i}$ must divide exactly one of $a$ or $a-1$. Hence for every $p_i^{e_i}\parallel m$, $$a\equiv 0 \pmod{p_i^{e_i}}\qquad\text{or}\qquad a\equiv 1 \pmod{p_i^{e_i}}.$$ These independent choices correspond to the unitary divisors of $m$. If we define $$d=\gcd(a,m),$$ then $d\mid m$ and $\gcd(d,m/d)=1$, so $d$ is unitary. Conversely, every unitary divisor $d \lt m$ determines one admissible residue class for $a$....

Detailed mathematical approach

Problem Summary

For each modulus $m \gt 1$, consider linear maps

$$f(x)\equiv ax+b \pmod{m},\qquad 0 \lt a \lt m,\quad 0\le b \lt m.$$

The map is a retraction when it is idempotent on every residue class:

$$f(f(x))\equiv f(x)\pmod{m}\qquad\text{for all }x.$$

Let $R(m)$ be the number of such maps. Problem 446 asks for

$$F(N)=\sum_{n=1}^{N}R(n^4+4),$$

with the checkpoint $F(1024)=77532377300600$. The challenge is to evaluate $F(10^7)\bmod(10^9+7)$ without factoring each number $n^4+4$ independently.

Mathematical Approach

Step 1: Convert the retraction condition into arithmetic conditions

Composing the map once more gives

$$f(f(x))\equiv a(ax+b)+b\equiv a^2x+ab+b \pmod{m}.$$

For this to equal $ax+b$ for every $x$, both the coefficient of $x$ and the constant correction must vanish modulo $m$. Therefore

$$m\mid a(a-1),\qquad m\mid ab.$$

Write

$$m=\prod_{i=1}^{r}p_i^{e_i}.$$

Because $\gcd(a,a-1)=1$, each prime power $p_i^{e_i}$ must divide exactly one of $a$ or $a-1$. Hence for every $p_i^{e_i}\parallel m$,

$$a\equiv 0 \pmod{p_i^{e_i}}\qquad\text{or}\qquad a\equiv 1 \pmod{p_i^{e_i}}.$$

These independent choices correspond to the unitary divisors of $m$. If we define

$$d=\gcd(a,m),$$

then $d\mid m$ and $\gcd(d,m/d)=1$, so $d$ is unitary. Conversely, every unitary divisor $d \lt m$ determines one admissible residue class for $a$.

Now fix such an $a$, and write

$$a=d\,a_1,\qquad m=d\,m_1,\qquad \gcd(a_1,m_1)=1.$$

The condition $m\mid ab$ becomes

$$d\,m_1\mid d\,a_1b \iff m_1\mid a_1b \iff m_1\mid b,$$

because $a_1$ is invertible modulo $m_1$. Among the residues $0\le b\lt m=d\,m_1$, exactly $d$ are multiples of $m_1$. Therefore each admissible $a$ contributes exactly $d$ values of $b$.

Summing over all unitary divisors gives

$$R(m)=\sum_{\substack{d\parallel m\\ d \lt m}} d=\left(\sum_{d\parallel m}d\right)-m.$$

Thus

$$\boxed{R(m)=\sigma^*(m)-m},$$

where $\sigma^*(m)$ is the sum of unitary divisors. If $m=\prod p^e$, then

$$\sigma^*(m)=\prod_{p^e\parallel m}(1+p^e).$$

Step 2: Factor $n^4+4$ in a useful way

The key algebraic identity is the Sophie Germain factorization

$$n^4+4=(n^2-2n+2)(n^2+2n+2).$$

It is convenient to write

$$A_n=(n-1)^2+1,\qquad B_n=(n+1)^2+1,$$

so that

$$n^4+4=A_nB_n.$$

Define also

$$Q(m)=\sigma^*(m^2+1).$$

Then $A_n$ and $B_n$ are exactly the values $m^2+1$ at $m=n-1$ and $m=n+1$. The whole problem reduces to precomputing $Q(m)$ for $0\le m\le N+1$.

Step 3: Determine the gcd of the two factors

Any common divisor $g$ of $A_n$ and $B_n$ must divide their difference:

$$g\mid (B_n-A_n)=4n.$$

Also

$$A_n\equiv B_n\equiv 2 \pmod{n}.$$

If an odd prime $p$ divided both $A_n$ and $B_n$, then from $p\mid 4n$ we would get $p\mid n$, but then $A_n\equiv 2\pmod p$, which is impossible. So the gcd can only be a power of $2$.

If $n$ is odd, then $n-1$ and $n+1$ are even, hence $A_n$ and $B_n$ are odd, so

$$\gcd(A_n,B_n)=1.$$

If $n$ is even, then $n-1$ and $n+1$ are odd, so

$$A_n\equiv B_n\equiv 2 \pmod{8}.$$

Therefore each factor contains exactly one power of $2$, and

$$\gcd(A_n,B_n)=2.$$

Step 4: Derive the odd/even formulas for $\sigma^*(n^4+4)$

When $n$ is odd, the factors $A_n$ and $B_n$ are coprime, and $\sigma^*$ is multiplicative on coprime inputs. Hence

$$\sigma^*(n^4+4)=\sigma^*(A_n)\sigma^*(B_n)=Q(n-1)Q(n+1).$$

When $n$ is even, write

$$A_n=2u,\qquad B_n=2v,$$

with $u$ and $v$ odd and coprime. Then

$$Q(n-1)=\sigma^*(2u)=(1+2)\sigma^*(u)=3\sigma^*(u),$$

$$Q(n+1)=\sigma^*(2v)=3\sigma^*(v).$$

But now

$$n^4+4=A_nB_n=4uv,$$

and since $4$, $u$, and $v$ are pairwise coprime,

$$\sigma^*(n^4+4)=\sigma^*(4)\sigma^*(u)\sigma^*(v)=5\sigma^*(u)\sigma^*(v).$$

Comparing the two expressions yields the constant correction

$$\sigma^*(n^4+4)=\frac{5}{9}Q(n-1)Q(n+1)\qquad(n\text{ even}).$$

So the final closed form is

$$\boxed{\sigma^*(n^4+4)= \begin{cases} Q(n-1)Q(n+1), & n\text{ odd},\\[4pt] \dfrac{5}{9}Q(n-1)Q(n+1), & n\text{ even}. \end{cases}}$$

Step 5: Sieve all values of $Q(m)=\sigma^*(m^2+1)$ at once

For an odd prime $p$ to divide $m^2+1$, we must have

$$m^2\equiv -1 \pmod{p}.$$

That means $-1$ is a quadratic residue modulo $p$, which happens only for

$$p\equiv 1 \pmod{4}.$$

The prime $2$ is special: $2\mid m^2+1$ exactly when $m$ is odd, and then $m^2+1\equiv 2\pmod 8$, so the exponent of $2$ is exactly $1$.

This leads to a quadratic-congruence sieve:

For $p=2$, update all odd $m$.

For each odd prime $p\equiv 1\pmod 4$, find one square root $r$ of $-1$ modulo $p$. Then the only possible indices are the two residue classes

$$m\equiv r \pmod p,\qquad m\equiv -r \pmod p.$$

For every such $m$, divide $m^2+1$ by $p$ repeatedly to recover the exact exponent $p^e\parallel (m^2+1)$, and multiply the running value of $Q(m)$ by $1+p^e$.

After all relevant primes up to $N+1$ have been processed, any remaining factor larger than $1$ must itself be prime. Indeed, two unfactored primes both exceeding $N+1$ would have product larger than $(N+1)^2+1$, but $m^2+1\le (N+1)^2+1$. So one final factor $1+\text{leftover}$ completes $Q(m)$.

Step 6: Worked examples

For $n=2$, we have

$$n^4+4=20,\qquad A_2=2,\qquad B_2=10.$$

Thus

$$Q(1)=\sigma^*(2)=3,\qquad Q(3)=\sigma^*(10)=(1+2)(1+5)=18.$$

Since $n$ is even,

$$\sigma^*(20)=\frac{5}{9}\cdot 3\cdot 18=30,$$

and therefore

$$R(20)=30-20=10.$$

For $n=3$, the factors are coprime:

$$n^4+4=85,\qquad A_3=5,\qquad B_3=17.$$

Hence

$$\sigma^*(85)=Q(2)Q(4)=\sigma^*(5)\sigma^*(17)=6\cdot 18=108,$$

$$R(85)=108-85=23.$$

How the Code Works

The C++, Python, and Java implementations allocate arrays for the unfactored remainder of $m^2+1$ and for the running value of $Q(m)$ modulo $10^9+7$, for every $0\le m\le N+1$. Including $m=0$ makes the first summand $n=1$ work without a special case, since $0^2+1=1$ and $Q(0)=1$.

They next build a prime table up to $N+1$. The prime $2$ is handled by scanning odd indices. For each odd prime $p\equiv 1\pmod 4$, the implementation computes a square root of $-1$ modulo $p$ with the Tonelli-Shanks algorithm, then walks both arithmetic progressions $r,r+p,r+2p,\dots$ and $-r,-r+p,-r+2p,\dots$.

Whenever an index is hit, the implementation strips the full $p$-adic exponent from the current remainder and multiplies the running unitary-divisor product by $1+p^e$. After the sieve phase, any leftover prime factor is inserted once more. Finally, for each $n$, the code combines the two precomputed values at $n-1$ and $n+1$, applies the $5/9$ adjustment when $n$ is even, subtracts $n^4+4$ modulo $10^9+7$, and accumulates the result.

Complexity Analysis

The prime preprocessing and the arrays over $0,\dots,N+1$ use $O(N)$ memory. The residue-class sweeps visit about $2N/p$ indices for each prime $p\equiv 1\pmod 4$, so the total work is near

$$\sum_{p\le N+1,\ p\equiv 1\!\!\!\pmod 4}\frac{N}{p}=O(N\log\log N)$$

in practice. The remaining arithmetic per visited index is constant apart from stripping a prime power, so the overall method runs in near $O(N\log\log N)$ time and $O(N)$ memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=446
Unitary divisor: Wikipedia - Unitary divisor
Sophie Germain identity: Wikipedia - Sophie Germain's identity
Quadratic residue of $-1$: Wikipedia - Quadratic residue
Tonelli-Shanks algorithm: Wikipedia - Tonelli-Shanks algorithm

Problem 446 source code

C++

#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u32 kMod = 1'000'000'007U;

struct Options {
    u32 n = 10'000'000U;
    bool run_checkpoints = true;
};

bool parse_u32_after_prefix(const std::string& arg, const std::string& prefix, u32& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        const unsigned long long parsed = std::stoull(tail);
        if (parsed > static_cast<unsigned long long>(std::numeric_limits<u32>::max())) {
            return false;
        }
        out = static_cast<u32>(parsed);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_u32_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1U;
}

u32 mod_pow(u64 base, u64 exp, const u32 mod = kMod) {
    u64 result = 1ULL;
    base %= mod;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = (result * base) % mod;
        }
        base = (base * base) % mod;
        exp >>= 1ULL;
    }
    return static_cast<u32>(result);
}

u32 tonelli_shanks(u32 n, u32 p) {
    if (p == 2U) {
        return n & 1U;
    }
    if (n == 0U) {
        return 0U;
    }
    if (mod_pow(n, (p - 1U) / 2U, p) != 1U) {
        return 0U;
    }
    if ((p & 3U) == 3U) {
        return mod_pow(n, (p + 1U) / 4U, p);
    }

    u32 q = p - 1U;
    u32 s = 0U;
    while ((q & 1U) == 0U) {
        q >>= 1U;
        ++s;
    }

    u32 z = 2U;
    while (mod_pow(z, (p - 1U) / 2U, p) != p - 1U) {
        ++z;
    }

    u64 m = s;
    u64 c = mod_pow(z, q, p);
    u64 t = mod_pow(n, q, p);
    u64 r = mod_pow(n, (q + 1U) / 2U, p);

    while (t != 1U) {
        u64 tt = t;
        u64 i = 0U;
        while (tt != 1U && i < m) {
            tt = (tt * tt) % p;
            ++i;
        }

        const u64 shift = m - i - 1U;
        const u64 b = mod_pow(static_cast<u32>(c), 1ULL << shift, p);
        r = (r * b) % p;
        const u64 b2 = (b * b) % p;
        t = (t * b2) % p;
        c = b2;
        m = i;
    }

    return static_cast<u32>(r);
}

u64 brute_f_1024_exact() {
    auto r_exact = [](u64 n) -> u64 {
        u64 m = n;
        u64 q = 1ULL;
        for (u64 p = 2ULL; p * p <= m; ++p) {
            if (m % p != 0ULL) {
                continue;
            }
            u64 pe = 1ULL;
            while (m % p == 0ULL) {
                m /= p;
                pe *= p;
            }
            q *= (1ULL + pe);
        }
        if (m > 1ULL) {
            q *= (1ULL + m);
        }
        return q - n;
    };

    u64 sum = 0ULL;
    for (u64 n = 1ULL; n <= 1024ULL; ++n) {
        const u64 v = n * n * n * n + 4ULL;
        sum += r_exact(v);
    }
    return sum;
}

u32 solve_mod(const u32 n_limit) {
    const u32 max_m = n_limit + 1U;
    std::vector<u64> remaining(static_cast<std::size_t>(max_m) + 1U, 0ULL);
    std::vector<u32> q_value(static_cast<std::size_t>(max_m) + 1U, 1U);

    for (u32 m = 0U; m <= max_m; ++m) {
        remaining[m] = static_cast<u64>(m) * static_cast<u64>(m) + 1ULL;
    }

    std::vector<u32> spf(static_cast<std::size_t>(max_m) + 1U, 0U);
    std::vector<u32> primes;
    primes.reserve(static_cast<std::size_t>(max_m / 10U));

    for (u32 i = 2U; i <= max_m; ++i) {
        if (spf[i] == 0U) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (const u32 p : primes) {
            const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
            if (v > max_m || p > spf[i]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }

    const auto apply_prime_to_index = [&](const u32 p, const u32 idx) {
        u64 value = remaining[idx];
        if (value % p != 0ULL) {
            return;
        }
        u64 p_power_mod = 1ULL;
        do {
            value /= p;
            p_power_mod = (p_power_mod * p) % kMod;
        } while (value % p == 0ULL);
        remaining[idx] = value;
        const u32 term = static_cast<u32>((1ULL + p_power_mod) % kMod);
        q_value[idx] = static_cast<u32>((static_cast<u64>(q_value[idx]) * term) % kMod);
    };

    for (const u32 p : primes) {
        if (p == 2U) {
            for (u32 idx = 1U; idx <= max_m; idx += 2U) {
                apply_prime_to_index(2U, idx);
            }
            continue;
        }
        if ((p & 3U) != 1U) {
            continue;
        }

        const u32 root = tonelli_shanks(p - 1U, p);
        if (root == 0U) {
            continue;
        }
        const u32 root2 = (p - root) % p;

        for (u32 idx = root; idx <= max_m; idx += p) {
            apply_prime_to_index(p, idx);
        }
        if (root2 != root) {
            for (u32 idx = root2; idx <= max_m; idx += p) {
                apply_prime_to_index(p, idx);
            }
        }
    }

    for (u32 m = 0U; m <= max_m; ++m) {
        const u64 leftover = remaining[m];
        if (leftover > 1ULL) {
            const u32 term = static_cast<u32>((1ULL + (leftover % kMod)) % kMod);
            q_value[m] = static_cast<u32>((static_cast<u64>(q_value[m]) * term) % kMod);
        }
    }

    const u32 inv9 = mod_pow(9U, static_cast<u64>(kMod) - 2ULL);
    const u32 even_adjust = static_cast<u32>((5ULL * inv9) % kMod);

    u64 answer = 0ULL;
    for (u32 n = 1U; n <= n_limit; ++n) {
        const u32 q1 = q_value[n - 1U];
        const u32 q2 = q_value[n + 1U];
        u64 q_total = (static_cast<u64>(q1) * q2) % kMod;
        if ((n & 1U) == 0U) {
            q_total = (q_total * even_adjust) % kMod;
        }

        const u64 n_mod = n % kMod;
        const u64 n2 = (n_mod * n_mod) % kMod;
        const u64 n4_plus_4 = (n2 * n2 + 4ULL) % kMod;

        const u64 r_mod = (q_total + kMod - n4_plus_4) % kMod;
        answer += r_mod;
        answer %= kMod;
    }

    return static_cast<u32>(answer);
}

bool run_checkpoints() {
    if (brute_f_1024_exact() != 77'532'377'300'600ULL) {
        std::cerr << "Checkpoint failed: exact F(1024)\n";
        return false;
    }
    if (solve_mod(1024U) != 376'757'876U) {
        std::cerr << "Checkpoint failed: F(1024) mod 1e9+7\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve_mod(options.n) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000000007
    n_limit = 10000000

    def mod_pow(base, exp, mod=MOD):
        r = 1; base %= mod
        while exp > 0:
            if exp & 1: r = r * base % mod
            base = base * base % mod
            exp >>= 1
        return r

    def tonelli_shanks(n, p):
        if p == 2: return n & 1
        if n == 0: return 0
        if mod_pow(n, (p-1)//2, p) != 1: return 0
        if p % 4 == 3: return mod_pow(n, (p+1)//4, p)
        q, s = p-1, 0
        while q % 2 == 0: q //= 2; s += 1
        z = 2
        while mod_pow(z, (p-1)//2, p) != p-1: z += 1
        m, c, t, r = s, mod_pow(z,q,p), mod_pow(n,q,p), mod_pow(n,(q+1)//2,p)
        while t != 1:
            tt, i = t, 0
            while tt != 1 and i < m: tt = tt*tt%p; i += 1
            b = mod_pow(c, 1<<(m-i-1), p)
            r = r*b%p; t = t*b%p*b%p; c = b*b%p; m = i
        return r

    max_m = n_limit + 1
    remaining = [0]*(max_m+1)
    q_value = [1]*(max_m+1)
    for m in range(max_m+1): remaining[m] = m*m+1

    spf = list(range(max_m+1))
    for p in range(2, int(max_m**0.5)+1):
        if spf[p] == p:
            for q in range(p*p, max_m+1, p):
                if spf[q] == q: spf[q] = p
    primes = [i for i in range(2, max_m+1) if spf[i] == i]

    def apply_prime(p, idx):
        v = remaining[idx]
        if v % p != 0: return
        pp = 1
        while v % p == 0: v //= p; pp = pp*p%MOD
        remaining[idx] = v
        q_value[idx] = q_value[idx] * ((1+pp)%MOD) % MOD

    for p in primes:
        if p == 2:
            for idx in range(1, max_m+1, 2): apply_prime(2, idx)
            continue
        if p % 4 != 1: continue
        root = tonelli_shanks(p-1, p)
        if root == 0: continue
        root2 = (p-root)%p
        for idx in range(root, max_m+1, p): apply_prime(p, idx)
        if root2 != root:
            for idx in range(root2, max_m+1, p): apply_prime(p, idx)

    for m in range(max_m+1):
        if remaining[m] > 1:
            q_value[m] = q_value[m] * ((1+remaining[m]%MOD)%MOD) % MOD

    inv9 = mod_pow(9, MOD-2)
    even_adj = 5*inv9%MOD
    answer = 0
    for n in range(1, n_limit+1):
        q1 = q_value[n-1]; q2 = q_value[n+1]
        qt = q1*q2%MOD
        if n % 2 == 0: qt = qt*even_adj%MOD
        nm = n%MOD; n4p4 = (nm*nm%MOD*nm%MOD*nm+4)%MOD
        answer = (answer + qt + MOD - n4p4) % MOD

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler446 {
    static final int MOD = 1000000007;

    static long modPow(long base, long exp, int mod) {
        long result = 1;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                result = (result * base) % mod;
            base = (base * base) % mod;
            exp >>= 1;
        }
        return result;
    }

    static int tonelliShanks(int n, int p) {
        if (p == 2)
            return n & 1;
        if (n == 0)
            return 0;
        if (modPow(n, (p - 1) / 2, p) != 1)
            return 0;
        if ((p & 3) == 3)
            return (int) modPow(n, (p + 1) / 4, p);

        int q = p - 1;
        int s = 0;
        while ((q & 1) == 0) {
            q >>= 1;
            s++;
        }

        int z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1) {
            z++;
        }

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

        while (t != 1) {
            long tt = t;
            long i = 0;
            while (tt != 1 && i < m) {
                tt = (tt * tt) % p;
                i++;
            }

            long shift = m - i - 1;
            long b = modPow((int) c, 1L << shift, p);
            r = (r * b) % p;
            long b2 = (b * b) % p;
            t = (t * b2) % p;
            c = b2;
            m = i;
        }

        return (int) r;
    }

    static void applyPrimeToIndex(int p, int idx, long[] remaining, int[] qValue) {
        long value = remaining[idx];
        if (value % p != 0)
            return;
        long pPowerMod = 1;
        do {
            value /= p;
            pPowerMod = (pPowerMod * p) % MOD;
        } while (value % p == 0);
        remaining[idx] = value;
        int term = (int) ((1 + pPowerMod) % MOD);
        qValue[idx] = (int) (((long) qValue[idx] * term) % MOD);
    }

    public static String solve() {
        int nLimit = 10000000;
        int maxM = nLimit + 1;
        long[] remaining = new long[maxM + 1];
        int[] qValue = new int[maxM + 1];

        for (int m = 0; m <= maxM; m++) {
            remaining[m] = (long) m * m + 1;
            qValue[m] = 1;
        }

        int[] spf = new int[maxM + 1];
        List<Integer> primes = new ArrayList<>(maxM / 10);

        for (int i = 2; i <= maxM; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > maxM || p > spf[i]) {
                    break;
                }
                spf[(int) v] = p;
            }
        }

        for (int p : primes) {
            if (p == 2) {
                for (int idx = 1; idx <= maxM; idx += 2) {
                    applyPrimeToIndex(2, idx, remaining, qValue);
                }
                continue;
            }
            if ((p & 3) != 1)
                continue;

            int root = tonelliShanks(p - 1, p);
            if (root == 0)
                continue;

            int root2 = (p - root) % p;

            for (int idx = root; idx <= maxM; idx += p) {
                applyPrimeToIndex(p, idx, remaining, qValue);
            }

            if (root2 != root) {
                for (int idx = root2; idx <= maxM; idx += p) {
                    applyPrimeToIndex(p, idx, remaining, qValue);
                }
            }
        }

        for (int m = 0; m <= maxM; m++) {
            long leftover = remaining[m];
            if (leftover > 1) {
                int term = (int) ((1 + (leftover % MOD)) % MOD);
                qValue[m] = (int) (((long) qValue[m] * term) % MOD);
            }
        }

        long inv9 = modPow(9, MOD - 2, MOD);
        long evenAdjust = (5 * inv9) % MOD;

        long answer = 0;
        for (int n = 1; n <= nLimit; n++) {
            int q1 = qValue[n - 1];
            int q2 = qValue[n + 1];
            long qTotal = ((long) q1 * q2) % MOD;
            if ((n & 1) == 0) {
                qTotal = (qTotal * evenAdjust) % MOD;
            }

            long nMod = n % MOD;
            long n2 = (nMod * nMod) % MOD;
            long n4Plus4 = (n2 * n2 + 4) % MOD;

            long rMod = (qTotal + MOD - n4Plus4) % MOD;
            answer = (answer + rMod) % MOD;
        }

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}