Problem 304: Primonacci

Project Euler Problem 304 Solution

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

Problem Summary Let $p_1,p_2,\dots,p_m$ be the first $m$ primes strictly larger than a huge starting value $A$. The task is to compute $$\sum_{i=1}^{m}F_{p_i}\pmod{M},$$ where $F_n$ is the Fibonacci sequence and the official parameters are extremely large. Mathematical Approach 1) Why prime generation must be segmented The primes we need live near $A\approx 10^{14}$. A naive primality test for each number in that region would be far too slow, and a full sieve up to $A$ is impossible in memory. The right tool is a segmented sieve. We scan blocks $$[L,H],$$ and mark composites inside that block using all base primes up to $\sqrt{H}$. 2) Why primes up to $\sqrt{H}$ are enough If $x\in [L,H]$ is composite, then $x=ab$ with $1<a\le b$. Hence $$a\le \sqrt{x}\le \sqrt{H}.$$ So every composite in the segment has a prime divisor not exceeding $\sqrt{H}$. Therefore, after marking multiples of all primes $p\le \sqrt{H}$, every remaining unmarked number in the segment is prime. 3) Where marking starts inside a segment For a fixed base prime $p$, the first multiple of $p$ inside the current segment is $$\left\lceil\frac{L}{p}\right\rceil p.$$ But we must also avoid marking the prime $p$ itself when $p\in[L,H]$....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, the next prime after $n$ is the smallest prime $p$ with $p \gt n$. The sequence $a$ walks from prime to prime above $10^{14}$: $a(1)$ is the next prime after $10^{14}$, and $a(n)$ the next prime after $a(n-1)$ for $n \gt 1$. The Fibonacci sequence is $f(0)=0$, $f(1)=1$, $f(n)=f(n-1)+f(n-2)$ for $n \gt 1$, and $b(n)=f(a(n))$. We must find

$$\sum_{n=1}^{100000} b(n) \pmod{1234567891011}.$$

So the indices are the first $100\,000$ primes above $10^{14}$, and we add the Fibonacci numbers there modulo $M=1234567891011$. The statement gives no sample value; every checkpoint below is one we built. ^[1]

The whole solution in one sentence: list those primes by sieving one window of four million integers instead of everything below $10^{14}$, turn each prime $p$ into $f(p) \bmod M$ with an identity that builds $f(2k)$ and $f(2k+1)$ from $f(k)$ and $f(k+1)$, and add the residues.

Mathematical Approach

A tiny version first: the first three primes after $5$ are $7$, $11$, $13$, their Fibonacci numbers $13$, $89$, $233$, and modulo $10$ the sum is $3+9+3\equiv5$. The real problem has the same two jobs, but neither may write a Fibonacci number down in full.

1. Two quantities that cannot be written down

Binet's formula gives ^[3]

$$f(n)=\frac{\varphi^{n}-\psi^{n}}{\sqrt5},$$

where $\varphi=(1+\sqrt5)/2$ and $\psi=1-\varphi$. As $|\psi|\lt1$, $f(n)\approx\varphi^{n}/\sqrt5$ has about $n\log_{10}\varphi$ decimal digits — for $n=10^{14}$ about $2.09\cdot10^{13}$ of them, some $21$ terabytes of text for one summand. A sieve over every integer up to $10^{14}$, one bit each, would need $12.5$ terabytes. The two obstacles get separate treatment: Fibonacci numbers only modulo $M$, and a sieve over a short interval above $10^{14}$ using only the small primes that mark its composites.

2. Why the residues are enough

Congruence modulo $M$ survives addition, subtraction and multiplication: if $x\equiv x'$ and $y\equiv y'$, then $x\pm y\equiv x'\pm y'$ and $xy\equiv x'y'$. ^[8] The Fibonacci identities and the running sum use only those three, so their values may be reduced at every step; indices and sieve positions stay ordinary integers.

That decides one implementation question. Here $M\lt2^{41}$ and

$$(M-1)^2=1524157877512879238820100,$$

about $8.3\cdot10^{4}$ times larger than $2^{64}$: a product of two residues does not fit in a 64-bit word, and the C++ and Python code adds two such squares before reducing.

3. Where the hundred thousand primes are

