Problem 492: Exploding Sequence

Project Euler Problem 492 Solution

EulerSolve provides an optimized solution for Project Euler Problem 492, Exploding Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The task is to evaluate $$B(x,y,n)=\sum_{\substack{p\ \text{prime}\\ x\le p\le x+y}} \left(a_n \bmod p\right),$$ for the specific target $B(10^9,10^7,10^{15})$, where the sequence is defined by $$a_1=1,\qquad a_{t+1}=6a_t^2+10a_t+3.$$ The recurrence is explosive: every step squares the previous value, so computing $a_n$ directly is hopeless when $n=10^{15}$. The solution therefore avoids building the integer $a_n$ itself and instead works modulo each prime in the interval. Mathematical Approach The central idea is to turn the quadratic recurrence into a Lucas-sequence problem whose enormous index can be reduced modulo a small group order. Step 1: Affine Change of Variables Introduce $$b_t=6a_t+5.$$ Then the recurrence becomes $$b_{t+1}=6a_{t+1}+5=6(6a_t^2+10a_t+3)+5=(6a_t+5)^2-2=b_t^2-2,$$ with initial value $$b_1=6\cdot 1+5=11.$$ This is the key simplification. Instead of a quadratic polynomial with three terms, we now have repeated application of the much cleaner map $u\mapsto u^2-2$....

Detailed mathematical approach

Problem Summary

The task is to evaluate

$$B(x,y,n)=\sum_{\substack{p\ \text{prime}\\ x\le p\le x+y}} \left(a_n \bmod p\right),$$

for the specific target $B(10^9,10^7,10^{15})$, where the sequence is defined by

$$a_1=1,\qquad a_{t+1}=6a_t^2+10a_t+3.$$

The recurrence is explosive: every step squares the previous value, so computing $a_n$ directly is hopeless when $n=10^{15}$. The solution therefore avoids building the integer $a_n$ itself and instead works modulo each prime in the interval.

Mathematical Approach

The central idea is to turn the quadratic recurrence into a Lucas-sequence problem whose enormous index can be reduced modulo a small group order.

Step 1: Affine Change of Variables

Introduce

$$b_t=6a_t+5.$$

Then the recurrence becomes

$$b_{t+1}=6a_{t+1}+5=6(6a_t^2+10a_t+3)+5=(6a_t+5)^2-2=b_t^2-2,$$

with initial value

$$b_1=6\cdot 1+5=11.$$

This is the key simplification. Instead of a quadratic polynomial with three terms, we now have repeated application of the much cleaner map $u\mapsto u^2-2$.

Step 2: Recognize the Lucas Structure

Let

$$\alpha,\beta=\frac{11\pm\sqrt{117}}{2},\qquad \alpha+\beta=11,\qquad \alpha\beta=1.$$

Define the Lucas $V$-sequence with parameters $P=11$ and $Q=1$ by

$$V_m=\alpha^m+\beta^m.$$

Then

$$V_0=2,\qquad V_1=11,\qquad V_{m+1}=11V_m-V_{m-1}.$$

Because $\alpha\beta=1$, we also have the doubling identity

$$V_{2m}=\alpha^{2m}+\beta^{2m}=(\alpha^m+\beta^m)^2-2(\alpha\beta)^m=V_m^2-2.$$

Since $b_1=V_1$ and both sequences obey the same doubling rule, induction gives

$$b_t=V_{2^{t-1}}.$$

Therefore the original problem is equivalent to evaluating

$$a_n \bmod p \quad \text{from} \quad V_{2^{n-1}} \bmod p.$$

Step 3: Reduce the Giant Index

For the primes relevant to the target sum, $p>13$, so the discriminant

$$\Delta=11^2-4=117$$

is nonzero modulo $p$, and $6$ is invertible modulo $p$.

If $117$ is a quadratic residue modulo $p$, then $\sqrt{117}\in\mathbb{F}_p$, so $\alpha,\beta\in\mathbb{F}_p^\times$. Since $\beta=\alpha^{-1}$, the value

$$V_k=\alpha^k+\alpha^{-k}$$

depends only on $k\bmod(p-1)$.

If $117$ is a quadratic nonresidue modulo $p$, then $\alpha$ and $\beta$ live in $\mathbb{F}_{p^2}$, and the Frobenius map swaps them:

$$\alpha^p=\beta=\alpha^{-1}.$$

Hence

$$\alpha^{p+1}=1,$$

so $V_k$ depends only on $k\bmod(p+1)$.

Thus the enormous index

$$k=2^{n-1}$$

is reduced by

$$k\equiv 2^{n-1}\pmod{p-1}\quad\text{if }\left(\frac{117}{p}\right)=1,$$

$$k\equiv 2^{n-1}\pmod{p+1}\quad\text{if }\left(\frac{117}{p}\right)=-1.$$

This turns an index of size roughly $2^{10^{15}}$ into one smaller than $p+1$.

Step 4: Fast Doubling for the Lucas Term

After the index has been reduced, the Lucas value is computed with standard doubling identities:

$$V_{2m}=V_m^2-2,$$

$$V_{2m+1}=V_mV_{m+1}-11,$$

$$V_{2m+2}=V_{m+1}^2-2.$$

Processing the bits of $k$ from top to bottom keeps a pair of consecutive Lucas values and reaches $V_k \bmod p$ in $O(\log k)$ arithmetic steps.

Once $b_n \equiv V_k \pmod p$ is known, we recover the original sequence via

$$a_n\equiv (b_n-5)\cdot 6^{-1}\pmod p.$$

Step 5: Sum Over the Prime Interval

The interval $[x,x+y]$ is handled with a segmented sieve. First, all base primes up to $\sqrt{x+y}$ are generated. Then their multiples are marked inside the target interval, leaving exactly the primes that contribute to the sum.

For each prime $p$ in the segment, the algorithm performs three modular tasks:

$$\text{determine whether }117\text{ is a residue mod }p,$$

$$\text{compute }2^{n-1}\text{ modulo }p-1\text{ or }p+1,$$

$$\text{evaluate }V_k\text{ and convert it back to }a_n\bmod p.$$

Adding those residues yields the required value of $B(x,y,n)$.

Worked Example

A small example shows the reduction clearly. Take $n=4$ and the interval $[5,7]$, so the contributing primes are $5$ and $7$.

For $p=5$, $117\equiv 2\pmod 5$, which is a nonresidue, so the index is reduced modulo $p+1=6$:

$$k\equiv 2^{4-1}=8\equiv 2\pmod 6.$$

Then

$$V_2=11^2-2=119\equiv 4\pmod 5,$$

$$a_4\equiv (4-5)\cdot 6^{-1}\equiv (-1)\cdot 1\equiv 4\pmod 5.$$

For $p=7$, $117\equiv 5\pmod 7$, also a nonresidue, so the index is reduced modulo $p+1=8$:

$$k\equiv 8\equiv 0\pmod 8.$$

Therefore

$$V_0=2,\qquad a_4\equiv (2-5)\cdot 6^{-1}\equiv (-3)\cdot 6\equiv 3\pmod 7.$$

Hence

$$B(5,2,4)=4+3=7.$$

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first build a simple sieve up to $\sqrt{x+y}$, then use those base primes to mark a segmented interval $[x,x+y]$. Every unmarked entry is a prime that contributes one term to the final sum.

For each such prime, the implementation computes a Legendre-style residue test for $117$, chooses $p-1$ or $p+1$ as the period modulus, and evaluates $2^{n-1}$ modulo that period with fast modular exponentiation. It then performs Lucas fast doubling on a pair of consecutive values until it reaches the reduced index $k$, converts the result back from $b_n$ to $a_n$, and adds the residue to the running total.

No large integer value of $a_n$ is ever constructed. Only modular values are propagated, which is why the method stays fast even when $n=10^{15}$.

Complexity Analysis

Let $H=x+y$, and let $\pi(x,y)$ denote the number of primes in $[x,x+y]$. Building the base sieve up to $\sqrt{H}$ costs $O(\sqrt{H}\log\log H)$ time and $O(\sqrt{H})$ memory. Marking the segment costs $O(y\log\log H)$ time and $O(y)$ memory.

For each prime in the segment, the modular residue test, the reduction of $2^{n-1}$, and the Lucas doubling stage together require $O(\log n+\log p)$ arithmetic operations. The overall running time is therefore

