Problem 437: Fibonacci Primitive Roots

Project Euler Problem 437 Solution

EulerSolve provides an optimized solution for Project Euler Problem 437, Fibonacci Primitive Roots, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want the sum of all primes $p \lt 10^8$ for which there exists an element $g \in \mathbb{F}_p^{\times}$ that is both a primitive root modulo $p$ and satisfies the Fibonacci-type recurrence $$g^{n+2}\equiv g^{n+1}+g^n \pmod{p}.$$ Such an element is called a Fibonacci primitive root. A naive strategy would search through many primitive roots for every prime, but the recurrence collapses to one quadratic congruence, so each prime has at most two relevant candidates. Mathematical Approach Step 1: Reduce the recurrence to a quadratic congruence Because $g$ is a primitive root modulo a prime, it is nonzero modulo $p$. Dividing the recurrence by $g^n$ gives $$g^2\equiv g+1\pmod{p}.$$ Therefore a Fibonacci primitive root is exactly a primitive root that solves $$x^2-x-1\equiv 0\pmod{p}.$$ This is the key simplification: instead of checking infinitely many recurrence steps, we only need to study one fixed polynomial....

Detailed mathematical approach

Problem Summary

We want the sum of all primes $p \lt 10^8$ for which there exists an element $g \in \mathbb{F}_p^{\times}$ that is both a primitive root modulo $p$ and satisfies the Fibonacci-type recurrence

$$g^{n+2}\equiv g^{n+1}+g^n \pmod{p}.$$

Such an element is called a Fibonacci primitive root. A naive strategy would search through many primitive roots for every prime, but the recurrence collapses to one quadratic congruence, so each prime has at most two relevant candidates.

Mathematical Approach

Step 1: Reduce the recurrence to a quadratic congruence

Because $g$ is a primitive root modulo a prime, it is nonzero modulo $p$. Dividing the recurrence by $g^n$ gives

$$g^2\equiv g+1\pmod{p}.$$

Therefore a Fibonacci primitive root is exactly a primitive root that solves

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

This is the key simplification: instead of checking infinitely many recurrence steps, we only need to study one fixed polynomial.

Step 2: Decide when solutions exist

The discriminant of $x^2-x-1$ is

$$\Delta = 1+4 = 5.$$

For $p\neq 5$, the quadratic has solutions if and only if $5$ is a quadratic residue modulo $p$, equivalently

$$\left(\frac{5}{p}\right)=1.$$

Since $5\equiv 1\pmod 4$, quadratic reciprocity yields

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

so for odd primes $p\neq 5$ the only possible residue classes are

$$p\equiv 1 \text{ or } 4 \pmod{5}.$$

This explains the early residue-class filter in the implementation. The small primes are exceptional: $p=2$ and $p=3$ give no solutions, while for $p=5$ the polynomial becomes $(x-3)^2$, so $3$ is the unique candidate and it is primitive modulo $5$.

Step 3: Recover the two candidates from $\sqrt{5}$

If $s^2\equiv 5\pmod p$, then the inverse of $2$ modulo an odd prime is

$$2^{-1}\equiv \frac{p+1}{2}\pmod p,$$

and the quadratic roots are

$$g_{\pm}\equiv \frac{1\pm s}{2}\pmod p.$$

For $p\neq 5$ there are exactly two such roots. Their sum is $1$ and their product is $-1$ modulo $p$, but that algebraic relation does not guarantee primitiveness, so both candidates must still be tested.

Step 4: Primitive-root criterion

The multiplicative group $\mathbb{F}_p^{\times}$ is cyclic of order $p-1$. An element is primitive if and only if its order is exactly $p-1$. Write

$$p-1=\prod_{i=1}^{k} q_i^{e_i}$$

for the prime factorization of $p-1$. Only the distinct prime divisors $q_i$ are needed. The standard criterion is

$$g^{(p-1)/q_i}\not\equiv 1\pmod p\qquad \text{for every } i=1,\dots,k.$$