Primes near $x$ are spaced by about $\ln x$ on average, ^[6] and $\ln(10^{14})=32.236\dots$, so the first $100\,000$ primes above $10^{14}$ should fit in some $3.22$ million integers; computing them confirms it. Writing each as $10^{14}+r$ for readability, the sequence begins

$$a(1)=10^{14}+31,\quad a(2)=10^{14}+67,$$

$$a(3)=10^{14}+97,\quad a(4)=10^{14}+99,$$

and ends at $a(100000)=10^{14}+3235443$. The mean gap is $32.354$; $4133$ of the $99\,999$ gaps equal $2$, and the largest is $380$, between $10^{14}+1217707$ and $10^{14}+1218087$. A window of $4\,000\,000$ integers above $10^{14}$ holds $124\,000$ primes, so one suffices.

4. Sieving a window far from zero

Let $L=10^{14}+1$ and $H=L+4\,000\,000-1$. We want the primes in $[L,H]$ without sieving everything below $L$. The tool is the segmented sieve of Eratosthenes, ^[2] resting on one observation.

Claim. If $x\in[L,H]$ is composite, then $x$ has a prime factor $q\le\sqrt H$.

Proof. Let $q$ be the smallest prime factor of $x$ and write $x=qy$. Every prime factor of $y$ is at least $q$, and $y\gt1$ because $x$ is composite, so $q\le y$. Hence $q^{2}\le qy=x\le H$. $\square$

So crossing out, for every prime $p\le\sqrt H$, the multiples of $p$ in $[L,H]$ leaves exactly the primes of the window. Here $\sqrt H=10000000.199\dots$, so the base primes are the $664\,579$ primes up to $10^{7}$.

5. Where to start crossing out, and what it costs

The smallest multiple of a base prime $p$ at or above $L$ is $\lceil L/p\rceil\,p$, yet the sieve may begin later, at

$$\max\!\left(p^{2},\ \left\lceil \frac{L}{p}\right\rceil p\right),$$

which is what all three programs compute. The proof above is the reason: a composite $x$ need only be struck once, and the pass of its smallest prime factor $q$ does strike it, starting at $q^{2}\le x$; composite multiples of $p$ below $p^{2}$ fall to an earlier pass. Starting at $p^{2}$ also leaves $p$ unstruck when $p$ lies in the window.

Each base prime writes about $\Delta/p$ marks into a window of width $\Delta$, about $\Delta\sum_{p\le\sqrt H}1/p$ in all. Mertens' second theorem puts that sum at $\ln\ln\sqrt H$ plus a constant, ^[7] and it comes to $3.0414\dots$ for $p\le10^{7}$; the window duly receives $12\,165\,859$ marks, $3.04$ per position.

6. The Fibonacci addition formula

Indices of size $10^{14}$ put the recurrence itself out of reach. The way in is an identity that jumps:

$$f(m+n)=f(m)f(n+1)+f(m-1)f(n)$$

for $m\ge1$, $n\ge0$. ^[3]

Proof. Fix $m$ and induct on $n$. At $n=0$ the right side is $f(m)f(1)+f(m-1)f(0)=f(m)$; at $n=1$ it is $f(m)+f(m-1)=f(m+1)$. For $n\ge2$ both sides obey the same recurrence in $n$: the left because

$$f(m+n)=f(m+n-1)+f(m+n-2),$$

the right because it is a fixed linear combination of $f(n+1)$ and $f(n)$, which both obey it; and two sequences agreeing at $n=0,1$ under the same recurrence agree everywhere. $\square$

7. From the addition formula to doubling

Put $m=n=k$. Then $f(2k)=f(k)f(k+1)+f(k-1)f(k)$, and $f(k-1)=f(k+1)-f(k)$ turns this into

$$f(2k)=f(k)\bigl(2f(k+1)-f(k)\bigr).$$

Putting $m=k+1$, $n=k$ gives

$$f(2k+1)=f(k+1)^{2}+f(k)^{2}.$$

These are the standard fast-doubling identities, ^[4] and both also hold at $k=0$.

8. Fast doubling

Carry the pair $\bigl(f(k),f(k+1)\bigr)$ instead of one value. Section 7 turns it into $\bigl(f(2k),f(2k+1)\bigr)$, and one addition into the pair at $2k+1$,

