Generate Fibonacci Primes Quickly

Question

Unsurprisingly, fibonacci primes are primes that are also Fibonacci numbers. There are currently 34 known Fibonacci primes and an additional 15 probable Fibonacci primes. For the purpose of this challenge, the Fibonacci numbers are the sequence \$F_n\$ defined as \$F_1 = 1\$, \$F_2 = 1\$, and \$F_n = F_{n-1} + F_{n-2}\$, and a number is considered prime if it passes a probabilistic prime test with a probability of being incorrect of less than \$2^{-32}\$. For example, since a \$k\$ round Miller-Rabin test has an error probability of \$4^{-k}\$, a 16 round Miller-Rabin test is sufficient to prove primality for the purpose of this challenge.

Submissions:

The goal of this challenge is to write a full program that calculates every Fibonacci prime and its index in the Fibonacci series as fast as possible.

Submissions shall be a full program, complete with instructions for building and running it. Submissions must be in a language freely available for noncommercial use and capable of running on Windows 10, and users must be prepared to provide installation instructions for that language. External libraries are permitted, with the same caveats that apply to languages.

Primes will be output by writing them as base 10 integers to stdout in ascending order in the format

index,prime\n

where \n is the newline character. The numbers can have extraneous leading/trailing whitespace other than the newline character.

Scoring

The programs will be run on an Intel(R) Core(TM) i5-8365U CPU with 8 threads, avx-2 support, and 24 Gigabytes of ram. The largest prime that can be correctly reached in one minute wins. Tiebreaker is the time taken to reach the largest value. Programs that tamper with my computer or the testing program will be disqualified. Programs that error or otherwise fail to produce the correct output will be judged based on the furthest Fibonacci prime reached before they failed.

Results

Sorry for the delay, getting Anders Kaseorg's entry to build properly wound up being way harder than it needed to be.

Results:
  anders_kaseorg:
    num: 27
    time: 44.6317962s
  aleph_alpha:
    num: 24
    time: 22.6188601s
  example:
    num: 12
    time: 23.2418081s

see also: A005478, A001605

The test program can be found here. Additionally, there is an example program here.

Answer 1

Rust, \$F_{14431}\$ in ≈ 58 s

(Unofficial estimate based on the expected CPU performance ratio.)

After special-casing \$F_3, F_4\$, this tests every prime-indexed Fibonacci number with the GMP mpz_probab_prime_p function, which performs a Baillie–PSW test and rep - 24 Miller–Rabin iterations. Since Baillie-PSW has no proven error bound, I pass 24 + 16 to ensure 16 Miller–Rabin iterations.

Build with cargo build --release and run target/release/fibprimes.

Cargo.toml

[package]
name = "fibprimes"
version = "0.1.0"
edition = "2021"

[dependencies]
pariter = "0.5.1"
rug = { version = "1.16.0", features = ["integer"], default-features = false }

src/main.rs

use pariter::IteratorExt;
use rug::{integer::IsPrime, Integer};
use std::iter;

fn main() {
    println!("3,2");
    println!("4,3");

    let mut fib0 = Integer::from(1);
    let mut fib1 = Integer::from(2);
    let mut n = Integer::from(3);

    for (n, fib) in iter::from_fn(move || {
        let n1 = Integer::from(n.next_prime_ref());
        while n != n1 {
            fib0 += &fib1;
            fib1 += &fib0;
            n += 2;
        }
        Some((n.clone(), fib1.clone()))
    })
    .parallel_filter(|(_, fib)| fib.is_probably_prime(24 + 16) != IsPrime::No)
    {
        println!("{n},{fib}");
    }
}

Answer 2

C + PARI/GP's C library, \$F_{9677}\$ in 25s on my computer

#include <pari/pari.h>

void f(long n)
{
    GEN a = fibo(n);
    if (ispseudoprime(a, 16))
        pari_printf("%d,%Ps\n", n, a);
    return;
}

int main()
{
    pari_init(8000000, 500000);

    f(3);
    f(4);

    long n;
    forprime_t iter;
    u_forprime_init(&iter, 5, ULONG_MAX);

    while (n = u_forprime_next(&iter))
        f(n);

    return 0;
}

If you are using Ubuntu or Debian, you need to install the package libpari-dev. For Arch Linux, pari-gp is enough.

Compiles with gcc -O3 filename.c -lpari.

Answer 3

Wolfram Language (Mathematica), F₉₆₇₇ in <1 s

PrimeQ[(MatrixPower[{{1,1},{1,0}},#].{1,0})[[2]]]&

Try it online!

If you are looking for fast code, then I strongly suggest using a closed form expression for the Fibonacci numbers, not a recursive definition. For the code above, I used the definition $$\begin{bmatrix}F_{n+1}\\F_n\end{bmatrix} = \begin{bmatrix}1&1\\1&0\end{bmatrix}^n\begin{bmatrix}1\\0\end{bmatrix}.$$ The second entry in the resulting matrix is the nth Fibonacci number, with the sequence starting at $$F_1 = F_2 = 1.$$

I suspect that this will be a lot faster to implement, especially if you use a language optimized for matrix algebra.

After generating the nth Fibonacci number, I checked whether it was prime using a Mathematica function PrimeQ.

Aside: While we are on the subject, there is another closed-form expression for the nth Fibonacci number. This might be helpful in a different context. $$ F_n = \frac{(1+\sqrt{5})^n - (1-\sqrt{5})^n}{2^n\sqrt{5}}.$$ The indexing is the same as the matrix formula above. For generating a list of Fibonacci numbers, this formula is very fast, but for n>75, there was roundoff error caused by the fact that the square root of 5 is stored as a floting point real. Nevertheless, if you have a language that handles numbers with great precision, then this formula might be best for you.

As for your restrictions on testing and scoring this program: you'd have to use WolframScript to run this from the command line in Windows 10. Alternatively you can just use the tio.run link here (or above) to see the first 9677 Fibonacci numbers tested for prime-ness. Wouldn't tio.run be a more fair test of speed?