If even one of these powers is $1$, then the order of $g$ is a proper divisor of $p-1$. If none is $1$, the order cannot drop, so $g$ is a primitive root.

Worked Example: $p=11$

The prime $11$ passes the residue filter because $11\equiv 1\pmod 5$. A square root of $5$ modulo $11$ is $4$, since $4^2=16\equiv 5\pmod{11}$. Therefore

$$g_{\pm}\equiv \frac{1\pm 4}{2}\equiv 8,\ 4 \pmod{11}.$$

Now $11-1=10=2\cdot 5$. For $g=8$, the primitive-root checks are

$$8^{10/2}=8^5\equiv 10\not\equiv 1\pmod{11},\qquad 8^{10/5}=8^2\equiv 9\not\equiv 1\pmod{11}.$$

So $8$ has order $10$, hence it is a primitive root and $11$ contributes to the final sum. This also shows that the residue test is only a necessary condition: it tells us when the quadratic has roots, not whether those roots generate the full group.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. They first build an odd-only smallest-prime-factor sieve. That single preprocessing step both identifies primes up to the limit and later factorizes $p-1$ efficiently. The scan skips $2$ and $3$, adds $5$ directly, and rejects all other primes with $p\bmod 5\notin\{1,4\}$.

For the remaining primes, the implementation computes $\sqrt{5}\pmod p$. When $p\equiv 3\pmod 4$, a short exponentiation formula gives the square root immediately; otherwise it uses Tonelli-Shanks. From that square root it constructs the two roots $(1\pm \sqrt{5})/2$, extracts the distinct prime divisors of $p-1$, and applies the primitive-root criterion. If either candidate passes, the prime is added to the running total.

Complexity Analysis

Let $N=10^8$. The odd-only sieve uses $O(N)$ memory asymptotically, with roughly half the storage of a full sieve, and it is built in $O(N\log\log N)$ time. After that, each surviving prime needs only a small number of modular exponentiations, one modular square-root computation, and factorization of $p-1$ through the precomputed SPF table.

In practice the sieve and the prime scan dominate the runtime. The per-prime work stays modest because only two candidates are ever examined, and the primitive-root test uses only the distinct prime divisors of $p-1$. This is dramatically faster than searching through all primitive roots or checking the Fibonacci recurrence term by term.

Footnotes and References

Problem page: https://projecteuler.net/problem=437
Primitive root modulo $n$: Wikipedia — Primitive root modulo n
Legendre symbol and quadratic reciprocity: Wikipedia — Legendre symbol
Tonelli-Shanks algorithm: Wikipedia — Tonelli-Shanks algorithm
Finite fields and cyclic multiplicative groups: Wikipedia — Finite field

Problem 437 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using u64 = std::uint64_t;

struct Options {
    int limit = 100000000;
    int threads = 0;
    bool run_checkpoints = true;
};

struct Result {
    u64 sum = 0;
    int count = 0;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } 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_int_after_prefix(arg, "--limit=", options.limit) ||
            parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 2 && options.threads >= 0;
}

int choose_thread_count(int requested, std::size_t work_items) {
    if (work_items <= 1) {
        return 1;
    }
    int threads = requested;
    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
    }
    if (threads <= 0) {
        threads = 4;
    }
    if (threads > static_cast<int>(work_items)) {
        threads = static_cast<int>(work_items);
    }
    if (threads < 1) {
        threads = 1;
    }
    return threads;
}

std::vector<std::uint16_t> build_odd_spf(const int limit) {
    std::vector<std::uint16_t> spf(static_cast<std::size_t>(limit / 2 + 1), 0);
    const int root = static_cast<int>(std::sqrt(static_cast<double>(limit)));
    for (int i = 3; i <= root; i += 2) {
        if (spf[static_cast<std::size_t>(i >> 1)] != 0) {
            continue;
        }
        const int step = i << 1;
        for (int j = i * i; j <= limit; j += step) {
            std::uint16_t& slot = spf[static_cast<std::size_t>(j >> 1)];
            if (slot == 0) {
                slot = static_cast<std::uint16_t>(i);
            }
        }
    }
    return spf;
}

