Programming Praxis

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, aⁿ ≡ a (mod n), which is equivalent to aⁿ⁻¹ ≡ 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 = 2^r · s + 1 with s odd. Then choose a random integer a with 1 < a < n. If a^s ≡ 1 (mod n) or a^2j·s ≡ −1 (mod n) for some 0 ≤ j ≤ r−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 a^s ≡ 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 2^j·s and returning #t if any a^s ≡ −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 b^e (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 4^−k; 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 X_k+1 = X_k · a (mod m); Park and Miller suggested the random number generator with a = 7⁵ = 16807 and m = 2³¹ − 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