examples/euler002.mlw

   1 (* Euler Project, problem 2
   2
   3 Each new term in the Fibonacci sequence is generated by adding the
   4 previous two terms. By starting with 1 and 2, the first 10 terms will
   5 be:
   6
   7 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
   8
   9 By considering the terms in the Fibonacci sequence whose values do not
  10 exceed four million, find the sum of the even-valued terms. *)
  11
  12 (** {2 Sum of even-valued Fibonacci numbers} *)
  13
  14 theory FibSumEven
  15
  16   use int.Int
  17   use int.Fibonacci
  18   use int.ComputerDivision
  19
  20   (* [fib_sum_even m n] is the sum of even-valued terms of the
  21       Fibonacci sequence from index 0 to n-1, that do not exceed m *)
  22   function fib_sum_even int int : int
  23
  24   axiom SumZero: forall m:int. fib_sum_even m 0 = 0
  25
  26   axiom SumEvenLe: forall n m:int.
  27      n >= 0 -> fib n <= m -> mod (fib n) 2 = 0 ->
  28        fib_sum_even m (n+1) = fib_sum_even m n + fib n
  29
  30   axiom SumEvenGt: forall n m:int.
  31      n >= 0 -> fib n > m -> mod (fib n) 2 = 0 ->
  32        fib_sum_even m (n+1) = fib_sum_even m n
  33
  34   axiom SumOdd: forall n m:int.
  35      n >= 0 -> mod (fib n) 2 <> 0 ->
  36        fib_sum_even m (n+1) = fib_sum_even m n
  37
  38   predicate is_fib_sum_even (m:int) (sum:int) =
  39     exists n:int.
  40       sum = fib_sum_even m n /\ fib n > m
  41    (* Note: we take for granted that [fib] is an
  42         increasing sequence *)
  43
  44 end
  45
  46 module FibOnlyEven
  47
  48   use int.Int
  49   use int.ComputerDivision
  50   use int.Fibonacci
  51
  52   let rec lemma fib_even_3n (n:int)
  53     requires { n >= 0 }
  54     variant { n }
  55     ensures { mod (fib n) 2 = 0 <-> mod n 3 = 0 }
  56   = if n > 2 then fib_even_3n (n-3)
  57
  58   function fib_even (n: int) : int = fib (3 * n)
  59
  60   lemma fib_even0: fib_even 0 = 0
  61   lemma fib_even1: fib_even 1 = 2
  62
  63   lemma fib_evenn: forall n:int [fib_even n].
  64      n >= 2 -> fib_even n = 4 * fib_even (n-1) + fib_even (n-2)
  65
  66 end
  67
  68 module Solve
  69
  70   use int.Int
  71   use ref.Ref
  72   use int.Fibonacci
  73   use FibSumEven
  74   use FibOnlyEven
  75
  76   let f m : int
  77     requires { m >= 0 }
  78     ensures  { exists n:int. result = fib_sum_even m n /\ fib n > m }
  79   = let x = ref 0 in
  80     let y = ref 2 in
  81     let sum = ref 0 in
  82     let ghost n = ref 0 in
  83     let ghost k = ref 0 in
  84     while !x <= m do
  85       invariant { !n >= 0 }
  86       invariant { !k >= 0 }
  87       invariant { !x = fib_even !n }
  88       invariant { !x = fib !k }
  89       invariant { !y = fib_even (!n+1) }
  90       invariant { !y = fib (!k+3) }
  91       invariant { !sum = fib_sum_even m !k }
  92       invariant { 0 <= !x < !y }
  93       variant { m - !x }
  94       let tmp = !x in
  95       x := !y;
  96       y := 4 * !y + tmp;
  97       sum := !sum + tmp;
  98       n := !n + 1;
  99       k := !k + 3
 100     done;
 101     !sum
 102
 103   let run () = f 4_000_000 (* should be 4613732 *)
 104
 105   exception BenchFailure
 106
 107   let bench () raises { BenchFailure -> true } =
 108     let x = run () in
 109     if x <> 4613732 then raise BenchFailure;
 110     x
 111
 112 end