std::vector<int> collect_primes(const int limit, const std::vector<std::uint16_t>& spf) {
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(limit / 12));
    if (limit >= 2) {
        primes.push_back(2);
    }
    for (int p = 3; p <= limit; p += 2) {
        if (spf[static_cast<std::size_t>(p >> 1)] == 0) {
            primes.push_back(p);
        }
    }
    return primes;
}

inline u64 mul_mod(u64 a, u64 b, int mod) {
    return (a * b) % static_cast<u64>(mod);
}

int mod_pow(int base, u64 exp, int mod) {
    u64 result = 1ULL;
    u64 cur = static_cast<u64>(base % mod + mod) % static_cast<u64>(mod);
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = mul_mod(result, cur, mod);
        }
        cur = mul_mod(cur, cur, mod);
        e >>= 1ULL;
    }
    return static_cast<int>(result);
}

int tonelli_sqrt_5(int p) {
    if (p == 2) {
        return 1;
    }
    if (mod_pow(5, static_cast<u64>((p - 1) / 2), p) != 1) {
        return -1;
    }
    if ((p & 3) == 3) {
        return mod_pow(5, static_cast<u64>((p + 1) / 4), p);
    }

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

    int z = 2;
    while (mod_pow(z, static_cast<u64>((p - 1) / 2), p) != p - 1) {
        ++z;
    }

    u64 c = static_cast<u64>(mod_pow(z, static_cast<u64>(q), p));
    u64 x = static_cast<u64>(mod_pow(5, static_cast<u64>((q + 1) / 2), p));
    u64 t = static_cast<u64>(mod_pow(5, static_cast<u64>(q), p));
    int m = s;

    while (t != 1ULL) {
        int i = 1;
        u64 tt = mul_mod(t, t, p);
        while (i < m && tt != 1ULL) {
            tt = mul_mod(tt, tt, p);
            ++i;
        }

        const int shift = m - i - 1;
        const u64 b = static_cast<u64>(mod_pow(static_cast<int>(c), 1ULL << shift, p));
        x = mul_mod(x, b, p);
        t = mul_mod(t, mul_mod(b, b, p), p);
        c = mul_mod(b, b, p);
        m = i;
    }

    return static_cast<int>(x);
}

int factor_unique(int n, const std::vector<std::uint16_t>& spf, std::array<int, 12>& out) {
    int cnt = 0;
    int x = n;

    if ((x & 1) == 0) {
        out[static_cast<std::size_t>(cnt++)] = 2;
        while ((x & 1) == 0) {
            x >>= 1;
        }
    }

    while (x > 1) {
        int f = 0;
        if ((x & 1) == 1) {
            f = static_cast<int>(spf[static_cast<std::size_t>(x >> 1)]);
            if (f == 0) {
                f = x;
            }
        } else {
            f = 2;
        }

        out[static_cast<std::size_t>(cnt++)] = f;
        while (x % f == 0) {
            x /= f;
        }
    }

    return cnt;
}

bool is_primitive_root(int g, int p, const std::array<int, 12>& factors, int factor_count) {
    const int phi = p - 1;
    for (int i = 0; i < factor_count; ++i) {
        const int q = factors[static_cast<std::size_t>(i)];
        if (mod_pow(g, static_cast<u64>(phi / q), p) == 1) {
            return false;
        }
    }
    return true;
}