$$O\!\left(\sqrt{H}\log\log H+y\log\log H+\pi(x,y)(\log n+\log H)\right),$$

with memory usage

$$O(\sqrt{H}+y).$$

Footnotes and References

Problem 492 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

namespace {

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

struct Options {
    u64 x = 1'000'000'000ULL;
    u64 y = 10'000'000ULL;
    u64 n = 1'000'000'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(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_u64_after_prefix(arg, "--x=", options.x)) {
            continue;
        }
        if (parse_u64_after_prefix(arg, "--y=", options.y)) {
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

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

u64 mod_inv_prime(const u64 a, const u64 p) { return mod_pow(a, p - 2ULL, p); }

u64 add_mod(const u64 a, const u64 b, const u64 p) {
    const u64 c = a + b;
    return (c >= p || c < a) ? (c % p) : c;
}

u64 sub_mod(const u64 a, const u64 b, const u64 p) {
    return (a >= b) ? (a - b) : (a + p - b);
}

u64 mul_mod(const u64 a, const u64 b, const u64 p) {
    return static_cast<u64>((static_cast<u128>(a) * b) % p);
}

std::vector<int> simple_primes_up_to(const int limit) {
    std::vector<bool> is_prime(static_cast<std::size_t>(limit + 1), true);
    if (limit >= 0) {
        is_prime[0] = false;
    }
    if (limit >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; static_cast<int64_t>(p) * p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int x = p * p; x <= limit; x += p) {
            is_prime[static_cast<std::size_t>(x)] = false;
        }
    }
    std::vector<int> primes;
    for (int i = 2; i <= limit; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

std::vector<u64> primes_in_segment(const u64 lo, const u64 hi) {
    const int root = static_cast<int>(std::sqrt(static_cast<long double>(hi))) + 1;
    const std::vector<int> base_primes = simple_primes_up_to(root);

    std::vector<bool> is_prime(static_cast<std::size_t>(hi - lo + 1ULL), true);
    for (const int p : base_primes) {
        u64 start = (lo + static_cast<u64>(p) - 1ULL) / static_cast<u64>(p) * static_cast<u64>(p);
        if (start < static_cast<u64>(p) * static_cast<u64>(p)) {
            start = static_cast<u64>(p) * static_cast<u64>(p);
        }
        for (u64 x = start; x <= hi; x += static_cast<u64>(p)) {
            is_prime[static_cast<std::size_t>(x - lo)] = false;
        }
    }
    if (lo == 0ULL) {
        if (!is_prime.empty()) {
            is_prime[0] = false;
        }
        if (is_prime.size() > 1U) {
            is_prime[1] = false;
        }
    } else if (lo == 1ULL) {
        is_prime[0] = false;
    }

    std::vector<u64> out;
    out.reserve(static_cast<std::size_t>((hi - lo + 1ULL) / std::log(static_cast<long double>(hi))));
    for (u64 x = lo; x <= hi; ++x) {
        if (is_prime[static_cast<std::size_t>(x - lo)]) {
            out.push_back(x);
        }
    }
    return out;
}

u64 lucas_V_mod(const u64 k, const u64 p, const u64 P = 11ULL) {
    // Fast doubling on V_n(P,1):
    // V_{2m}   = V_m^2 - 2
    // V_{2m+1} = V_m V_{m+1} - P
    // V_{2m+2} = V_{m+1}^2 - 2
    if (k == 0ULL) {
        return 2ULL % p;
    }

    u64 vm = 2ULL % p;   // V_0
    u64 vmp1 = P % p;    // V_1

    int msb = 63 - __builtin_clzll(k);
    for (int bit = msb; bit >= 0; --bit) {
        const u64 v2m = sub_mod(mul_mod(vm, vm, p), 2ULL % p, p);
        const u64 v2m1 = sub_mod(mul_mod(vm, vmp1, p), P % p, p);
        const u64 v2m2 = sub_mod(mul_mod(vmp1, vmp1, p), 2ULL % p, p);

        if (((k >> bit) & 1ULL) == 0ULL) {
            vm = v2m;
            vmp1 = v2m1;
        } else {
            vm = v2m1;
            vmp1 = v2m2;
        }
    }
    return vm;
}

u64 a_n_mod_p(const u64 n, const u64 p) {
    if (n == 1ULL) {
        return 1ULL % p;
    }

    // b_n = 6 a_n + 5, b_{n+1}=b_n^2-2, b_1=11.
    // b_n = V_{2^{n-1}}(11,1) in F_p / F_{p^2}.
    // The index can be reduced modulo p-1 if 117 is a residue, otherwise p+1.
    const u64 leg = mod_pow(117ULL % p, (p - 1ULL) / 2ULL, p);
    const bool residue = (leg == 1ULL);
    const u64 order_mod = residue ? (p - 1ULL) : (p + 1ULL);

    const u64 idx = mod_pow(2ULL, n - 1ULL, order_mod);
    const u64 b = lucas_V_mod(idx, p, 11ULL);

    const u64 inv6 = mod_inv_prime(6ULL, p);
    const u64 a = mul_mod(sub_mod(b, 5ULL % p, p), inv6, p);
    return a;
}

u64 B(const u64 x, const u64 y, const u64 n) {
    const u64 lo = x;
    const u64 hi = x + y;
    const std::vector<u64> primes = primes_in_segment(lo, hi);

    u64 total = 0ULL;
    for (const u64 p : primes) {
        total += a_n_mod_p(n, p);
    }
    return total;
}

bool run_checkpoints() {
    if (a_n_mod_p(6ULL, 1'000'000'007ULL) != 203'064'689ULL) {
        std::cerr << "Checkpoint failed: a_6 mod 1e9+7\n";
        return false;
    }
    if (a_n_mod_p(100ULL, 1'000'000'007ULL) != 456'482'974ULL) {
        std::cerr << "Checkpoint failed: a_100 mod 1e9+7\n";
        return false;
    }
    if (B(1'000'000'000ULL, 1'000ULL, 1'000ULL) != 23'674'718'882ULL) {
        std::cerr << "Checkpoint failed: B(1e9,1e3,1e3)\n";
        return false;
    }
    if (B(1'000'000'000ULL, 1'000ULL, 1'000'000'000'000'000ULL) != 20'731'563'854ULL) {
        std::cerr << "Checkpoint failed: B(1e9,1e3,1e15)\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 1;
    }

    std::cout << B(options.x, options.y, options.n) << '\n';
    return 0;
}

Python

import math

def simple_primes_up_to(limit):
    is_prime = [True] * (limit + 1)
    if limit >= 0: is_prime[0] = False
    if limit >= 1: is_prime[1] = False
    p = 2
    while p * p <= limit:
        if is_prime[p]:
            for x in range(p * p, limit + 1, p):
                is_prime[x] = False
        p += 1
    return [i for i in range(2, limit + 1) if is_prime[i]]

def primes_in_segment(lo, hi):
    root = int(math.sqrt(hi)) + 1
    base_primes = simple_primes_up_to(root)
    is_prime = [True] * (hi - lo + 1)
    for p in base_primes:
        start = (lo + p - 1) // p * p
        if start < p * p:
            start = p * p
        for x in range(start, hi + 1, p):
            is_prime[x - lo] = False
            
    if lo == 0:
        if len(is_prime) > 0: is_prime[0] = False
        if len(is_prime) > 1: is_prime[1] = False
    elif lo == 1:
        is_prime[0] = False
        
    out = []
    for x in range(lo, hi + 1):
        if is_prime[x - lo]:
            out.append(x)
    return out

def lucas_V_mod(k, p, P=11):
    if k == 0: return 2 % p
    vm = 2 % p
    vmp1 = P % p
    msb = k.bit_length() - 1
    for bit in range(msb, -1, -1):
        v2m = (vm * vm - 2) % p
        v2m1 = (vm * vmp1 - P) % p
        v2m2 = (vmp1 * vmp1 - 2) % p
        if ((k >> bit) & 1) == 0:
            vm, vmp1 = v2m, v2m1
        else:
            vm, vmp1 = v2m1, v2m2
    return vm

def a_n_mod_p(n, p):
    if n == 1: return 1 % p
    leg = pow(117 % p, (p - 1) // 2, p)
    order_mod = (p - 1) if leg == 1 else (p + 1)
    idx = pow(2, n - 1, order_mod)
    b = lucas_V_mod(idx, p, 11)
    inv6 = pow(6, p - 2, p)
    a = ((b - 5) % p * inv6) % p
    return a

def B(x, y, n):
    lo = x
    hi = x + y
    primes = primes_in_segment(lo, hi)
    total = 0
    for p in primes:
        total += a_n_mod_p(n, p)
    return total

def solve():
    x = 1000000000
    y = 10000000
    n = 1000000000000000
    return str(B(x, y, n))

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

Java

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

public class Euler492 {

    private static long modPow(long base, long exp, long mod) {
        if (mod == 1L)
            return 0L;
        long out = 1L % mod;
        base %= mod;
        while (exp > 0L) {
            if ((exp & 1L) != 0L) {
                out = mulMod(out, base, mod);
            }
            base = mulMod(base, base, mod);
            exp >>= 1L;
        }
        return out;
    }

    private static long modInvPrime(long a, long p) {
        return modPow(a, p - 2L, p);
    }

    private static long subMod(long a, long b, long p) {
        return (a >= b) ? (a - b) : (a + p - b);
    }

    private static long mulMod(long a, long b, long mod) {
        return ((a % mod) * (b % mod)) % mod;
    }

    private static List<Integer> simplePrimesUpTo(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrime[i] = true;

        for (long p = 2; p * p <= limit; ++p) {
            if (!isPrime[(int) p])
                continue;
            for (long x = p * p; x <= limit; x += p) {
                isPrime[(int) x] = false;
            }
        }

        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; ++i) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    private static List<Long> primesInSegment(long lo, long hi) {
        int root = (int) Math.sqrt(hi) + 1;
        List<Integer> basePrimes = simplePrimesUpTo(root);

        int len = (int) (hi - lo + 1);
        boolean[] isPrime = new boolean[len];
        for (int i = 0; i < len; i++)
            isPrime[i] = true;

        for (int p : basePrimes) {
            long start = (lo + p - 1) / p * p;
            if (start < (long) p * p)
                start = (long) p * p;
            for (long x = start; x <= hi; x += p) {
                isPrime[(int) (x - lo)] = false;
            }
        }

        if (lo == 0) {
            if (len > 0)
                isPrime[0] = false;
            if (len > 1)
                isPrime[1] = false;
        } else if (lo == 1) {
            if (len > 0)
                isPrime[0] = false;
        }

        List<Long> out = new ArrayList<>();
        for (long x = lo; x <= hi; ++x) {
            if (isPrime[(int) (x - lo)])
                out.add(x);
        }
        return out;
    }

    private static long lucasVMod(long k, long p, long P) {
        if (k == 0L)
            return 2L % p;

        long vm = 2L % p;
        long vmp1 = P % p;

        int msb = 63 - Long.numberOfLeadingZeros(k);
        for (int bit = msb; bit >= 0; --bit) {
            long v2m = subMod(mulMod(vm, vm, p), 2L % p, p);
            long v2m1 = subMod(mulMod(vm, vmp1, p), P % p, p);
            long v2m2 = subMod(mulMod(vmp1, vmp1, p), 2L % p, p);

            if (((k >> bit) & 1L) == 0L) {
                vm = v2m;
                vmp1 = v2m1;
            } else {
                vm = v2m1;
                vmp1 = v2m2;
            }
        }
        return vm;
    }

    private static long aNModP(long n, long p) {
        if (n == 1L)
            return 1L % p;

        long leg = modPow(117L % p, (p - 1L) / 2L, p);
        boolean residue = (leg == 1L);
        long orderMod = residue ? (p - 1L) : (p + 1L);

        long idx = modPow(2L, n - 1L, orderMod);
        long b = lucasVMod(idx, p, 11L);

        long inv6 = modInvPrime(6L, p);
        long a = mulMod(subMod(b, 5L % p, p), inv6, p);
        return a;
    }

    private static long B(long x, long y, long n) {
        long lo = x;
        long hi = x + y;
        List<Long> primes = primesInSegment(lo, hi);

        long total = 0L;
        for (long p : primes) {
            total += aNModP(n, p);
        }
        return total;
    }

    public static void main(String[] args) {
        long x = 1000000000L;
        long y = 10000000L;
        long n = 1000000000000000L;
        System.out.println(B(x, y, n));
    }
}