Problem 399: Squarefree Fibonacci Numbers

Project Euler Problem 399 Solution

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

Problem Summary Let $F_1=F_2=1$ and $F_{n+2}=F_{n+1}+F_n$. The task is to locate the index $n$ of the $N$-th Fibonacci number that is squarefree, meaning that no prime square divides $F_n$. After that index is known, the program prints $F_n$ in a compressed decimal format: the last 16 digits and a scientific-notation leading value with one decimal place. The crucial observation is that the search is carried out on the index line $n=1,2,3,\dots$, not on the enormous Fibonacci values themselves. Mathematical Approach 1. Move the Squarefree Test to the Index Line A positive integer is squarefree exactly when it is not divisible by $p^2$ for any prime $p$. Therefore $$\forall p,\ p^2 \nmid F_n.$$ So instead of examining the factorization of $F_n$ directly, we ask a sharper question: for a fixed prime $p$, which indices $n$ force $p^2\mid F_n$? Once those forbidden indices are understood, the problem becomes a sieve over the natural numbers. 2....

Detailed mathematical approach

Problem Summary

The Fibonacci numbers are $F_1=F_2=1$ and $F_{n+2}=F_{n+1}+F_n$, so the first fifteen of them are

$$\begin{gathered}1,1,2,3,5,8,13,21,34,\\55,89,144,233,377,610.\end{gathered}$$

A positive integer is squarefree when no prime square divides it. Two of those fifteen fail: $8=2^3$ is divisible by $4$, and $144=2^4\cdot3^2$ by $4$ and by $9$. Dropping them leaves the first thirteen squarefree Fibonacci numbers, $1,1,2,3,5,13,21,34,55,89,233,377$ and $610$ — the two copies of $1$ counting separately, so the task is to count indices $n$ with $F_n$ squarefree.

We must find the $100\,000\,000$-th of them and report its last sixteen digits, then a comma, then the number in scientific notation rounded to one digit after the decimal point. The statement prints the $200$-th squarefree Fibonacci number in full, a $54$-digit number whose required output is 1608739584170445,9.7e53. It also tells us to assume part of Wall's conjecture: for every prime $p$, the first Fibonacci number divisible by $p$ is not divisible by $p^2$. That holds for every $p\le3\cdot10^{15}$ but is unproved, and were it false the computed value would only be a lower bound. ^[1]

The whole solution in one sentence: for each prime $p$ the Fibonacci numbers divisible by $p^2$ are exactly those whose index is a multiple of one fixed number $b_p$, so we sieve the indices $1,2,3,\dots$ with the $b_p$, take the $100\,000\,000$-th survivor, and read the two requested pieces off a modular computation and a logarithm.

Mathematical Approach

1. Why the index, and not the number

Factorising the target is hopeless: $F_n$ has about $0.209\,n$ decimal digits (Section 10), so the number asked for has tens of millions of them. Turn the question around instead. Squarefreeness fails exactly when some $p^2$ divides $F_n$, so for each prime $p$ we ask which indices $n$ make $p^2\mid F_n$. The answer is the multiples of a single number, which makes the problem a sieve over the index line.

2. The entry point of a prime

For a prime $p$, the entry point (or rank of apparition) is

$$z(p)=\min\{m\ge1:\ p\mid F_m\}.$$

Running the recurrence modulo $p$ gives $z(2)=3$, $z(3)=4$, $z(5)=5$, $z(7)=8$, $z(11)=10$, $z(13)=7$, $z(17)=9$, $z(19)=18$. ^[3] Since $F_1=F_2=1$, every prime has $z(p)\ge3$, with equality only at $p=2$. The classical fact we need is

$$p\mid F_n\iff z(p)\mid n.$$

One direction follows from Section 3. For the other use $\gcd(F_a,F_b)=F_{\gcd(a,b)}$: ^[2] if $p\mid F_n$ then $p$ divides $\gcd(F_n,F_{z(p)})=F_{\gcd(n,z(p))}$, and since $\gcd(n,z(p))\le z(p)$, minimality forces $\gcd(n,z(p))=z(p)$.

3. One identity that does all the work

