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.

About these ads

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 [...]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2010/10/29/programming-praxis-fibonacci-primes/ for a version with comments):

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 630 other followers

%d bloggers like this: