6
\$\begingroup\$

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.

\$\endgroup\$
4
  • \$\begingroup\$ 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). \$\endgroup\$
    – Arnauld
    Jun 18 at 15:02
  • \$\begingroup\$ (But any primality test is likely to reject the other ones immediately anyway.) \$\endgroup\$
    – Arnauld
    Jun 18 at 15:09
  • 1
    \$\begingroup\$ @Arnauld Those are already rejected by testing the index (\$n\$) for primality. \$\endgroup\$ Jun 18 at 17:40
  • \$\begingroup\$ @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. \$\endgroup\$
    – Aiden4
    Jun 19 at 4:36

3 Answers 3

7
\$\begingroup\$

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}");
    }
}
\$\endgroup\$
2
  • \$\begingroup\$ 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? \$\endgroup\$
    – Aiden4
    Jun 22 at 6:36
  • \$\begingroup\$ @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. \$\endgroup\$ Jun 22 at 16:35
3
\$\begingroup\$

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.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – Aiden4
    Jun 22 at 6:31
1
\$\begingroup\$

Wolfram Language (Mathematica), F9677 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?

\$\endgroup\$
4
  • \$\begingroup\$ The challenge is to calculate every Fibonacci prime, not just to test one of them. (Also, there are reasons we don’t score fastest-code questions using TIO.) \$\endgroup\$ Jun 18 at 23:13
  • \$\begingroup\$ Duly noted. But I did test more than one Fibonacci number for prime-ness. Select[Table[(MatrixPower[{{1, 1}, {1, 0}}, n] . {1, 0})[[2]], {n, 1, 9677}], PrimeQ] will generate a list of the Fibonacci numbers that are prime, up to F_9677. \$\endgroup\$
    – dirvine
    Jun 19 at 0:20
  • 2
    \$\begingroup\$ That’s neither the code you listed nor the code you benchmarked. The tiebreaker score is supposed to be the total time from the start of the sequence, not an individual time for one number. Also: there doesn’t seem to be a known error bound for PrimeQ as required by the challenge; and the specific binary output format, although I think it’s silly, is part of the challenge too unless the OP decides to amend it. \$\endgroup\$ Jun 19 at 0:40
  • 2
    \$\begingroup\$ For the reasons listed above, this answer is invalid. I need a single runnable program that produces output complying with the now text-based format, and a command-line invocation that can run it. Also, see this meta question for why I'm doing this on my own computer rather than tio. \$\endgroup\$
    – Aiden4
    Jun 19 at 4:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.