Set $F_0=0$, let $I$ be the identity matrix and let $A$ have rows $(1,1)$ and $(1,0)$. For $n\ge1$ the rows of $A^n$ are $(F_{n+1},F_n)$ and $(F_n,F_{n-1})$: true at $n=1$, and multiplying by $A$ advances both rows by the recurrence. ^[2] Hence $A^m=F_m\,A+F_{m-1}\,I$ for $m\ge1$. The two summands commute, so the binomial theorem expands the $k$-th power, its term $j$ carrying the scalar $\binom{k}{j}F_m^{\,j}F_{m-1}^{\,k-j}$. The upper-right entry of $A^j$ is $F_j$ for $j\ge1$ and $0$ for $j=0$, so reading off that entry gives, for $m,k\ge1$,

$$F_{km}=\sum_{j=1}^{k}\binom{k}{j}F_m^{\,j}F_{m-1}^{\,k-j}F_j.$$

Two consequences matter. Every term has a factor $F_m$, so $F_m\mid F_{km}$: index divisibility carries over to the numbers. And pulling out one $F_m$,

$$F_{km}=F_m\bigl(k\,F_{m-1}^{\,k-1}+F_m\,T\bigr)$$

for some integer $T$, since every term with $j\ge2$ keeps a spare factor $F_m$.

4. The key claim

Consecutive Fibonacci numbers are coprime, $\gcd(F_m,F_{m-1})=1$, because the recurrence turns any common divisor into a divisor of $F_1=1$; so $p\mid F_m$ implies $p\nmid F_{m-1}$. Write $v_p(x)$ for the exponent of $p$ in $x$ and put $z=z(p)$. Wall's assumption says exactly $v_p(F_z)=1$. Then, for every prime $p$,

$$p^2\mid F_n\iff p\,z(p)\mid n.$$

Proof. Let $p^2\mid F_n$. Then $p\mid F_n$, so $n=kz$ by Section 2, and the second identity of Section 3 gives $F_n=F_z(k\,F_{z-1}^{\,k-1}+F_zT)$. If $p\nmid k$ the bracket is nonzero modulo $p$, since $p\nmid F_{z-1}$ while $p\mid F_z$, so $v_p(F_n)=v_p(F_z)=1$ — a contradiction. Hence $pz\mid n$. Conversely take $k=p$, $m=z$ in the sum: its term $j=1$ is $p\,F_z\,F_{z-1}^{\,p-1}$, divisible by $p^2$; for $2\le j\le p-1$ the coefficient $\binom{p}{j}$ is divisible by $p$ and $F_z^{\,j}$ by $p^2$; for $j=p$, $F_z^{\,p}$ alone is. So $p^2\mid F_{pz}$, and $pz\mid n$ gives $F_{pz}\mid F_n$, hence $p^2\mid F_n$. □

So the first Fibonacci number divisible by $p^2$ sits at index $b_p=p\,z(p)$, and the bad indices for $p$ are precisely the multiples of $b_p$. ^[5]

5. A worked example up to index 60

The first values are $b_2=2\cdot3=6$, $b_3=3\cdot4=12$, $b_5=25$, $b_7=7\cdot8=56$; the next, $b_{13}=13\cdot7=91$, already exceeds $60$. Striking out their multiples up to $60$ removes

$$\begin{gathered}6,12,18,24,25,30,36,\\42,48,50,54,56,60,\end{gathered}$$

thirteen indices. Factorising $F_1,\dots,F_{60}$ completely gives exactly the same thirteen, the first thirteen terms of the published list of indices whose Fibonacci number has a square factor. ^[6] For $n\le15$ only $6$ and $12$ are struck out, reproducing the thirteen values quoted in the statement.

6. Removing redundant periods

Here $12$ is wasted work: $6\mid12$, so its multiples are already multiples of $6$. In general, if $b\mid b'$ then $b'$ may be deleted. Sorting the $b_p$ and keeping only those that are not multiples of a smaller kept one leaves

$$\begin{gathered}6,25,56,91,110,153,\\406,703,752,820,915,979,\dots\end{gathered}$$

which is precisely the published sequence of indices $n$ at which $F_n$ has a square factor but no $F_k$ does for a proper divisor $k$ of $n$. ^[5] Note also that inside a window $1\le n\le K$ only primes $p\le K/3$ can matter: since $z(p)\ge3$, a larger prime has $b_p\ge3p\gt K$.

7. How wide the window must be

