-- Compute 9 rightmost decimal digits of n!, disregarding trailing zeroes module FactLimit where import Data.Bits -- The (executable) specification spec :: Int -> Int spec = fromIntegral . lastnz9 . fact . fromIntegral fact :: Integer -> Integer fact n = product [1..n] lastnz9 :: Integral a => a -> a lastnz9 n = drop0s n `mod` 10^9 drop0s n | (q,0) <- quotRem n 10 = drop0s q | otherwise = n -- The optimal version -- Accumulating factorial fact_accum' :: Int -> Int -> Integer -> Integer fact_accum' n i acc | i > n = acc | otherwise = fact_accum' n (i+1) (acc * fromIntegral i) fact_accum n = fact_accum' n 1 1 fact_accum_9' :: Int -> Int -> Int -> Int fact_accum_9' n i acc | i > n = acc | otherwise = fact_accum_9' n (i+1) (lastnz9 (acc * i)) fact_accum_9 n = fact_accum_9' n 1 1 -- Testing test1 = [spec 37 == fact_accum_9 37, spec 53 == fact_accum_9 53, spec 100 == fact_accum_9 100] -- All True testsmall = all (\n -> spec n == fact_accum_9 n) [0..21] testbad = all (\n -> spec n == fact_accum_9 n) [0..25] tlz90 = lastnz9 (25 * (4*10^9+16)) == lastnz9 (25 * lastnz9 (4*10^9+16)) tlz91 = lastnz9 (25 * (10^9+16)) == lastnz9 (25 * lastnz9 (10^9+16)) -- drop0s (fact 24) is 62044840173323943936 -- drop0s (fact 24) `div` 10^9 is 62044840173 -- and it is indeed not divisible by 4. -- downward counting fact_accum_9d :: Int -> Int -> Int -> Int fact_accum_9d n i acc | i < 2 = acc | otherwise = fact_accum_9d n (i-1) (drop0s (acc * i) `mod` 10^9) -- Correct answer -- return (n',k) such that n == n'*f^k factors n f = go 0 n where go k n | (q,0) <- quotRem n f = go (k+1) q | otherwise = (n,k) modulus = 10^9 -- We assume that Int is 64-bit long. It is left as an exercise -- to the reader how to modify the program for 32-bit Ints. factl :: Int -> Int factl n = go 2 0 1 where go i twos acc | i > n = shifts acc twos go i twos acc = let (i1,k1) = factors i 2 (i2,k2) = factors i1 5 in go (i+1) (twos + k1 - k2) (acc * i2 `mod` modulus) -- compute (acc * 2^twos) `mod` modulus shifts acc 0 = acc shifts acc twos | twos >= 4 = shifts ((acc `shift` 4) `mod` modulus) (twos - 4) | otherwise = (shift acc twos) `mod` modulus test_mod = all (\n -> spec n == factl n) [0..45]