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.

• Arnauld: It may not be worth the effort, but \$F_n\bmod 6\$ has a period of \24ドル\$ and only half of them are \1ドル\$ or \5ドル\$ (i.e. prime candidates).
• Arnauld: (But any primality test is likely to reject the other ones immediately anyway.)
• Anders Kaseorg: @Arnauld Those are already rejected by testing the index (\$n\$) for primality.
• Aiden4: @AndersKaseorg I changed the output to a simple text one. Do note that this invalidates your answer (not that it's a hard fix). Also, I definitely meant in ascending order, good catch.

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}");
 }
}

• Aiden4: I just ran the tests, and your code came up with 27 primes in 44 seconds. Just out of curiosity, is there a reason you are using the fairly obscure pariter crate instead of rayon?
• Anders Kaseorg: @Aiden4 While Rayon is great and I’ve used it on several challenges here, it doesn’t currently support parallelizing infinite iterators or extracting the results sequentially. This is a fairly obscure problem, but pariter is exactly the right tool for it.

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.

• Aiden4: I just ran the test, your program came up with 24 primes in 22 seconds. I think there is some sort of buffering issue with your program because all of the primes were printed simultaneously.

Answer 3

Python, gmpy2, \$F_{9677}\$ in 27s on my computer.

Same trick as @Anders Kaseorg's Rust answer.

import time
import gmpy2
from gmpy2 import mpz
# from itertools import count
import time
# set the precision for probable prime test
gmpy2.get_context().precision = 24 + 16
def next_prime(n):
 candidate = n + 1 if n % 2 == 0 else n + 2
 while not gmpy2.is_prime(candidate):
 candidate += 2
 return candidate
def fib_prime_gen():
 fib0, fib1 = mpz(3), mpz(5)
 n = mpz(5)
 yield n, fib1
 while True:
 fib0, fib1 = fib0 + fib1, (fib1 << 1) + fib0
 n += 2
 # print(f"{n=} {fib1=}")
 if gmpy2.is_prime(fib1):
 yield n, fib1
print("3,2")
print("4,3")
start_time=time.time()
for n, fib in fib_prime_gen():
 print(f"{n},{fib}")
 if n>=9677:
 break
print(f"F9677 in {time.time()-start_time} seconds.")