The program searches a fixed window $1\le n\le K$ with $K=220\,000\,000$. Is that enough for the $100\,000\,000$-th survivor? There are $2\,994$ distinct values $b_p\le K$, from primes up to $3\,010\,349$ (entry point only $62$), and $2\,139$ remain after the reduction above. At most $\sum_p\lfloor K/b_p\rfloor=78\,554\,623$ indices are struck out, so at least $141\,445\,377$ survive — already more than enough. Counting exactly gives $168\,227\,174$ survivors, a density of $0.764669$. As a cross-check, the reciprocal of the bad-index proportion is $4.2493$, putting the $n$-th bad index near $4.2493\,n$ in this window; the published entry asks whether such a constant exists between $4$ and $4.5$ and argues that, if it does, it is below $4.3$. ^[6]

8. Computing $z(p)$ for millions of primes

Running the recurrence until $p$ appears costs $O(p)$ steps per prime — far too slow. Instead use the classical congruence $p\mid F_{p-(5/p)}$, where the Legendre symbol $(5/p)$ is $+1$ when $5$ is a square modulo $p$ and $-1$ when it is not, ^[2] which with Section 2 says $z(p)\mid r$ for

$$r=p-\left(\tfrac{5}{p}\right),\qquad p\ne2,5.$$

Euler's criterion evaluates the symbol by one modular power, ^[7]

$$5^{(p-1)/2}\equiv\left(\tfrac{5}{p}\right)\pmod p.$$

Now shrink $r$ to $z(p)$: for each prime $q$ dividing $r$, repeatedly test whether $F_{r/q}\equiv0\pmod p$ and, if so, replace $r$ by $r/q$. The divisors $d$ of $r$ with $p\mid F_d$ are exactly the multiples of $z(p)$ dividing $r$, so dropping a factor $q$ stays inside that set precisely while $v_q(r)-1\ge v_q(z(p))$; the loop thus lowers $v_q(r)$ to $v_q(z(p))$ and stops, in any order of the $q$. Take $p=13$: Euler's criterion gives $5^6\equiv-1\pmod{13}$, so $r=14=2\cdot7$. Dividing by $2$ succeeds, since $F_7=13$, leaving $r=7$; dividing by $7$ fails, since $F_1=1$; so $z(13)=7$. The exceptions are set directly: $z(2)=3$, $z(5)=5$.

9. Fast doubling

Every Fibonacci residue above is produced by fast doubling: from the pair $(F_k,F_{k+1})$,

$$F_{2k}=F_k\,(2F_{k+1}-F_k),$$

$$F_{2k+1}=F_k^2+F_{k+1}^2,$$

and $F_{2k+2}=F_{2k}+F_{2k+1}$. ^[8] Both come from squaring the matrix $A^k$ of Section 3: the upper-right entry of $A^{2k}$ is $F_k(F_{k+1}+F_{k-1})=F_k(2F_{k+1}-F_k)$ and its upper-left entry is $F_{k+1}^2+F_k^2$. For instance $(F_3,F_4)=(2,3)$ gives $F_6=2(6-2)=8$ and $F_7=4+9=13$. Reading the bits of $n$ from the top, one pair update per bit computes $F_n$ modulo any $M$ in $O(\log n)$ modular multiplications.

10. The two pieces of the answer

The last sixteen digits are simply $F_n\bmod10^{16}$, one fast-doubling run with $M=10^{16}$. For the scientific notation use Binet's formula $F_n=(\varphi^n-\psi^n)/\sqrt5$ with $\varphi=(1+\sqrt5)/2$ and $\psi=(1-\sqrt5)/2$. ^[2] Since $\psi/\varphi=-\varphi^{-2}$,

$$\log_{10}F_n=n\log_{10}\varphi-\log_{10}\sqrt5+\varepsilon_n$$

with $\varepsilon_n=\log_{10}(1-(-\varphi^{-2})^{n})$. Put $t=\varphi^{-2n}$: then $|\ln(1\pm t)|\le-\ln(1-t)\le t/(1-t)$, and dividing by $\ln10$ gives $|\varepsilon_n|\lt t$, because $(1-\varphi^{-2})\ln10=1.42\gt1$ bounds $1/((1-t)\ln10)$ by $1$ for $n\ge1$. At the index we need $t$ is below $10^{-10^{7}}$, so the term is dropped. Writing $x$ for the result — note $\log_{10}\varphi=0.2089876\ldots$, the digit count of Section 1 — the exponent is $e=\lfloor x\rfloor$ and the mantissa $10^{\,x-e}\in[1,10)$ is rounded to one decimal; should that give $10.0$, it becomes $1.0$ and $e$ grows by one. For $n=260$ this returns $9.7\mathrm{e}53$ and $F_{260}\bmod10^{16}=1608739584170445$, the statement's output — and $260$ is indeed the $200$-th surviving index.

How the Code Works

The C++ program ^[10] reads --target= (default $100\,000\,000$), --k-upper= (default $220\,000\,000$) and --skip-checkpoints. parse_arguments rejects any other argument and demands a target of at least $1$ and a window of at least $1000$, but its std::stoll never checks that the whole value was consumed, so --target=200garbage passes as $200$. Unless skipped, run_checkpoints solves target $200$ in a window of $2000$ and aborts unless the result is 1608739584170445,9.7e53.

periods_up_to_k implements Sections 6 and 8. It sets the prime limit to $K/3$ and calls build_spf, a linear sieve recording the smallest prime factor of every integer up to that limit and collecting the primes in one $O(K)$ pass. ^[9] Two spare entries leave room for the index $r=p+1$; at the default bound the largest prime in the table is $73\,333\,319$. The prime list is cut into contiguous blocks, each building its own period list; every prime contributes once and the lists are merged and sorted afterwards, so the result is independent of how the blocks are cut.

Per prime, rank_of_apparition_prime returns $3$ or $5$ for the exceptions, otherwise calls legendre_5_mod_p — Euler's criterion through mod_pow_u64 — forms $r$, factors it with unique_prime_factors off the sieve, and runs the shrinking loop of Section 8 through fibonacci_mod. That wraps fib_pair_mod, a recursive fast doubling returning $(F_n,F_{n+1})$ modulo $m$ at depth $O(\log n)$, with every product widened to __uint128_t before reduction, so a modulus of $10^{16}$ cannot overflow. reduce_periods then sorts, deduplicates and drops each period that is a multiple of a kept smaller one.

nth_squarefree_fib_index allocates one byte per index up to $K$, marks the multiples of every reduced period and returns the index where the survivor count reaches the target, or $-1$ so that main reports too small a window. last16_and_scientific does the two computations of Section 10: fibonacci_mod with modulus $10^{16}$, then Boost's cpp_dec_float_50, $50$ decimal digits plus guard digits, ^[4] for $\sqrt5$, $\varphi$ and the logarithms — a vast margin here, though no error interval is computed. floor gives the exponent, floor(mant*10+0.5)/10 the rounded mantissa, renormalised if it reaches $10$, and setw(16) with setfill('0') keeps a leading zero in the digit block. Dropping $\varepsilon_n$ unconditionally costs something at the bottom of the range: with --target=1 --k-upper=1000 the program prints 0000000000000001,7.2e-1 where $F_1=1$ calls for 1.0e0, and targets $2$, $3$, $4$ fail likewise; from $5$ up it is right.

The Python and Java files hold no mathematics: each wraps the same C++ source. Python's solve rebuilds solutionsCpp/Euler399.cpp at -O2 -std=c++17 when the cached binary is missing or stale, runs it with no option but --skip-checkpoints, and takes the last non-empty output line unless an earlier one matches answer: or contains =; Java's Euler399 mirrors this, but finds the project root in the current working directory rather than beside its own file, so it must be started from that root. Both pass --skip-checkpoints whenever the C++ source text mentions that flag, which it always does, so the checkpoint runs only when the C++ binary is invoked directly.

Complexity and Verification

The linear sieve runs in $O(K)$ time and yields $4\,302\,429$ primes below $K/3=73\,333\,333$. The entry points then dominate: each prime needs $O(\log p)$ factorisation lookups and $O(\log p)$ Fibonacci tests of $O(\log p)$ steps each, so $O(\pi(K)\log^2K)$ word operations. Measured, the shrinking loops make $18\,239\,259$ tests, $4.24$ per prime, of $18.74$ doubling steps on average: $1\,025\,259\,738$ products there and $156\,829\,048$ in the Legendre powers, about $1.2\cdot10^{9}$ modular multiplications. Reducing the $D=2\,994$ collected periods costs $O(D\log D)$ to sort plus $O(D^2)$ divisibility tests, and marking takes $\sum_b\lfloor K/b\rfloor=58\,581\,698$ writes with the reduced set against $78\,554\,623$ without it, then a scan of at most $K$ indices.

