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

ここからがいよいよ証明

Require Import ZArith.
Open Scope Z_scope.

Definition pent2 (k:Z) := k * (3*k - 1).
Definition D (i n:Z) := pent2 (n+i) - pent2 n.

Theorem D_increase_n : forall i m n:Z, 1 <= i -> 1 <= m -> 1 <= n -> m <= n -> D i m <= D i n.
Proof.
intros i m n H1 H2 H3 H4.
apply Zle_0_minus_le.
cut (D i n - D i m = 6 * i * (n-m)).
(* sub1 *)
intro eq; rewrite eq.
apply Zmult_le_0_compat; omega.
(* sub2 *)
unfold D; unfold pent2.
ring.
Qed.

Theorem D_increase_i : forall i j n:Z, i >= 1 -> j >= 1 -> n >= 1 -> i <= j -> D i n <= D j n.
Proof.
intros i j n H1 H2 H3 H4.
unfold D.
apply (Zplus_le_compat_r (pent2(n+i)) (pent2(n+j)) (- pent2 n)).
unfold pent2.
 apply Zmult_le_compat; omega.
Qed.

ソースコード