Result solve(int limit, int requested_threads) {
    const std::vector<std::uint16_t> spf = build_odd_spf(limit);
    const std::vector<int> primes = collect_primes(limit, spf);

    const int thread_count = choose_thread_count(requested_threads, primes.size());
    std::vector<u64> partial_sum(static_cast<std::size_t>(thread_count), 0ULL);
    std::vector<int> partial_count(static_cast<std::size_t>(thread_count), 0);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(thread_count));

    for (int t = 0; t < thread_count; ++t) {
        const std::size_t begin =
            primes.size() * static_cast<std::size_t>(t) / static_cast<std::size_t>(thread_count);
        const std::size_t end =
            primes.size() * static_cast<std::size_t>(t + 1) / static_cast<std::size_t>(thread_count);

        workers.emplace_back([&, begin, end, t]() {
            u64 local_sum = 0ULL;
            int local_count = 0;
            std::array<int, 12> factors{};

            for (std::size_t idx = begin; idx < end; ++idx) {
                const int p = primes[idx];
                if (p == 2 || p == 3) {
                    continue;
                }
                if (p == 5) {
                    local_sum += 5ULL;
                    ++local_count;
                    continue;
                }

                const int mod5 = p % 5;
                if (mod5 != 1 && mod5 != 4) {
                    continue;
                }

                const int sqrt5 = tonelli_sqrt_5(p);
                if (sqrt5 < 0) {
                    continue;
                }

                const int inv2 = (p + 1) / 2;
                const int r1 = static_cast<int>((static_cast<u64>((1 + sqrt5) % p) * inv2) % p);
                const int r2 = static_cast<int>((static_cast<u64>((1 + p - sqrt5) % p) * inv2) % p);

                const int fc = factor_unique(p - 1, spf, factors);
                if (is_primitive_root(r1, p, factors, fc) || is_primitive_root(r2, p, factors, fc)) {
                    local_sum += static_cast<u64>(p);
                    ++local_count;
                }
            }

            partial_sum[static_cast<std::size_t>(t)] = local_sum;
            partial_count[static_cast<std::size_t>(t)] = local_count;
        });
    }

    for (auto& worker : workers) {
        worker.join();
    }

    Result res;
    for (int t = 0; t < thread_count; ++t) {
        res.sum += partial_sum[static_cast<std::size_t>(t)];
        res.count += partial_count[static_cast<std::size_t>(t)];
    }
    return res;
}

bool run_checkpoints(int requested_threads) {
    const Result small = solve(10000, requested_threads);
    if (small.count != 323) {
        std::cerr << "Checkpoint failed: count below 10000\n";
        return false;
    }
    if (small.sum != 1480491ULL) {
        std::cerr << "Checkpoint failed: sum below 10000\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(options.threads)) {
        return 2;
    }

    const Result ans = solve(options.limit, options.threads);
    std::cout << ans.sum << '\n';
    return 0;
}

Python

import math

