euler001.mlw 4.58 KB
Newer Older
1
(** {1 Euler Project, problem 1}
2 3 4 5

If we list all the natural numbers below 10 that are multiples of 3 or
5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

6
Find the sum of all the multiples of 3 or 5 below 1000.*)
7

MARCHE Claude's avatar
MARCHE Claude committed
8 9
theory DivModHints

10 11
  use int.Int
  use int.ComputerDivision
12 13

  lemma mod_div_unique :
MARCHE Claude's avatar
MARCHE Claude committed
14
    forall x y q r:int. x >= 0 /\ y > 0 /\ x = q*y + r /\ 0 <= r < y ->
MARCHE Claude's avatar
MARCHE Claude committed
15
      q = div x y /\ r = mod x y
MARCHE Claude's avatar
MARCHE Claude committed
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

  lemma mod_succ_1 :
    forall x y:int. x >= 0 /\ y > 0 ->
      mod (x+1) y <> 0 -> mod (x+1) y = (mod x y) + 1

  lemma mod_succ_2 :
    forall x y:int. x >= 0 /\ y > 0 ->
      mod (x+1) y = 0 -> mod x y = y-1

  lemma div_succ_1 :
    forall x y:int. x >= 0 /\ y > 0 ->
      mod (x+1) y = 0 -> div (x+1) y = (div x y) + 1

  lemma div_succ_2 :
    forall x y:int. x >= 0 /\ y > 0 ->
      mod (x+1) y <> 0 -> div (x+1) y = (div x y)

MARCHE Claude's avatar
MARCHE Claude committed
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
  lemma mod2_mul2:
    forall x:int. mod (2 * x) 2 = 0

  lemma mod2_mul2_aux:
    forall x y:int. mod (y * (2 * x)) 2 = 0

  lemma mod2_mul2_aux2:
    forall x y z t:int. mod (y * (2 * x) + z * (2 * t)) 2 = 0

  lemma div2_mul2:
    forall x:int. div (2 * x) 2 = x

  lemma div2_mul2_aux:
    forall x y:int. div (y * (2 * x)) 2 = y * x

  lemma div2_add:
    forall x y:int. mod x 2 = 0 /\ mod y 2 = 0 ->
      div (x+y) 2 = div x 2 + div y 2

  lemma div2_sub:
    forall x y:int. mod x 2 = 0 /\ mod y 2 = 0 ->
      div (x-y) 2 = div x 2 - div y 2

end

theory TriangularNumbers

60 61 62
  use int.Int
  use int.ComputerDivision
  use int.Div2
63
  use DivModHints as DMH
MARCHE Claude's avatar
MARCHE Claude committed
64 65 66 67 68 69 70 71 72 73 74 75

  lemma tr_mod_2:
    forall n:int. n >= 0 -> mod (n*(n+1)) 2 = 0

  function tr (n:int) : int =  div (n*(n+1)) 2

  lemma tr_repr:
    forall n:int. n >= 0 -> n*(n+1) = 2 * tr n

  lemma tr_succ:
    forall n:int. n >= 0 -> tr (n+1) = tr n + n + 1

MARCHE Claude's avatar
MARCHE Claude committed
76 77
end

78 79 80

theory SumMultiple

81 82
  use int.Int
  use int.ComputerDivision
83

84
  (* [sum_multiple_3_5_lt n] is the sum of all the multiples
85
     of 3 or 5 below n] *)
Andrei Paskevich's avatar
Andrei Paskevich committed
86
  function sum_multiple_3_5_lt int : int
87 88 89 90

  axiom SumEmpty: sum_multiple_3_5_lt 0 = 0

  axiom SumNo : forall n:int. n >= 0 ->
Andrei Paskevich's avatar
Andrei Paskevich committed
91
    mod n 3 <> 0 /\ mod n 5 <> 0 ->
92 93 94
    sum_multiple_3_5_lt (n+1) = sum_multiple_3_5_lt n

  axiom SumYes: forall n:int. n >= 0 ->
Andrei Paskevich's avatar
Andrei Paskevich committed
95
    mod n 3 = 0 \/ mod n 5 = 0 ->
96 97
    sum_multiple_3_5_lt (n+1) = sum_multiple_3_5_lt n + n

98
  use TriangularNumbers
MARCHE Claude's avatar
MARCHE Claude committed
99 100

  function closed_formula_aux (n:int) : int =
101 102 103
    let n3 = div n 3 in
    let n5 = div n 5 in
    let n15 = div n 15 in
MARCHE Claude's avatar
MARCHE Claude committed
104
    3 * tr n3 + 5 * tr n5 - 15 * tr n15
105

