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 ==  = 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 (,) = 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 (>).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 ==  = (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 ==  = 1
| even n = powMod (mod (a*a) p) (div n 2) p
| otherwise = (a * powMod a (n-1) p)`mod` p