def solve():
    limit = 100000000

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

    # SPF sieve for odd numbers
    half = limit // 2 + 1
    spf = [0] * half
    root = int(limit**0.5)
    for i in range(3, root + 1, 2):
        if spf[i >> 1] == 0:
            step = 2 * i
            for j in range(i*i, limit + 1, step):
                if spf[j >> 1] == 0: spf[j >> 1] = i

    primes = [2]
    for p in range(3, limit + 1, 2):
        if spf[p >> 1] == 0: primes.append(p)

    def tonelli_sqrt5(p):
        if p == 2: return 1
        if mod_pow(5, (p-1)//2, p) != 1: return -1
        if p % 4 == 3: return mod_pow(5, (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
        c = mod_pow(z, q, p); x = mod_pow(5, (q+1)//2, p)
        t = mod_pow(5, q, p); m = s
        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)
            x = x*b%p; t = t*b%p*b%p; c = b*b%p; m = i
        return x

    def factor_unique(n):
        factors = []
        if n % 2 == 0:
            factors.append(2)
            while n % 2 == 0: n //= 2
        while n > 1:
            if n & 1:
                f = spf[n >> 1]
                if f == 0: f = n
            else: f = 2
            factors.append(f)
            while n % f == 0: n //= f
        return factors

    def is_prim_root(g, p, factors):
        phi = p - 1
        for q in factors:
            if mod_pow(g, phi // q, p) == 1: return False
        return True

    total = 0
    for p in primes:
        if p <= 3: continue
        if p == 5: total += 5; continue
        if p % 5 != 1 and p % 5 != 4: continue
        s5 = tonelli_sqrt5(p)
        if s5 < 0: continue
        inv2 = (p + 1) // 2
        r1 = (1 + s5) % p * inv2 % p
        r2 = (1 + p - s5) % p * inv2 % p
        factors = factor_unique(p - 1)
        if is_prim_root(r1, p, factors) or is_prim_root(r2, p, factors):
            total += p

    return str(total)

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

Java

public class Euler437 {
    static int[] spfArr;

    public static String solve() {
        int limit = 100000000;
        int half = limit / 2 + 1;
        spfArr = new int[half];
        int root = (int) Math.sqrt(limit);
        for (int i = 3; i <= root; i += 2)
            if (spfArr[i >> 1] == 0)
                for (int j = i * i; j <= limit; j += 2 * i)
                    if (spfArr[j >> 1] == 0)
                        spfArr[j >> 1] = i;
        int[] primes;
        {
            java.util.List<Integer> pl = new java.util.ArrayList<>();
            pl.add(2);
            for (int p = 3; p <= limit; p += 2)
                if (spfArr[p >> 1] == 0)
                    pl.add(p);
            primes = pl.stream().mapToInt(Integer::intValue).toArray();
        }
        long total = 0;
        for (int p : primes) {
            if (p <= 3)
                continue;
            if (p == 5) {
                total += 5;
                continue;
            }
            if (p % 5 != 1 && p % 5 != 4)
                continue;
            long s5 = tonelliSqrt5(p);
            if (s5 < 0)
                continue;
            long inv2 = (p + 1L) / 2;
            long r1 = (1 + s5) % p * inv2 % p, r2 = (1 + p - s5) % p * inv2 % p;
            int[] factors = factorUnique(p - 1);
            if (isPrimRoot(r1, p, factors) || isPrimRoot(r2, p, factors))
                total += p;
        }
        return String.valueOf(total);
    }

    static long modPow(long b, long e, long m) {
        long r = 1;
        b %= m;
        while (e > 0) {
            if ((e & 1) != 0)
                r = r * b % m;
            b = b * b % m;
            e >>= 1;
        }
        return r;
    }

    static long tonelliSqrt5(int p) {
        if (p == 2)
            return 1;
        if (modPow(5, (p - 1) / 2, p) != 1)
            return -1;
        if (p % 4 == 3)
            return modPow(5, (p + 1L) / 4, p);
        long q = p - 1;
        int s = 0;
        while (q % 2 == 0) {
            q /= 2;
            s++;
        }
        int z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1)
            z++;
        long c = modPow(z, q, p), x = modPow(5, (q + 1) / 2, p), t = modPow(5, q, p);
        int m = s;
        while (t != 1) {
            long tt = t;
            int i = 0;
            while (tt != 1 && i < m) {
                tt = tt * tt % p;
                i++;
            }
            long b = modPow(c, 1L << (m - i - 1), p);
            x = x * b % p;
            t = t * b % p * b % p;
            c = b * b % p;
            m = i;
        }
        return x;
    }

    static int[] factorUnique(int n) {
        java.util.List<Integer> f = new java.util.ArrayList<>();
        if (n % 2 == 0) {
            f.add(2);
            while (n % 2 == 0)
                n /= 2;
        }
        while (n > 1) {
            int ff = (n & 1) != 0 ? (spfArr[n >> 1] == 0 ? n : spfArr[n >> 1]) : 2;
            f.add(ff);
            while (n % ff == 0)
                n /= ff;
        }
        return f.stream().mapToInt(Integer::intValue).toArray();
    }

    static boolean isPrimRoot(long g, int p, int[] factors) {
        long phi = p - 1;
        for (int q : factors)
            if (modPow(g, phi / q, p) == 1)
                return false;
        return true;
    }

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