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 ->
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)

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
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] *)
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 ->
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 ->
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
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
104
    3 * tr n3 + 5 * tr n5 - 15 * tr n15
105

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

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:
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:
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:
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:
129 130
    forall n:int. n > 0 -> p (n-1) ->
      mod n 3 = 0 /\ mod n 5 = 0 -> p n
131

132
  constant b : int = 0
133

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

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