Project Euler 44

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

真面目に因数分解すると速いです。

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 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 div_f f g =
   if f = [] then
      []
   else if g = [] then
      f
   else
      let f0 = List.head f
      let g0 = List.head g
      if fst f0 = fst g0 then
         if snd f0 = snd g0 then
            div_f (List.tail f) (List.tail g)
         else
            (fst f0, (snd f0) - (snd g0)) :: (div_f (List.tail f) (List.tail g))
      else
         div_f (List.tail f) g

let count n = Seq.initInfinite ((+) n)

let P' n = n * (3 * n - 1)
let P n = (P' n) / 2
let is_penta n =
   let m = 1 + 24 * n
   let l = int (sqrt (float m))
   l * l = m && (1 + l) % 6 = 0

// 差が小さい方から順に五角数のペアを出す
let seq_penta_pairs = seq {
   for l in count 1 do
      let v = P' l
      let f = (factorize l) @ (factorize (3 * l - 1))
      for d1 in List.map value (divisors f) do
         let d2 = v / d1
         if d2 % 3 = 2 then
            let d3 = (d2 + 1) / 3
            if d3 > d1 && (d1 + d3) % 2 = 0 then
               let k = (d1 + d3) / 2
               let j = d3 - k
               yield (P k, P j)
}

printfn "%d" (Seq.head (Seq.map (fun x -> (fst x) - (snd x))
                  (Seq.filter (fun x -> is_penta ((fst x) + (snd x)))
                                                   seq_penta_pairs)))