Memory is $O(K)$ and peaks in the prime phase, not at the end: the smallest-prime-factor table takes about $293$ MB and the prime list $17$ MB, both released before the $220\,000\,001$-byte mark array is allocated. Here the peak resident set was $314$ MB and the run took about one second — both machine-dependent.

The program checks itself on the one case the statement supplies; the following were verified independently by computation for this article:

the identity of Section 3 for $1\le m,k\le24$, $\gcd(F_m,F_{m-1})=1$ for $m\le50$, $F_a\mid F_b$ whenever $a\mid b$ for $a,b\le60$, and the doubling identities for $k\lt190$;
$v_p(F_{z(p)})=1$ — Wall's assumption — for every prime $p\lt2000$, and $z(p)\mid p-(5/p)$ for every prime $2\lt p\lt20\,000$, $p\ne5$;
the equivalence of Section 4 for every prime $p\lt120$ and $n\le600$, and that the struck-out set for $n\le200$ matches full factorisation of $F_1,\dots,F_{200}$: $46$ bad indices, $154$ survivors;
the bound $|\varepsilon_n|\lt\varphi^{-2n}$ for $n\le2000$, and that the Binet route gives the same rounded mantissa and exponent as exact integer evaluation for every $5\le n\le399$, failing only for $n\le4$;
the window of Section 7: a separate program recomputed the $2\,994$ periods for $K=220\,000\,000$ and counted $168\,227\,174$ survivors.

The published C++ program invoked with --target=200 --k-upper=2000 prints 1608739584170445,9.7e53.

One caveat remains: the sieve is exact given Wall's assumption. Without it every marked index is still bad, but bad indices may be missed, so the answer can only be too small or right — a lower bound, as the statement says, equality still possible if no extra exclusion falls among the first $100\,000\,000$ survivors. The same gap makes the published sequence of Section 6 conjectural beyond its twelfth term: one prime square refutes squarefreeness at once, but certifying that a large $F_n$ is squarefree needs its factorisation. ^[5]

Footnotes and References

The references supply the background facts used above; the identity of Section 3, the equivalence of Section 4 with its proof, the window estimate of Section 7 and the operation counts are derived or measured here.

Project Euler 399 — Squarefree Fibonacci Numbers. Prints the first fifteen Fibonacci numbers and the 200th squarefree one with its output 1608739584170445,9.7e53, asks for the 100,000,000th in the same format, and says to assume — checked up to 3 · 10^15 — that the first Fibonacci number divisible by a prime p is never divisible by p squared.
Wikipedia — Fibonacci sequence. Gives gcd(F_a, F_b) = F_gcd(a, b) and that F_kn is divisible by F_n, the congruence p | F_{p − (5/p)} with the Legendre symbol, the matrix A with rows (1, 1) and (1, 0) whose powers have entries F_{n+1}, F_n, F_n, F_{n−1}, and Binet's formula.
OEIS A001177 — Fibonacci entry points. Here a(n) is the least k ≥ 1 with n | F_k, beginning 1, 3, 4, 6, 5, 12, 8, and the comments state that all solutions of F_m ≡ 0 (mod n) are the multiples of a(n) — the equivalence of Section 2. The companion A001602 gives the same over the primes: 3, 4, 5, 8, 10, 7, 9, 18.
Boost.Multiprecision — cpp_dec_float. A header-only radix-10 software float; cpp_dec_float_50 carries 50 decimal digits plus internal guard digits and supplies the square root and logarithms of Section 10.
OEIS A065069. Numbers n with F_n not squarefree while F_k is squarefree for every proper divisor k of n — the reduced list of Section 6: 6, 25, 56, 91, 110, 153, 406, 703, 752, 820, 915, 979, 1431, 1892. It states that the first Fibonacci number divisible by the square of the n-th prime has index A001602(n) times that prime, drops products that are multiples of an earlier term such as 12, and calls later terms conjectural pending factorisations.
OEIS A037917. Numbers n with F_n divisible by a square: 6, 12, 18, 24, 25, 30, 36, 42, 48, 50, 54, 56, 60, … — the thirteen indices of Section 5. These are the positive multiples of A065069; the entry asks whether a(n) ~ C·n with 4 < C < 4.5, and argues that such a C would be below 4.3.
Wikipedia — Euler's criterion. For an odd prime p and a coprime to p, (a/p) ≡ a^((p−1)/2) (mod p), the power being 1 for a quadratic residue and −1 otherwise — the single exponentiation that supplies (5/p).
Nayuki — Fast Fibonacci algorithms. Gives the doubling identities F(2k) = F(k)·(2F(k+1) − F(k)) and F(2k+1) = F(k+1)² + F(k)², needing Θ(log n) arithmetic operations.
cp-algorithms — Linear sieve. Records the minimum prime factor of every integer up to n in O(n) time, each composite written once because the inner loop multiplies i only by primes up to the minimum prime factor of i — the sieve behind build_spf.
C++ — Euler399.cpp. The linked immutable revision contains periods_up_to_k, rank_of_apparition_prime, reduce_periods, nth_squarefree_fib_index and last16_and_scientific with the options above; it is the implementation record for this article, while the measurements above come from separate instrumented programs.

