Project Euler 26

http://projecteuler.net/index.php?section=problems&id=26

この問題はオイラーの定理を使うと計算が速くなるのですが、コードがかなり長くなってしまいますね。本当は大きいほうから調べるとさらに速いです。

let rec calc_exp n p =
   if n % p = 0 then
      let t = calc_exp (n / p) p
      (fst t + 1, snd t)
   else
      (0, n)

let rec fac n p =
   if n = 1 then
      []
   else if p * p > n then
      [(n, 1)]
   else
      let t = calc_exp n p
      if fst t > 0 then
         (p, fst t) :: (fac (snd t) (p + 1))
      else
         fac n (p + 1)

let rec pow n e = 
    if e = 0 then 1
    else
        let m = pow n (e / 2)
        if e % 2 = 1 then m * m * n else m * m

let factorize n = fac n 2
let value f = List.fold (fun x y -> x * (pow (fst y) (snd y))) 1 f

let phi n = 
    let f = factorize n
    List.fold (fun x p -> x * (p - 1) / p) n [ for g in f -> fst g ]

let rec divisors = function
      | [] -> [[]]
      | (p,e)::fs -> 
            [ for t in (divisors fs) do
                  for e' in [0..e]
                     -> if e' = 0 then t else (p, e')::t ]
let divs n = divisors (factorize n)

let rec div25 n b =
    if b then
        if n % 2 = 0 then div25 (n / 2) b else div25 n false
    else
        if n % 5 = 0 then div25 (n / 5) b else n

let rec pow_mod n e d =
    if e = 0 then 1
    else
        let m = pow_mod n (e / 2) d
        if e % 2 = 1 then (m * m * n) % d else m * m % d

let cycle n =
    let m = div25 n true
    let s = List.toSeq (List.sort (List.map value (divs (phi m))))
    if m = 1 then
        1
    else
        Seq.head (Seq.filter (fun e -> pow_mod 10 e m = 1) s)

let longer x y = if snd x >= snd y then x else y
let N = 1000
printfn "%d" (fst (List.fold longer (0, 0)
                    (List.zip [1..N-1] (List.map cycle [1..N-1]))))