CoqでProject Euler (problem 44)

Problem44が解けた? - みずぴー日記に触発されて、Project EulerのProblem 44に挑戦してみた。

問題44:
五角数は Pk = k(3k-1)/2で生成される. 
五角数のペア PmとPnについて, 差と和が五角数になるものを考える. 
差を D = |Pm - Pn| と書く. 
差 D の最小値を求めよ.

この問題は D の最小性がポイントになる。見付けたほとんどのプログラムは、扱う五角数の上限を固定して計算していたが、それでは厳密には D 自体の最小性を保証することはできない。
今回は D の最小性を保証するプログラムと、その根拠となる単調性をCoqで証明した。

アイデアは次のとおり:
m = n+i とおき、五角数の差D を考えると、Dはi, nそれぞれにおいて単調増加になる。
この単調増加性を使えば、ある基準値Mに対して、D(i,n)がMを越えないようなすべての五角数の組を調べることができる。

厳密さを保つために任意精度演算ライブラリNumsのbigなintを使用したが、わりと現実的な時間内に応答を得ることができた。

実行結果:

Answer: 5482660

real    1m52.673s
user    1m52.563s
sys     0m0.059s

(**
   solve.ml [projecto euler problem 44]
   2008/06/23 Y. Imai wrote.
   compile $ ocamlopt nums.cmxa solve.ml
**)

open List
open Big_int
let ( @@ ) f x = f x

let big = big_int_of_int
let ( =. ), ( >. ), ( *. ), ( +. ), ( -. ), ( /. ) =
  eq_big_int,
  gt_big_int,
  mult_big_int,
  add_big_int,
  sub_big_int,
  div_big_int
let pent n = (big 3 *. n *. n -. n) /. big 2
let is_pent y =
  let s = big 1 +. big 24 *. y in
    square_big_int (sqrt_big_int s) =. s &&
      mod_big_int (sqrt_big_int s +. big 1) (big 6) =. zero_big_int
      
let solve m =
  let store = ref [] in
  let rec aux i n =
    let diff = pent (big (n+i)) -. pent (big n) in
    let sum  = pent (big (n+i)) +. pent (big n) in
    match diff >. m, n = 1 with
      | true, true -> !store
      | true, false -> aux (i+1) 1
      | false, _ ->
	  if is_pent sum && is_pent diff then
	    store := diff :: !store;
	  aux i (n+1)
  in aux 1 1
    
let main =
  print_endline @@ "Answer: " ^ string_of_big_int @@ hd
  @@ sort compare_big_int @@ solve @@ big 5500000

ここからがいよいよ証明