Problem 399 source code

C++

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <mutex>
#include <sstream>
#include <string>
#include <thread>
#include <vector>

#include <boost/multiprecision/cpp_dec_float.hpp>

namespace {

using u64 = std::uint64_t;
using i64 = long long;
using boost::multiprecision::cpp_dec_float_50;

constexpr u64 kLast16Mod = 10000000000000000ULL;

struct Options {
    i64 target = 100000000;
    i64 k_upper = 220000000;
    bool run_checkpoints = true;
};

bool parse_i64_after_prefix(const std::string& arg, const std::string& prefix, i64& 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::stoll(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_i64_after_prefix(arg, "--target=", options.target) ||
            parse_i64_after_prefix(arg, "--k-upper=", options.k_upper)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.target >= 1 && options.k_upper >= 1000;
}

u64 mod_pow_u64(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0) {
        if (e & 1ULL) {
            result = static_cast<u64>((__uint128_t)result * cur % mod);
        }
        cur = static_cast<u64>((__uint128_t)cur * cur % mod);
        e >>= 1ULL;
    }
    return result;
}

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

u64 fibonacci_mod(u64 n, u64 mod) {
    return fib_pair_mod(n, mod).first;
}

int legendre_5_mod_p(const int p) {
    if (p == 5) {
        return 0;
    }
    const u64 ls = mod_pow_u64(5ULL, static_cast<u64>((p - 1) / 2), static_cast<u64>(p));
    if (ls == 1ULL) {
        return 1;
    }
    if (ls == static_cast<u64>(p - 1)) {
        return -1;
    }
    return 0;
}

std::vector<int> build_spf(const int limit, std::vector<int>& primes) {
    std::vector<int> spf(static_cast<std::size_t>(limit + 1), 0);
    primes.clear();
    primes.reserve(limit / 10);
    for (int i = 2; i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (const int p : primes) {
            const i64 x = 1LL * i * p;
            if (x > limit) {
                break;
            }
            spf[static_cast<std::size_t>(x)] = p;
            if (p == spf[static_cast<std::size_t>(i)]) {
                break;
            }
        }
    }
    return spf;
}

std::vector<int> unique_prime_factors(int x, const std::vector<int>& spf) {
    std::vector<int> factors;
    int n = x;
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        factors.push_back(p);
        while (n % p == 0) {
            n /= p;
        }
    }
    return factors;
}

int rank_of_apparition_prime(const int p, const std::vector<int>& spf) {
    if (p == 2) {
        return 3;
    }
    if (p == 5) {
        return 5;
    }

    const int leg = legendre_5_mod_p(p);
    int r = p - leg;
    const std::vector<int> factors = unique_prime_factors(r, spf);
    for (const int q : factors) {
        while (r % q == 0) {
            const int candidate = r / q;
            if (fibonacci_mod(static_cast<u64>(candidate), static_cast<u64>(p)) == 0ULL) {
                r = candidate;
            } else {
                break;
            }
        }
    }
    return r;
}

std::vector<int> reduce_periods(std::vector<int> periods) {
    std::sort(periods.begin(), periods.end());
    periods.erase(std::unique(periods.begin(), periods.end()), periods.end());

    std::vector<int> reduced;
    reduced.reserve(periods.size());
    for (const int b : periods) {
        bool redundant = false;
        for (const int kept : reduced) {
            if (kept > b) {
                break;
            }
            if (b % kept == 0) {
                redundant = true;
                break;
            }
        }
        if (!redundant) {
            reduced.push_back(b);
        }
    }
    return reduced;
}