$$\bigl(f(2k+1),\,f(2k)+f(2k+1)\bigr).$$

So from the pair at $k$ we reach the pair at $2k$ or at $2k+1$ — the two moves that append a binary digit $0$ or $1$ to $k$. Reading the digits of $p$ from the top down, starting from $(0,1)$ at index $0$, reaches the pair at $p$ in as many steps as $p$ has binary digits.

Every prime here satisfies $10^{14}\lt p\lt2^{47}$, so each index needs $47$ steps, each using three modular multiplications — $f(k)\cdot\bigl(2f(k+1)-f(k)\bigr)$, $f(k)^{2}$, $f(k+1)^{2}$ — so $141$ per prime and $14\,100\,000$ in all.

9. The whole task worked out small

Replace $10^{14}$ by $100$, the count by $5$, the modulus by $1000$. For the first stage, sieve $[101,120]$: the base primes are $2$, $3$, $5$, $7$, the first multiples crossed out $102$, $102$, $105$, $105$, and the survivors exactly $101$, $103$, $107$, $109$, $113$. Their Fibonacci residues are

$$101,\quad 277,\quad 913,\quad 9,\quad 333 \pmod{1000},$$

with sum $1633\equiv633\pmod{1000}$, matching the last digits of the exact integers $f(101)=573147844013817084101$ and $f(103)=1500520536206896083277$.

The second stage is as easy by hand. For $k=10$, binary digits $1,0,1,0$, the pairs are

$$(0,1)\to(1,1)\to(1,2)\to(5,8)\to(55,89),$$

the pairs at indices $0,1,2,5,10$, so $f(10)=55$ — the first checkpoint built into the C++ solution.

10. The Pisano period saves less than one might hope

Modulo $3$ the residues run $0,1,1,2,0,2,2,1,0,1,\dots$: the pair $(0,1)$ recurs after $8$ steps, and a pair determines all later ones, so everything repeats. The step $(a,b)\mapsto(b,a+b)$ is reversible, $(u,v)$ coming from $(v-u,u)$, so the finitely many pairs must cycle back to the start; the first return of $(0,1)$ is the Pisano period $\pi(M)$, and the periods of coprime factors combine by least common multiple. ^[5] Factoring

$$M=3\cdot7\cdot13\cdot67\cdot107\cdot630803,$$

and running the recurrence for each factor gives the periods $8$, $16$, $28$, $136$, $72$, $1\,261\,608$, with least common multiple

$$\pi(M)=900\,788\,112.$$

Replacing $p$ by $p\bmod\pi(M)$ shortens the index to $30$ binary digits, so $90$ multiplications instead of $141$ — a third saved for the price of factoring $M$. The published solutions skip it; the next section is worth more.

11. Stepping along the prime gaps

Fast doubling treats the indices as unrelated, but consecutive ones differ by a gap $g=a(n+1)-a(n)$ of at most $380$ (Section 3), and the addition formula moves a pair forward by $g$ in one shot:

$$\begin{aligned} f(N+g)&=f(g)f(N+1)+f(g-1)f(N),\\ f(N+g+1)&=f(g+1)f(N+1)+f(g)f(N). \end{aligned}$$

With $f(0),\dots,f(382)\bmod M$ tabulated by the plain recurrence, one fast-doubling evaluation for $a(1)$ and four multiplications per gap suffice: $141+4\cdot99\,999=400\,137$ in all, some $35$ times fewer than Section 8 needs. We implemented it as an independent check; it gives the same residue as all three published programs.

How the Code Works

All three solutions share the same three stages; Python and Java are ports of the C++ file. ^[10] The names below are its own, which Java spells sieveBasePrimes, firstPrimesAfter and fibPairMod.

sieve_base_primes(limit) is a textbook sieve on a byte array: for each $p$ with $p^{2}\le$ limit still marked prime it clears p*p, p*p+p, …, and collects the survivors. C++ uses limit = 20000000 and returns $1\,270\,607$ primes, Python and Java $11\,000\,000$ and $726\,517$. For the official input both are generous: the consumer stops at $p^{2}\gt H$ and uses only the $664\,579$ primes up to $10^{7}$. The list is fixed, though, so the sieve holds only while it reaches $\sqrt H$; the caveats below show what happens otherwise.

