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