euler002.mlw 2.7 KB
Newer Older
1
(* Euler Project, problem 2
2 3 4 5 6 7 8 9

Each new term in the Fibonacci sequence is generated by adding the
previous two terms. By starting with 1 and 2, the first 10 terms will
be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

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 13 14
(** {2 Sum of even-valued Fibonacci numbers} *)

theory FibSumEven
MARCHE Claude's avatar
MARCHE Claude committed
15

16 17 18
  use int.Int
  use int.Fibonacci
  use int.ComputerDivision
MARCHE Claude's avatar
MARCHE Claude committed
19 20 21 22 23 24 25 26

  (* [fib_sum_even m n] is the sum of even-valued terms of the
      Fibonacci sequence from index 0 to n-1, that do not exceed m *)
  function fib_sum_even int int : int

  axiom SumZero: forall m:int. fib_sum_even m 0 = 0

  axiom SumEvenLe: forall n m:int.
27
     n >= 0 -> fib n <= m -> mod (fib n) 2 = 0 ->
MARCHE Claude's avatar
MARCHE Claude committed
28 29 30
       fib_sum_even m (n+1) = fib_sum_even m n + fib n

  axiom SumEvenGt: forall n m:int.
31
     n >= 0 -> fib n > m -> mod (fib n) 2 = 0 ->
MARCHE Claude's avatar
MARCHE Claude committed
32 33 34
       fib_sum_even m (n+1) = fib_sum_even m n

  axiom SumOdd: forall n m:int.
35
     n >= 0 -> mod (fib n) 2 <> 0 ->
MARCHE Claude's avatar
MARCHE Claude committed
36 37
       fib_sum_even m (n+1) = fib_sum_even m n

38 39
  predicate is_fib_sum_even (m:int) (sum:int) =
    exists n:int.
MARCHE Claude's avatar
MARCHE Claude committed
40 41 42
      sum = fib_sum_even m n /\ fib n > m
   (* Note: we take for granted that [fib] is an
        increasing sequence *)
43

44 45
end

46
module FibOnlyEven
47

48 49 50
  use int.Int
  use int.ComputerDivision
  use int.Fibonacci
51

52 53 54 55 56
  let rec lemma fib_even_3n (n:int)
    requires { n >= 0 }
    variant { n }
    ensures { mod (fib n) 2 = 0 <-> mod n 3 = 0 }
  = if n > 2 then fib_even_3n (n-3)
57

58
  function fib_even (n: int) : int = fib (3 * n)
59

60 61
  lemma fib_even0: fib_even 0 = 0
  lemma fib_even1: fib_even 1 = 2
62

63
  lemma fib_evenn: forall n:int [fib_even n].
64 65 66 67 68 69
     n >= 2 -> fib_even n = 4 * fib_even (n-1) + fib_even (n-2)

end

module Solve

70 71 72 73 74
  use int.Int
  use ref.Ref
  use int.Fibonacci
  use FibSumEven
  use FibOnlyEven
75

76 77
  let f m : int
    requires { m >= 0 }
78 79 80
    ensures  { exists n:int. result = fib_sum_even m n /\ fib n > m }
  = let x = ref 0 in
    let y = ref 2 in
81
    let sum = ref 0 in
MARCHE Claude's avatar
MARCHE Claude committed
82
    let ghost n = ref 0 in
83
    let ghost k = ref 0 in
MARCHE Claude's avatar
MARCHE Claude committed
84 85
    while !x <= m do
      invariant { !n >= 0 }
86
      invariant { !k >= 0 }
MARCHE Claude's avatar
MARCHE Claude committed
87 88 89 90 91
      invariant { !x = fib_even !n }
      invariant { !x = fib !k }
      invariant { !y = fib_even (!n+1) }
      invariant { !y = fib (!k+3) }
      invariant { !sum = fib_sum_even m !k }
92
      invariant { 0 <= !x < !y }
93
      variant { m - !x }
94 95 96
      let tmp = !x in
      x := !y;
      y := 4 * !y + tmp;
MARCHE Claude's avatar
MARCHE Claude committed
97 98 99
      sum := !sum + tmp;
      n := !n + 1;
      k := !k + 3
100
    done;
101
    !sum
102

103
  let run () = f 4_000_000 (* should be 4613732 *)
104

MARCHE Claude's avatar
MARCHE Claude committed
105 106 107
  exception BenchFailure

  let bench () raises { BenchFailure -> true } =
108 109 110
    let x = run () in
    if x <> 4613732 then raise BenchFailure;
    x
MARCHE Claude's avatar
MARCHE Claude committed
111

112
end