読者です 読者をやめる 読者になる 読者になる

にわとり小屋でのプログラミング ブログ

名古屋のCoqエンジニアです。日記は勘で書いてるところがあるので細かいとこが違うことがあります。

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.

ソースコード

関連記事
http://projecteuler.net/index.php?section=problems&id=44
http://www.haskell.org/haskellwiki/Euler_problems/41_to_50#Problem_44
http://www.euphe.net/memo/euler03.html
http://d.hatena.ne.jp/tanakaBox/20080523/1211478949
http://d.hatena.ne.jp/rst76/20080421/1208787760