## Fibonacci Primes

### October 29, 2010

This is the second level of the Greplin Programming Challenge; we saw the first level in a previous exercise. It’s easy enough to find the solution without any programming if you know your way around the internet: Neil J. A. Sloane gives the prime fibonacci numbers at A005478 and Dario Alpern’s factorization applet finds the factors. And it’s also easy if you’ve been following Programming Praxis for a while, because we have written functions to compute fibonacci numbers, test for primality, and factor large integers. So we’ll look at this exercise from the point of view of a college sophomore who has been given this question on the final exam of his first algorithms course.

We’ll begin at the end with the top-level function that computes the sequence of fibonacci numbers, each number the sum of the two prior numbers in the sequence, starting with 0 and 1. As we compute the fibonacci numbers, we test each one for primality and stop when we find the desired result:

```(define (greplin2 n)   (let loop ((f-2 0) (f-1 1) (f 1))     (if (and (< n f) (prime? f))         (sum (factors (+ f 1)))         (loop f-1 f (+ f-1 f)))))```

Here `f` is the current fibonacci number, `f-1` is the previous fibonacci number, and `f-2` is the fibonacci number before that. The function loops until it finds a prime fibonacci number larger than the input number, then computes the required factors, sums them, and reports the result.

Next we look at the integer factorization part of the problem. Most algorithms courses don’t go beyond Pollard’s rho algorithm, but we’ll use an even simpler algorithm, trial division, dividing by the integers in succession starting from 2 and collecting the factors in reverse order in the `fs` list; we could do better by testing 2 and then the odd integers starting from 3 and adding 2 at each step, but we’ll leave that improvement for another time. Note that `x` doesn’t advance when it is a factor, in case the same factor appears more than once:

```(define (factors n)   (let loop ((n n) (x 2) (fs '()))     (cond ((< n (* x x)) (reverse (cons n fs)))           ((zero? (modulo n x)) (loop (/ n x) x (cons x fs)))           (else (loop n (+ x 1) fs)))))```

Now we’ll turn our attention to the function that determines if a number is prime. Most algorithms texts explain the Miller-Rabin probabilistic method for checking primality. The method is based on Fermat’s little theorem, which states that if n is prime, then for any integer a, ana (mod n), which is equivalent to an−1 ≡ 1 (mod n) after dividing both sides by a. Thus, if we find an a for which Fermat’s little theorem fails, n is certainly composite; we say that a is a witness to the compositeness of n. Unfortunately, the opposite is not true; just because an a passes Fermat’s little theorem, it is not necessarily true that n is prime. In fact, there are some numbers, such as 561, for which there is no a that causes Fermat’s little theorem to fail; these numbers are called Carmichael numbers, because they were studied by Robert Carmichael, and though they are quite rare, they force us to use a somewhat more complicated method to check for primality.

Our procedure `witness?` works like this. Given an odd integer n, let n = 2r · s + 1 with s odd. Then choose a random integer a with 1 < a < n. If as ≡ 1 (mod n) or a2j·s ≡ −1 (mod n) for some 0 ≤ jr−1, then n passes the test. A prime n will pass the test for all a. Here we check if a is a witness to the compositeness of n:

```(define (witness? a n)   (let loop ((r 0) (s (- n 1)))     (if (even? s) (loop (+ r 1) (/ s 2))       (if (= (expm a s n) 1) #t         (let loop ((j 0) (s s))           (cond ((= j r) #f)                 ((= (expm a s n) (- n 1)) #t)                 (else (loop (+ j 1) (* s 2)))))))))```

`Loop1` computes r and s. When s becomes odd, `witness?` checks if as ≡ 1 (mod n), and if so, returns `#t` indicating that a is a witness to the compositeness of n. Otherwise, `witness?` enters `loop2`, which loops for j from 0 to r−1, doubling s at each step to make the calculation 2j·s and returning `#t` if any as ≡ −1 (mod n) or `#f` when j reaches r. Note that `#t` means that a is a witness to the compositeness of n, but `#f` only means that a is not a witness to the compositeness of n, not that n is prime.

