Fibonacci Primes
October 29, 2010
A recent programming challenge asked you to solve the following problem:
Find the smallest prime fibonacci number greater than a given input number. Add one to that fibonacci prime, find the factors of the result, and return the sum of the factors. For instance, if the input number is 10, the smallest fibonacci prime greater than 10 is 13, the factors of 13 + 1 = 14 are 2 and 7, and their sum is 9.
Your task is to write a function to solve that problem, then apply your function to the input number 227000. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] Praxis – Fibonacci Primes By Remco Niemeijer Today’s Programming Praxis exercise comes from the same challenge as the longest palindrome exercise from […]
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]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.
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
;; 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
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 227000ruby solution (http://codepad.org/zOatd6NO)
A version in Perl that works for bigints. E.g.
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; }