first_primes_after(start, count) implements Sections 4 and 5. It sets low = start + 1 and segment_size = 4000000, then loops until count primes are collected. Each round sets high = low + segment_size - 1, fills the window flags with ones, and for every base prime p breaks out once p*p > high, otherwise sets first to the larger of p*p and the first multiple of p at or above low, clearing flags in steps of p. A final upward scan appends low + i for each survivor, so the primes arrive sorted, as $a(1),a(2),\dots$ require. With the official input the body runs once, the first window holding $124\,000$ primes.

fib_pair_mod(n, mod) is Section 8 as a recursion rather than a digit loop. It returns $\bigl(f(n),f(n+1)\bigr)$ modulo mod, with $(0,1\bmod m)$ at n == 0; otherwise it recurses on n >> 1, forms c $=a(2b-a)$ and d $=a^{2}+b^{2}$ from the returned pair $(a,b)$, and returns (c, d) for even n, (d, c + d) for odd n; unwinding it traverses the binary digits from the top down. There are $47$ non-base levels, $48$ live calls with the base case, so little stack is used. The doubling term is (two_b + mod - a) % mod, not 2*b - a, because the C++ residues are unsigned: adding mod keeps the difference non-negative.

solve calls first_primes_after and adds the first component of fib_pair_mod(p, mod) for each prime, reducing with if (sum >= mod) sum %= mod instead of dividing every term — safe because the running sum and the addend are each below $M$, so their sum is below $2M\lt2^{63}$ and fits an unsigned C++ word as well as Java's signed long.

Integer widths, and the self-test

A product of two residues reaches $1.5\cdot10^{24}$ (Section 2). For the three products of fib_pair_mod the C++ file casts one operand to unsigned __int128, a GCC and Clang extension, and reduces back to 64 bits; it forms $a^{2}+b^{2}$ in that width and reduces the sum once, whereas Java reduces each square separately. ^[9] Doubling a residue and adding two reduced residues stay in 64 bits. Java has no 128-bit type and its long stops at $2^{63}-1$, wrapping silently, so the port supplies mulMod, which builds the product by doubling a and adding it whenever the current bit of b is set; with $M\lt2^{41}$ every intermediate stays below $2^{42}$. Python needs no such care, its integers being unbounded, and halves the index with an ordinary shift where Java needs the unsigned n >>> 1.

Only the C++ file carries self-tests. run_checkpoints checks $f(10)=55$ modulo $10^{9}+7$, then compares solve with solve_bruteforce on start = 1000, count = 30, mod = 1000000007, the latter finding its primes by trial division. That path calls the full first_primes_after, so the program builds its $20\,000\,000$-entry base sieve and its four-million window twice. Python and Java have no equivalent test.

Complexity and Verification

Let $B$ be the base-prime bound, $\Delta=4\,000\,000$ the window width, $m=100\,000$ the number of primes, $R=\lfloor\sqrt H\rfloor$ and $P$ the number of base primes at most $R$. The base sieve costs $O(B\log\log B)$ time and $B$ bytes: $20$ MB for C++, $11$ MB for Python and Java. The window costs $\Delta$ bytes, one set-up step per base prime and $\sum_{p\le R}\Delta/p$ marking writes, $12\,165\,859$ of them for $P=664\,579$ (Section 5), and the Fibonacci stage $O(m\log H)$, concretely $14\,100\,000$ modular multiplications. In total

$$\begin{aligned} O\bigl(&B\log\log B+\Delta\log\log R\\ &+P+m\log H\bigr), \end{aligned}$$

dominated by the base sieve and the multiplications. Counting a modular multiplication as one operation suits C++ but flatters Java, whose mulMod loops over up to $41$ bits. On one machine, C++ at -O2 took $0.16$ s of processor time, checkpoint sieve included, Java $1.5$ s and Python $2.5$ s.

Each of the following was verified by running code:

Fast doubling matches the plain recurrence for every $n\le3000$ modulo $M$ and $10^{9}+7$; against exact integers, the addition formula holds for $1\le m\le59$, $0\le n\le59$ and the doubling identities for $k\le199$.
The segmented sieve returns exactly the $6174$ primes an independent primality test finds in $[10^{14},10^{14}+200000]$.
The built-in checkpoint reproduces $f(10)=55$ and, modulo $10^{9}+7$, $f(1009)\equiv529241575$, $f(1013)\equiv613716502$, $f(1019)\equiv506824938$ with total $682050181$ over the $30$ primes after $1000$, as exact integers also give.
Section 9's small task gives $633$ by the fast method and from exact integers, and $\pi(M)=900\,788\,112$ satisfies $f(\pi)\equiv0$, $f(\pi+1)\equiv1$ modulo $M$ while no proper divisor of it does.
The three published programs agree with each other, with an independent implementation on a different sieve, and with Section 11's gap stepping; the last two are not packaged in the sources.

Two caveats for non-official inputs. Section 4 needs the base-prime list to reach $\sqrt H$, which --start= and --count= do not check: at --start=400000120000008 the C++ program reports $20000003^{2}$ as the next prime, and the Python function fails likewise past $11000027^{2}$. Second, residue doublings and additions are plain machine words, so a general modulus is safe to $2^{63}$ in unsigned C++ but only to $2^{62}$ in Java, whose mulMod returns $-1$ at $2^{62}+1$ instead of $1$. This modulus leaves a wide margin.

Footnotes and References

The references below supply the background facts; the two proofs, the operation counts, $\pi(M)$ and the gap-stepping method are derived here.

Project Euler 304 — Primonacci. The official statement defines the next prime after n, sets a(1) to the next prime after 10^14 and a(n) to the next prime after a(n−1), defines f by f(0) = 0, f(1) = 1, f(n) = f(n−1) + f(n−2) and b(n) = f(a(n)), and asks for the sum of b(n) over 1 ≤ n ≤ 100000 modulo 1234567891011, with no sample value.
Wikipedia — Sieve of Eratosthenes. Sieving a range in segments, with the multiples of every prime up to the square root of a segment's top value, each prime starting at its square: Sections 4 and 5.
Wikipedia — Fibonacci sequence. Binet's formula and F(n) ~ φ^n/√5 for Section 1, and the identity F(m)F(n+1) + F(m−1)F(n) = F(m+n) proved in Section 6.
Project Nayuki — Fast Fibonacci algorithms. F(2k) = F(k)[2F(k+1) − F(k)] and F(2k+1) = F(k+1)² + F(k)², presented as matrix exponentiation with the redundant work removed, at Θ(log n) operations.
Wikipedia — Pisano period. π(n), the period with which the Fibonacci numbers modulo n repeat, and π(mn) = lcm(π(m), π(n)) for coprime m and n, used in Section 10.
Wikipedia — Prime number theorem. The average gap between consecutive primes among the first N integers is roughly log N, the estimate of Section 3.
Wikipedia — Mertens' theorems. Its second theorem: the sum of 1/p over primes p ≤ n grows like log log n plus a constant, as used in Section 5; the value 3.0414… for p ≤ 10^7 we computed directly.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition, subtraction and multiplication — the licence used in Section 2.
GCC manual — 128-bit integers documents unsigned __int128 on targets with a wide enough integer mode. Java Language Specification §4.2 fixes long at 64 signed bits and states that the integer operators do not indicate overflow, hence mulMod. Python — numeric types gives integers unlimited precision.
C++ — Euler304.cpp. This immutable revision holds sieve_base_primes, first_primes_after, fib_pair_mod, solve, solve_bruteforce and run_checkpoints; the independent implementation and the gap stepping above are extra validation, not tests in it.

Problem 304 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    u64 start = 100000000000000ULL;
    int count = 100000;
    u64 modulo = 1234567891011ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0ULL;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

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;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    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, "--start=", options.start) ||
            parse_int_after_prefix(arg, "--count=", options.count) ||
            parse_u64_after_prefix(arg, "--modulo=", options.modulo)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.start >= 2ULL && options.count >= 1 && options.modulo >= 2ULL;
}

std::vector<int> sieve_base_primes(const int limit) {
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit + 1), 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;
    for (int p = 2; static_cast<u64>(p) * static_cast<u64>(p) <= static_cast<u64>(limit); ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (int q = p * p; q <= limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = 0U;
        }
    }
    std::vector<int> primes;
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(p);
        }
    }
    return primes;
}

