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
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