std::vector<int> periods_up_to_k(const i64 k_upper) {
    const int p_limit = static_cast<int>(k_upper / 3);
    std::vector<int> primes;
    const std::vector<int> spf = build_spf(p_limit + 2, primes);

    const unsigned hw = std::thread::hardware_concurrency();
    const int threads = std::max(1U, hw == 0 ? 4U : hw);
    std::vector<std::vector<int>> locals(static_cast<std::size_t>(threads));
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    auto worker = [&](const int tid) {
        const std::size_t total = primes.size();
        const std::size_t l = total * static_cast<std::size_t>(tid) / static_cast<std::size_t>(threads);
        const std::size_t r = total * static_cast<std::size_t>(tid + 1) / static_cast<std::size_t>(threads);

        std::vector<int>& out = locals[static_cast<std::size_t>(tid)];
        out.reserve((r - l) / 20 + 16);

        for (std::size_t i = l; i < r; ++i) {
            const int p = primes[i];
            const int z = rank_of_apparition_prime(p, spf);
            const i64 period = static_cast<i64>(p) * static_cast<i64>(z);
            if (period <= k_upper) {
                out.push_back(static_cast<int>(period));
            }
        }
    };

    for (int t = 0; t < threads; ++t) {
        workers.emplace_back(worker, t);
    }
    for (std::thread& th : workers) {
        th.join();
    }

    std::vector<int> periods;
    for (const auto& v : locals) {
        periods.insert(periods.end(), v.begin(), v.end());
    }
    return reduce_periods(std::move(periods));
}

i64 nth_squarefree_fib_index(const i64 target, const i64 k_upper, const std::vector<int>& periods) {
    std::vector<std::uint8_t> bad(static_cast<std::size_t>(k_upper + 1), 0U);
    for (const int b : periods) {
        for (i64 x = b; x <= k_upper; x += b) {
            bad[static_cast<std::size_t>(x)] = 1U;
        }
    }

    i64 count = 0;
    for (i64 i = 1; i <= k_upper; ++i) {
        if (bad[static_cast<std::size_t>(i)] == 0U) {
            ++count;
            if (count == target) {
                return i;
            }
        }
    }
    return -1;
}

std::string last16_and_scientific(const i64 index) {
    const u64 last16 = fibonacci_mod(static_cast<u64>(index), kLast16Mod);

    const cpp_dec_float_50 sqrt5 = sqrt(cpp_dec_float_50(5));
    const cpp_dec_float_50 phi = (cpp_dec_float_50(1) + sqrt5) / cpp_dec_float_50(2);
    const cpp_dec_float_50 lg = cpp_dec_float_50(index) * log10(phi) - log10(sqrt5);

    long long exponent = floor(lg).convert_to<long long>();
    cpp_dec_float_50 frac = lg - cpp_dec_float_50(exponent);
    cpp_dec_float_50 mant = pow(cpp_dec_float_50(10), frac);
    cpp_dec_float_50 rounded = floor(mant * 10 + cpp_dec_float_50("0.5")) / 10;
    if (rounded >= 10) {
        rounded /= 10;
        ++exponent;
    }

    std::ostringstream oss;
    oss << std::setfill('0') << std::setw(16) << last16 << ','
        << std::fixed << std::setprecision(1) << rounded << 'e' << exponent;
    return oss.str();
}

bool run_checkpoints() {
    const i64 sample_target = 200;
    const i64 sample_upper = 2000;
    const std::vector<int> sample_periods = periods_up_to_k(sample_upper);
    const i64 idx = nth_squarefree_fib_index(sample_target, sample_upper, sample_periods);
    if (idx <= 0) {
        std::cerr << "Checkpoint failed: could not locate 200th index\n";
        return false;
    }

    const std::string sample = last16_and_scientific(idx);
    if (sample != "1608739584170445,9.7e53") {
        std::cerr << "Checkpoint failed for 200th sample, got " << sample << '\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::vector<int> periods = periods_up_to_k(options.k_upper);
    i64 index = nth_squarefree_fib_index(options.target, options.k_upper, periods);

    if (index < 0) {
        std::cerr << "Upper bound too small, no index found up to " << options.k_upper << '\n';
        return 3;
    }

    std::cout << last16_and_scientific(index) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler399 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler399.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler399_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler399 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler399 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler399 C++ bridge produced empty output.");
        }
        return parsed;
    }

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