std::vector<u64> first_primes_after(u64 start, const int count) {
    const int base_limit = 20000000;  // sqrt(1e14 + margin) is around 1e7.
    const std::vector<int> base_primes = sieve_base_primes(base_limit);

    const u64 segment_size = 4000000ULL;
    std::vector<u64> out;
    out.reserve(static_cast<std::size_t>(count));

    u64 low = start + 1ULL;
    while (static_cast<int>(out.size()) < count) {
        const u64 high = low + segment_size - 1ULL;
        std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(segment_size), 1U);

        for (int p : base_primes) {
            const u64 pp = static_cast<u64>(p);
            if (pp * pp > high) {
                break;
            }
            u64 first = (low + pp - 1ULL) / pp * pp;
            if (first < pp * pp) {
                first = pp * pp;
            }
            for (u64 x = first; x <= high; x += pp) {
                is_prime[static_cast<std::size_t>(x - low)] = 0U;
            }
        }

        for (u64 i = 0; i < segment_size && static_cast<int>(out.size()) < count; ++i) {
            const u64 value = low + i;
            if (value >= 2ULL && is_prime[static_cast<std::size_t>(i)] != 0U) {
                out.push_back(value);
            }
        }
        low = high + 1ULL;
    }
    return out;
}

std::pair<u64, u64> fib_pair_mod(const u64 n, const u64 mod) {
    if (n == 0ULL) {
        return {0ULL, 1ULL % mod};
    }
    const auto [a, b] = fib_pair_mod(n >> 1U, mod);
    const u64 two_b = static_cast<u64>((2ULL * b) % mod);
    const u64 c = static_cast<u64>((static_cast<u128>(a) * ((two_b + mod - a) % mod)) % mod);
    const u64 d = static_cast<u64>((static_cast<u128>(a) * a + static_cast<u128>(b) * b) % mod);
    if ((n & 1ULL) == 0ULL) {
        return {c, d};
    }
    return {d, static_cast<u64>((c + d) % mod)};
}

u64 solve(const u64 start, const int count, const u64 mod) {
    const std::vector<u64> primes = first_primes_after(start, count);
    u64 sum = 0ULL;
    for (u64 p : primes) {
        sum += fib_pair_mod(p, mod).first;
        if (sum >= mod) {
            sum %= mod;
        }
    }
    return sum % mod;
}

u64 solve_bruteforce(const u64 start, const int count, const u64 mod) {
    auto is_prime_trial = [](u64 x) {
        if (x < 2ULL) {
            return false;
        }
        if ((x & 1ULL) == 0ULL) {
            return x == 2ULL;
        }
        for (u64 d = 3ULL; d * d <= x; d += 2ULL) {
            if (x % d == 0ULL) {
                return false;
            }
        }
        return true;
    };

    std::vector<u64> primes;
    for (u64 x = start + 1ULL; static_cast<int>(primes.size()) < count; ++x) {
        if (is_prime_trial(x)) {
            primes.push_back(x);
        }
    }
    u64 sum = 0ULL;
    for (u64 p : primes) {
        sum = (sum + fib_pair_mod(p, mod).first) % mod;
    }
    return sum;
}

bool run_checkpoints() {
    if (fib_pair_mod(10ULL, 1000000007ULL).first != 55ULL) {
        std::cerr << "Checkpoint failed for Fibonacci doubling F(10)" << '\n';
        return false;
    }
    if (solve(1000ULL, 30, 1000000007ULL) != solve_bruteforce(1000ULL, 30, 1000000007ULL)) {
        std::cerr << "Checkpoint failed for brute cross-check on small prime range" << '\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(options.start, options.count, options.modulo) << '\n';
    return 0;
}

Python

def sieve_base_primes(limit):
    is_prime = bytearray(limit + 1)
    for i in range(2, limit + 1):
        is_prime[i] = 1
        
    p = 2
    while p * p <= limit:
        if is_prime[p]:
            for q in range(p * p, limit + 1, p):
                is_prime[q] = 0
        p += 1
        
    primes = [p for p in range(2, limit + 1) if is_prime[p]]
    return primes

def first_primes_after(start, count):
    base_limit = 10000000 + 1000000  # sqrt(10^14) is 10^7
    base_primes = sieve_base_primes(base_limit)
    
    segment_size = 4000000
    out = []
    
    low = start + 1
    while len(out) < count:
        high = low + segment_size - 1
        is_prime = bytearray(segment_size)
        for i in range(segment_size):
            is_prime[i] = 1
            
        for p in base_primes:
            if p * p > high:
                break
                
            first = ((low + p - 1) // p) * p
            if first < p * p:
                first = p * p
                
            start_idx = first - low
            for x in range(start_idx, segment_size, p):
                is_prime[x] = 0
                
        for i in range(segment_size):
            if len(out) >= count:
                break
            if is_prime[i]:
                value = low + i
                if value >= 2:
                    out.append(value)
                    
        low = high + 1
        
    return out

def fib_pair_mod(n, mod):
    if n == 0:
        return (0, 1 % mod)
        
    a, b = fib_pair_mod(n >> 1, mod)
    two_b = (2 * b) % mod
    c = (a * ((two_b + mod - a) % mod)) % mod
    d = (a * a + b * b) % mod
    
    if (n & 1) == 0:
        return (c, d)
    return (d, (c + d) % mod)

def solve(start=100000000000000, count=100000, mod=1234567891011):
    primes = first_primes_after(start, count)
    total_sum = 0
    for p in primes:
        total_sum += fib_pair_mod(p, mod)[0]
        if total_sum >= mod:
            total_sum %= mod
            
    return str(total_sum % mod)

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

Java

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

public class Euler304 {
    static List<Integer> sieveBasePrimes(int limit) {
        byte[] isPrime = new byte[limit + 1];
        for (int i = 2; i <= limit; ++i)
            isPrime[i] = 1;

        for (int p = 2; (long) p * p <= limit; ++p) {
            if (isPrime[p] == 1) {
                for (int q = p * p; q <= limit; q += p) {
                    isPrime[q] = 0;
                }
            }
        }

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

    static List<Long> firstPrimesAfter(long start, int count) {
        int baseLimit = 11000000;
        List<Integer> basePrimes = sieveBasePrimes(baseLimit);

        int segmentSize = 4000000;
        List<Long> out = new ArrayList<>(count);

        long low = start + 1;
        byte[] isPrime = new byte[segmentSize];

        while (out.size() < count) {
            long high = low + segmentSize - 1;
            for (int i = 0; i < segmentSize; ++i)
                isPrime[i] = 1;

            for (int p : basePrimes) {
                long pp = (long) p;
                if (pp * pp > high)
                    break;

                long first = ((low + pp - 1) / pp) * pp;
                if (first < pp * pp) {
                    first = pp * pp;
                }

                int startIdx = (int) (first - low);
                for (int x = startIdx; x < segmentSize; x += p) {
                    isPrime[x] = 0;
                }
            }

            for (int i = 0; i < segmentSize && out.size() < count; ++i) {
                long value = low + i;
                if (value >= 2 && isPrime[i] == 1) {
                    out.add(value);
                }
            }

            low = high + 1;
        }

        return out;
    }

    static class Pair {
        long first, second;

        Pair(long first, long second) {
            this.first = first;
            this.second = second;
        }
    }

    static Pair fibPairMod(long n, long mod) {
        if (n == 0)
            return new Pair(0, 1 % mod);

        Pair half = fibPairMod(n >>> 1, mod);
        long a = half.first;
        long b = half.second;

        long twoB = (2L * b) % mod;

        long c = mulMod(a, (twoB + mod - a) % mod, mod);
        long d = (mulMod(a, a, mod) + mulMod(b, b, mod)) % mod;

        if ((n & 1) == 0) {
            return new Pair(c, d);
        }
        return new Pair(d, (c + d) % mod);
    }

    // multiply modulo handling overflow
    static long mulMod(long a, long b, long mod) {
        long res = 0;
        a %= mod;
        while (b > 0) {
            if ((b & 1) == 1) {
                res = (res + a) % mod;
            }
            a = (a * 2) % mod;
            b >>>= 1;
        }
        return res;
    }

    public static String solve() {
        long start = 100000000000000L;
        int count = 100000;
        long mod = 1234567891011L;

        List<Long> primes = firstPrimesAfter(start, count);
        long sum = 0;
        for (long p : primes) {
            sum += fibPairMod(p, mod).first;
            if (sum >= mod)
                sum %= mod;
        }

        return String.valueOf(sum % mod);
    }

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