The `witness?` function calls a function `expm` that calculates be (mod n). One way to compute that is by repeated multiplication, applying the modulo operator at each step so the intermediate results don’t get too big. A better approach for large e is the “square and multiply” algorithm that loops over the binary encoding of the exponent, squaring when the bit is zero (that is, when the exponent e is even) and doubling when the bit is one (that is, when the exponent e is odd):

```(define (expm b e m)   (define (m* x y) (modulo (* x y) m))   (cond ((zero? e) 1)         ((even? e) (expm (m* b b) (/ e 2) m))         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))```

Let’s return to the primality checker. By testing a large number of possible witnesses, chosen at random, we can determine to any desired degree of probability that a number is prime; if any of k numbers are witnesses to the compositeness of n, then n is certainly not prime, but if none of k numbers are witnesses to the compositeness of n, then n is prime with probability 4k; choosing k=50 is sufficient for all but the largest numbers. Our function first eliminates even numbers which can never be prime:

```(define (prime? n)   (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f)         (else (let loop ((k 50))                 (cond ((zero? k) #t)                       ((not (witness? (randint 1 n) n)) #f)                       (else (loop (- k 1))))))))```

The `prime?` function requires random integers a for the `witness?` test. Most programming languages have a function that returns random integers, but since our hypothetical algorithms student probably studied how to generate a sequence of random numbers, we’ll write that code ourselves. We use a linear congruential random number generator Xk+1 = Xk · a (mod m); Park and Miller suggested the random number generator with a = 75 = 16807 and m = 231 − 1 = 2147483647. Calling `(rand)` gets the next number in the sequence; calling `(rand `seed`)`, with 0 < seed < m, resets the random number seed before computing the next number in the sequence:

```(define rand   (let ((a 16807) (m 2147483647) (seed 1043618065))     (lambda args       (if (pair? args) (set! seed (car args)))       (set! seed (modulo (* a seed) m))       (/ seed m))))```

`Rand` returns a rational number strictly between 0 and 1. Given `rand`, it is easy to write a function that returns random integers:

```(define (randint first past)   (+ first (inexact->exact (floor (* (rand) (- past first))))))```

The only remaining piece is `sum`, which is simple:

`(define (sum xs) (apply + xs))`

And now we are ready to solve the problem:

```> (greplin2 227000) 352```

The smallest prime fibonacci number greater than 227000 is 514229. The factors of 514229 + 1 are 2, 3, 5, 61 and 281, and the sum of those factors is 352.

You can run the program at http://programmingpraxis.codepad.org/8i16Cbx7.

Pages: 1 2

### 8 Responses to “Fibonacci Primes”

1. […] Praxis – Fibonacci Primes By Remco Niemeijer Today’s Programming Praxis exercise comes from the same challenge as the longest palindrome exercise from […]

```fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

factors :: Integer -> [Integer]
factors = f 2 where
f d n | n < d        = []
| mod n d == 0 = d : f (d + 1) (div n d)
| otherwise    = f (d + 1) n

isPrime :: Integer -> Bool
isPrime n = factors n == [n]

greplin :: Integer -> Integer
greplin n = head [sum \$ factors (f + 1) | f <- fibs, f > n, isPrime f]
```
3. Graham said

I admit, I used Sage’s (see sagemath.org) sweet built-in routines when working on the Greplin Challenge.
It’s certainly cheating. Congrats to Remco on the impressive Haskell, and you for your Scheme prowess and
excellent blog.

```sage: def greplin(n):
....:     for f in fibonacci_xrange(n, 10*n):
....:         if is_prime(f):
....:             return sum([b for b, e in factor(1 + f)])
....:
sage: greplin(227000)
352
```
4. fengshaun said

Damn, my code is so similar to Remco’s:

``` import qualified Data.Numbers.Primes as P```

