# [Haskell-cafe] faster factorial function via FFI?

Dan Weston westondan at imageworks.com
Mon Apr 23 23:02:24 EDT 2007

```Why all the fuss? n! is in fact very easily *completely* factored into
prime numbers (and this never changes, so you can just read in a binary
lookup table). It took my machine only 3 seconds to find the prime
factors of the first 1000 factorials:

Then you can do for example:

import List

binomCoeffFactors :: Int -> Int -> [Int]
binomCoeffFactors n k = let pfs = primeFactorials in
(pfs !! n \\ pfs !! k)
\\ pfs !! (n-k)

binomCoeff :: Int -> Int -> Int
binomCoeff n k = foldr (*) 1 (binomCoeffFactors n k)

-- List of primes. This is reasonably fast but doesn't have to be,
-- since we only need it this once to generate the lookup table
primes = (2 : [x | x <- [3,5..], isPrime x])

-- Efficient test presupposes the existence of primes
-- This works because to determine whether p is prime you only need
-- to know the primes strictly less than p (except for 2 of course!)
isPrime x = null divisors
where divisors      = [y | y <- onlyUpToSqrtX primes, x `mod` y == 0]
onlyUpToSqrtX = fst . span (<= sqrtX)
sqrtX         = floor (sqrt (fromIntegral x))

primeFactor' :: [Int] -> Int -> [Int]
primeFactor' p@(h:t) n
| n < 2 = []
| otherwise = let (d,m) = n `divMod` h in
if m == 0 then h : primeFactor' p d
else     primeFactor' t n

primeFactors :: Int -> [Int]
primeFactors = let p = primes in primeFactor' p

primeFactorial' :: (Int, [Int]) -> (Int, [Int])
primeFactorial' (n,pfs) = (n1,ret)
where pfn = primeFactors n1
ret = merge pfn pfs
n1  = n + 1

primeFactorials :: [[Int]]
primeFactorials = map snd . iterate primeFactorial' \$ (0,[])

-- timingTest 10000 takes only 3 seconds!
timingTest :: Int -> Int
timingTest nMax = let pfs = primeFactorials in
sum . map sum . take nMax \$ pfs

-- PRE : Args  are sorted
-- POST: Result is sorted union
merge as []  = as
merge []  bs = bs
merge aas@(a:as) bbs@(b:bs) | a < b     = a : merge as  bbs
| otherwise = b : merge aas bs

```