[Haskell-cafe] faster factorial function via FFI?

[Dan Weston]

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