Problem 160 （階乗の非零末尾数字)

http://projecteuler.net/index.php?section=problems&id=160
階乗の非零末尾の数字をk個求める。
n!に含まれる5の数をp(n)とすると
n!/10^p(n) mod (10^k) が分かればよい。
まず、n!/10^p(n) = 0 (mod 2^k) であるから
n!/10^p(n) mod (5^k) を求めて、あとは連立合同式
x = 0 (mod 2^k), x = ? (mod 5^k)
を解けばよい。
n!/5^p(n) mod (5^k) から　n!/10^p(n) mod (5^k)　を求めている。

import Data.List
import Mod
-- p:odd prime, k: natural number
-- productMod (p^k) [ n | n <- [1..p^k-1], mod n p /= 0] = - 1 (mod p^k) (Generalization of Wilson's theorem)
-- 2^(4*5^(k-1)) = 1 (mod 5^k)  (Fermat's little theorem)

del5 n | mod n 5 == 0 = del5 $ div n 5
       | otherwise = n

-- solve n!/5^p(n) mod (5^k)
solve5 n k = f $ divMod n m
    where m = 5^k
          f (0,0) = 1
          f (q,r) | even q = proMod5 (q*m) (q*m+r) * g (q,r) `mod` m
                  | otherwise = - proMod5 (q*m) (q*m+r) * g (q,r) `mod` m
          g (q,r) = let (q',r') = divMod q 5
                    in f (q',r'*5^(k-1)) `mod` m
          proMod5 x y = productMod m.map del5$ [x+1..y]

-- n!/10^p(n) mod (5^k)
solve10 n k = let phi = 4*5^(k-1)
                  n5 = sum.takeWhile (>0).map (div n).iterate (*5)$ 5
                  r = mod n5 phi
                  a = powMod 2 (phi-r) (5^k)
                  b = solve5 n k
              in mod (a*b) $ 5^k

lastDigit n k = let a = invMod (2^k) (5^k)
                    b = solve10 n k
                in mod (a*b*2^k) (10^k)

main = print.lastDigit (10^12) $ 5

Modには剰余に関する関数が定義されている。

module Mod where
import Data.List
-- input; m,n:正整数
-- output; a,b s.t. am+bn=gcd(m,n)
gcd' m n | m < n = let (a,b) = gcd'  n m in (b,a)
         | r == 0 = (1,1-q)
         | otherwise = let (a,b) = gcd' n r
                       in (b,a-b*q)
         where (q,r) = divMod m n 

-- product xs (mod p)
productMod p = foldl' (mulMod p) 1
    where mulMod p a b = mod (a*b) p

-- xn=1(mod p) なる x を返す
invMod n = fst.gcd' n

-- a^n (mod p)
powMod a n p | n == 0 = 1
             | even n = powMod (mod (a*a) p) (div n 2) p
             | otherwise = (a * powMod a (n-1) p)`mod` p