``` fibs :: [Integer] fibs = 0:1:zipWith (+) fibs (tail fibs) smallestFibPrimeAfter :: Integer -> Integer smallestFibPrimeAfter lowBound = head . dropWhile (\a -> not . P.isPrime \$ a) . filter (>lowBound) \$ fibs factors :: Integer -> [Integer] factors x = 1:factors' x 2 where factors' n f | f >= (n `div` 2) = [n] | n `mod` f == 0 = f:factors' (n `div` f) (f+1) | otherwise = factors' n (f+1) ```

```main :: IO () main = print . show . sum . factors . (+ 1) . smallestFibPrimeAfter \$ 227000 ```

5. Axio said

;; The question’s specification is incomplete.
;; You want the (unique) *prime* factors, not the factors (divisors)
;; I read the specification as (factors 12) => (1 2 3 4 6 12)
;; while you expected (factors 12) => (2 3)
;; With regard to this issue,the example given is of no help at all
;; since 14 is the product of two primes.

;; Tail recursive
(define (fib n #!optional (g 0) (d 1))
(if (= 0 n)
g
(fib (- n 1) d (+ g d))))

;; Excluding 1 and oneself
;; Quite naive, but not completely
(define (factors n)
(let ((tmax (floor (sqrt n))))
(let loop ((i 2))
(if (> i tmax)
‘()
(if (zero? (modulo n i))
(cons i (cons (quotient n i) (loop (+ i 1))))
(loop (+ i 1)))))))

;; Unefficient, since we’ve already done most computations
(define (prime? n)
(and (not (= 1 n)) (null? (factors n))))

(define (run n)
(let loop ((i 2))
(let ((f (fib i)))
(if (>= f n)
(if (null? (factors f)) ;; prime
(apply + (filter prime? (factors (+ 1 f))) )
(loop (+ i 1)))
(loop (+ i 1))))))

(run 227000) ;; => 352

6. Khanh Nguyen said
```open System

let divisors (n:int) = [
let x = float n / 2.0 |> int
for i in 2..x do
if n % i = 0 then yield i
]

let isPrime (n:int) = List.length (divisors n)  = 0

let fib_primes (n:int) =
let smallest_fib =
let rec fin_fibs a b =
if isPrime a && a > n then a
else
fin_fibs b (a+b)
fin_fibs 1 1

(smallest_fib + 1) |> divisors |> List.filter (fun x -> isPrime x) |> List.sum

fib_primes 227000
```
7. danaj said

A version in Perl that works for bigints. E.g.

```perl greplin2.pl 227000
Using F(29) = 514229
Sum of factors: 352

perl greplin2.pl 2000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
Using F(449) = 3061719992484545030554313848083717208111285432353738497131674799321571238149015933442805665949
Sum of factors: 4966336310425739390388368472890698620
```

I reused some bits I wrote for a parallel Fibonacci prime finder. The Fn optimization is interesting but not that useful here. The fib generator could simpler for this code. You need the Math::Prime::Util::GMP module installed to make bigint factoring work at reasonable speeds (or use Math::Pari). Math::Primality and Math::Pari could be used for primality testing instead.

```#!/usr/bin/env perl
use strict;
use warnings;
use Math::BigInt;
use Math::Prime::Util ':all';
use List::Util qw/sum/;

my \$n = \$ARGV[0] || 227000;
my @fibstate;
my \$F;
my \$fn = 0;
do {
# All prime F(n) have n prime when n > 4.  Useful for very large n.
\$fn = (\$fn <= 4) ? \$fn+1 : next_prime(\$fn);
\$F = fib_n(\$fn, \@fibstate);
} while (\$F <= \$n || !is_prime(\$F));
print "Using F(\$fn) = \$F\n";
print "Sum of factors: ", sum(factor(\$F+1)), "\n";

sub fib_n {
my (\$n, \$fibstate) = @_;
@\$fibstate = (1, Math::BigInt->new(0), Math::BigInt->new(1))
unless defined \$fibstate->[0];
my (\$curn, \$a, \$b) = @\$fibstate;
die "fib_n only increases" if \$n < \$curn;
do { (\$a, \$b) = (\$b, \$a+\$b); } for (1 .. \$n-\$curn);
@\$fibstate = (\$n, \$a, \$b);
\$b;
}
```