MARCHE Claude's avatar
MARCHE Claude committed
106
  predicate p (n:int) = sum_multiple_3_5_lt (n+1) = closed_formula_aux n
107

108
  use DivModHints as DMH
MARCHE Claude's avatar
MARCHE Claude committed
109

110 111 112 113
  lemma mod_15:
    forall n:int.
      mod n 15 = 0 <-> mod n 3 = 0 /\ mod n 5 = 0

MARCHE Claude's avatar
MARCHE Claude committed
114
  lemma Closed_formula_0: p 0
MARCHE Claude's avatar
MARCHE Claude committed
115

MARCHE Claude's avatar
MARCHE Claude committed
116
  lemma Closed_formula_n:
MARCHE Claude's avatar
MARCHE Claude committed
117 118
    forall n:int. n > 0 -> p (n-1) ->
      mod n 3 <> 0 /\ mod n 5 <> 0 -> p n
119

MARCHE Claude's avatar
MARCHE Claude committed
120
  lemma Closed_formula_n_3:
MARCHE Claude's avatar
MARCHE Claude committed
121 122
    forall n:int. n > 0 -> p (n-1) ->
      mod n 3 = 0 /\ mod n 5 <> 0 -> p n
MARCHE Claude's avatar
MARCHE Claude committed
123 124

  lemma Closed_formula_n_5:
MARCHE Claude's avatar
MARCHE Claude committed
125 126
    forall n:int. n > 0 -> p (n-1) ->
      mod n 3 <> 0 /\ mod n 5 = 0 -> p n
MARCHE Claude's avatar
MARCHE Claude committed
127 128

  lemma Closed_formula_n_15:
MARCHE Claude's avatar
MARCHE Claude committed
129 130
    forall n:int. n > 0 -> p (n-1) ->
      mod n 3 = 0 /\ mod n 5 = 0 -> p n
131

MARCHE Claude's avatar
MARCHE Claude committed
132
  constant b : int = 0
133

MARCHE Claude's avatar
MARCHE Claude committed
134 135 136
  clone int.Induction as I with constant bound = b, predicate p = p

  lemma Closed_formula_ind:
137 138
    forall n:int. 0 <= n -> p n

MARCHE Claude's avatar
MARCHE Claude committed
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
  function closed_formula (n:int) : int =
    let n3 = div n 3 in
    let n5 = div n 5 in
    let n15 = div n 15 in
    div (3 * (n3 * (n3+1)) +
         5 * (n5 * (n5+1)) -
         15 * (n15 * (n15+1))) 2

  lemma div_15: forall n:int. 0 <= n -> div n 15 >= 0
  lemma div_5: forall n:int. 0 <= n -> div n 5 >= 0
  lemma div_3: forall n:int. 0 <= n -> div n 3 >= 0

  lemma Closed_Formula:
    forall n:int. 0 <= n -> sum_multiple_3_5_lt (n+1) = closed_formula n

154 155 156 157
end

module Euler001

158 159 160
  use SumMultiple
  use int.Int
  use mach.int.Int
161

162 163 164
  let solve n
    requires { n >= 1 }
    ensures  { result = sum_multiple_3_5_lt n }
Mário Pereira's avatar
Mário Pereira committed
165 166 167 168
  = let n3 = (n-1) / 3 in
    let n5 = (n-1) / 5 in
    let n15 = (n-1) / 15 in
    (3 * n3 * (n3+1) + 5 * n5 * (n5+1) - 15 * n15 * (n15+1)) / 2
169

170 171 172 173 174
  (** Small test. Run it with

     why3 examples/euler001.mlw --exec Euler001.run

  *)
175
  let run () = solve 1000
MARCHE Claude's avatar
MARCHE Claude committed
176 177
    (* should return 233168 *)

178
  (** for the Why3 bench *)
179
  exception BenchFailure
MARCHE Claude's avatar
MARCHE Claude committed
180 181

  let bench () raises { BenchFailure -> true } =
182 183 184
    let x = run () in
    if x <> 233168 then raise BenchFailure;
    x
MARCHE Claude's avatar
MARCHE Claude committed
185

186 187
  (** for extraction *)

188
(*
189 190 191
  use string.Char
  use io.StdIO
  use ref.Ref
192

193
  let go ()
MARCHE Claude's avatar
MARCHE Claude committed
194
   ensures { !cur_linenum = (old !cur_linenum) + 1 }
195
  =
196 197 198 199 200 201
    print_char (chr 71); (* G *)
    print_char (chr 79); (* O *)
    print_char (chr 58); (* : *)
    print_char (chr 32); (*   *)
    print_int (run ());
    print_newline ()
202
*)
203

204
end