examples/fibonacci.mlw

   1
   2 theory FibonacciTest
   3
   4   use int.Fibonacci
   5
   6   lemma isfib_2_1 : fib 2 = 1
   7   lemma isfib_6_8 : fib 6 = 8
   8
   9   lemma not_isfib_2_2 : fib 2 <> 2
  10
  11 end
  12
  13 module FibonacciLinear
  14
  15   use int.Fibonacci
  16   use int.Int
  17   use ref.Ref
  18
  19   let fib (n:int) : int
  20     requires { n >= 0 }
  21     ensures { fib n = result}
  22   = let y = ref 0 in
  23     let x = ref 1 in
  24     for i = 0 to n - 1 do
  25       invariant { 0 <= i <= n /\ fib (i+1) = !x /\ fib i = !y }
  26       let aux = !y in
  27       y := !x;
  28       x := !x + aux
  29     done;
  30     !y
  31
  32 end
  33
  34 module FibonacciTailRecList
  35
  36   use int.Fibonacci
  37   use int.Int
  38   use int.Power
  39   use list.List
  40   use list.Mem
  41
  42   let rec ghost function sum_fib (l:list int) : int
  43     requires { forall n. mem n l -> n >= 0 }
  44   =
  45     match l with
  46     | Nil -> 0
  47     | Cons x r -> fib x + sum_fib r
  48     end
  49
  50   let rec ghost function sum_pow (l:list int) : int
  51     requires { forall n. mem n l -> n >= 0 }
  52     ensures { result >= 0 }
  53   =
  54     match l with
  55     | Nil -> 0
  56     | Cons x r -> assert { mem x l };
  57                   assert { forall n. mem n r -> mem n l };
  58                   power 2 x + sum_pow r
  59     end
  60
  61   lemma pow_pos : forall x. x >= 0 -> power 3 x > 0
  62
  63   let rec sum_fib_acc (acc:int) (l:list int) : int
  64     requires { forall n. mem n l -> n >= 0 }
  65     variant { sum_pow l }
  66     ensures { result = acc + sum_fib l }
  67   = match l with
  68     | Nil -> acc
  69     | Cons n r ->
  70       assert { n >= 0 };
  71       if n <= 1 then sum_fib_acc (acc + n) r
  72       else begin
  73       let l1 = Cons (n-2) r in
  74       assert { forall u. mem u l1 -> u >= 0 };
  75       let l2 = Cons (n-1) l1 in
  76       assert { forall u. mem u l2 -> u >= 0 };
  77       sum_fib_acc acc l2
  78       end
  79     end
  80
  81   let fib (n:int) : int
  82     requires { n >= 0 }
  83     ensures { result = fib n }
  84   =
  85     sum_fib_acc 0 (Cons n Nil)
  86
  87 end
  88
  89 (** {2 Recursive version, using ghost code} *)
  90
  91 module FibRecGhost
  92
  93   use int.Fibonacci
  94   use int.Int
  95
  96   let rec fib_aux (ghost n: int) (a b k: int) : int
  97     requires { k >= 0 }
  98     requires { 0 <= n && a = fib n && b = fib (n+1) }
  99     variant  { k }
 100     ensures  { result = fib (n+k) }
 101   = if k = 0 then a else fib_aux (n+1) b (a+b) (k-1)
 102
 103   let fib (n: int) : int
 104     requires { 0 <= n }
 105     ensures  { result = fib n }
 106   = fib_aux 0 0 1 n
 107
 108   let test42 () = fib 42
 109
 110   exception BenchFailure
 111
 112   let bench () raises { BenchFailure } =
 113     if test42 () <> 267914296 then raise BenchFailure
 114
 115 end
 116
 117 (** {2 Recursive version, without ghost code} *)
 118
 119 module FibRecNoGhost
 120
 121   use int.Fibonacci
 122   use int.Int
 123
 124   let rec fib_aux (a b k: int) : int
 125     requires { k >= 0 }
 126     requires { exists n: int. 0 <= n && a = fib n && b = fib (n+1) }
 127     variant  { k }
 128     ensures  { forall n: int. 0 <= n && a = fib n && b = fib (n+1) ->
 129                               result = fib (n+k) }
 130   = if k = 0 then a else fib_aux b (a+b) (k-1)
 131
 132   let fib (n: int) : int
 133     requires { 0 <= n }
 134     ensures  { result = fib n }
 135   = fib_aux 0 1 n
 136
 137 end
 138
 139 module SmallestFibAbove
 140
 141   use int.Fibonacci
 142   use int.Int
 143   use int.MinMax
 144   use ref.Ref
 145
 146   let smallest_fib_above (x: int) : int
 147     requires { 0 <= x }
 148     ensures  { exists k. 0 <= k /\ fib k <= x < fib (k+1) = result }
 149   =
 150     let a = ref 0 in
 151     let b = ref 1 in
 152     let ghost k = ref 0 in
 153     while !b <= x do
 154       invariant { 0 <= !k /\ !a = fib !k <= x /\ !b = fib (!k+1) }
 155       invariant { 0 <= !a /\ 1 <= !b }
 156       variant   { 2*x - (!a + !b) }
 157       let f = !a + !b in
 158       a := !b;
 159       b := f;
 160       k := !k + 1
 161    done;
 162    !b
 163
 164 end
 165
 166 (**
 167 Zeckendorf's theorem states that every positive integer can be
 168 represented uniquely as the sum of one or more distinct Fibonacci
 169 numbers in such a way that the sum does not include any two
 170 consecutive Fibonacci numbers.
 171
 172 Cf https://en.wikipedia.org/wiki/Zeckendorf%27s_theorem
 173 *)
 174
 175 module Zeckendorf
 176
 177   use int.Fibonacci
 178   use int.Int
 179   use int.MinMax
 180   use ref.Ref
 181   use list.List
 182   use SmallestFibAbove
 183
 184   function sum (l: list int) : int = match l with
 185   | Nil -> 0
 186   | Cons k r -> fib k + sum r
 187   end
 188
 189   (* sorted in increasing order, above min, and non consecutive *)
 190   predicate wf (min: int) (l: list int) = match l with
 191   | Nil -> true
 192   | Cons k r -> min <= k /\ wf (k + 2) r
 193   end
 194
 195   let rec lemma fib_nonneg (n: int) : unit
 196     requires { 0 <= n }
 197     ensures  { 0 <= fib n }
 198     variant  { n }
 199   = if n > 1 then begin fib_nonneg (n-2); fib_nonneg (n-1) end
 200
 201   let rec lemma fib_increasing (k1 k2: int) : unit
 202     requires { 0 <= k1 <= k2 }
 203     ensures  { fib k1 <= fib k2 }
 204     variant  { k2 - k1 }
 205   = if k1 < k2 then fib_increasing (k1+1) k2
 206
 207   let greatest_fib (x: int) : (int, int)
 208     requires { 0 < x }
 209     ensures  { let k, fk = result in
 210                2 <= k /\ 1 <= fk = fib k <= x < fib (k + 1) }
 211   =
 212     let a = ref 1 in
 213     let b = ref 1 in
 214     let k = ref 1 in
 215     while !b <= x do
 216       invariant { 1 <= !k /\ !a = fib !k <= x /\ !b = fib (!k + 1) }
 217       invariant { 1 <= !a /\ 1 <= !b }
 218       variant   { 2*x - (!a + !b) }
 219       let f = !a + !b in
 220       a := !b;
 221       b := f;
 222       k := !k + 1
 223    done;
 224    !k, !a
 225
 226   let zeckendorf (n: int) : list int
 227     requires { 0 <= n }
 228     ensures  { wf 2 result }
 229     ensures  { n = sum result }
 230   =
 231     let x = ref n in
 232     let l = ref Nil in
 233     while !x > 0 do
 234       invariant { 0 <= !x <= n }
 235       invariant { wf 2 !l }
 236       invariant { !x + sum !l = n }
 237       invariant { match !l with Nil -> true | Cons k _ -> !x < fib (k-1) end }
 238       variant   { !x }
 239       let k, fk = greatest_fib !x in
 240       x := !x - fk;
 241       l := Cons k !l
 242     done;
 243     !l
 244
 245   (* a more efficient, linear implementation *)
 246
 247   let zeckendorf_fast (n: int) : list int
 248     requires { 0 <= n }
 249     ensures  { wf 2 result }
 250     ensures  { n = sum result }
 251   =
 252     if n = 0 then Nil else
 253     let a = ref 1 in
 254     let b = ref 1 in
 255     let k = ref 1 in
 256     while !b <= n do
 257       invariant { 1 <= !k /\ !a = fib !k <= n /\ !b = fib (!k + 1) }
 258       invariant { 1 <= !a /\ 1 <= !b }
 259       variant   { 2*n - (!a + !b) }
 260       let f = !a + !b in
 261       a := !b;
 262       b := f;
 263       k := !k + 1
 264     done;
 265     assert { 2 <= !k /\ 1 <= !a = fib !k <= n < fib (!k + 1) = !b };
 266     let l = ref (Cons !k Nil) in
 267     let x = ref (n - !a) in
 268     while !x > 0 do
 269       invariant { 1 <= !k /\ !a = fib !k <= n /\ !x < !b = fib (!k + 1) }
 270       invariant { 1 <= !a /\ 1 <= !b }
 271       invariant { 0 <= !x <= n }
 272       invariant { wf 2 !l }
 273       invariant { !x + sum !l = n }
 274       invariant { match !l with Nil -> true | Cons k _ -> !x < fib (k-1) end }
 275       variant   { !k }
 276       if !a <= !x then begin
 277         x := !x - !a;
 278         l := Cons !k !l
 279       end;
 280       k := !k - 1;
 281       let tmp = !b - !a in
 282       b := !a;
 283       a := tmp
 284    done;
 285    !l
 286
 287   (* unicity proof *)
 288
 289   function snoc (l:list int) (x:int) : list int = match l with
 290     | Nil -> Cons x Nil
 291     | Cons y q -> Cons y (snoc q x)
 292     end
 293
 294   let rec lemma zeckendorf_unique (l1 l2:list int) : unit
 295     requires { wf 2 l1 /\ wf 2 l2 }
 296     requires { sum l1 = sum l2 }
 297     ensures { l1 = l2 }
 298     variant { sum l1 }
 299   = let rec decomp (k acc:int) (lc lb:list int) : (x: int, p: list int)
 300       requires { wf k lc }
 301       requires { k >= 2 /\ lc <> Nil }
 302       requires { 0 <= acc = sum lb - sum lc < fib (k-1) }
 303       ensures { fib x <= sum lb = acc + fib x + sum p < fib (x+1) }
 304       ensures { wf k p /\ x >= k /\ lc = snoc p x }
 305       variant { lc }
 306     = match lc with
 307       | Nil -> absurd
 308       | Cons x Nil -> x,Nil
 309       | Cons x q -> let y,p = decomp (x+2) (acc+fib x) q lb in y,Cons x p
 310       end in
 311     match l1 , l2 with
 312     | Nil , Nil -> ()
 313     | Nil , l | l , Nil -> let _ = decomp 2 0 l l in absurd
 314     | _ , _ -> let _,q1 = decomp 2 0 l1 l1 in
 315       let _,q2 = decomp 2 0 l2 l2 in
 316       zeckendorf_unique q1 q2
 317     end
 318
 319 end
 320
 321 (** {2 2x2 integer matrices} *)
 322
 323 theory Mat22
 324
 325   use int.Int
 326
 327   type t = { a11: int; a12: int; a21: int; a22: int }
 328
 329   constant id : t = { a11 = 1; a12 = 0; a21 = 0; a22 = 1 }
 330
 331   function mult (x: t) (y: t) : t =
 332     {
 333     a11 = x.a11 * y.a11 + x.a12 * y.a21; a12 = x.a11 * y.a12 + x.a12 * y.a22;
 334     a21 = x.a21 * y.a11 + x.a22 * y.a21; a22 = x.a21 * y.a12 + x.a22 * y.a22;
 335     }
 336
 337   clone export
 338     int.Exponentiation with
 339       type t = t, function one = id, function (*) = mult,
 340       goal Assoc, goal Unit_def_l, goal Unit_def_r,
 341       axiom . (* FIXME: replace with "goal" and prove *)
 342
 343 end
 344
 345 module FibonacciLogarithmic
 346
 347   use int.Int
 348   use int.Fibonacci
 349   use int.ComputerDivision
 350   use Mat22
 351
 352   val constant m1110 : t
 353     ensures { result = { a11 = 1; a12 = 1;
 354                          a21 = 1; a22 = 0 } }
 355
 356   (* computes ((1 1) (1 0))^n in O(log(n)) time
 357
 358      since it is a matrix of the shape ((a+b b) (b a)),
 359      we only return the pair (a, b) *)
 360
 361   let rec logfib (n:int) variant { n }
 362     requires { n >= 0 }
 363     ensures  { let a, b = result in
 364       power m1110 n = { a11 = a+b; a12 = b; a21 = b; a22 = a } }
 365   = if n = 0 then
 366       1, 0
 367     else begin
 368       assert { 0 <= div n 2 };
 369       let a, b = logfib (div n 2) in
 370       let c = a + b in
 371       if mod n 2 = 0 then
 372         begin
 373         assert { 2 * (div n 2) = (div n 2) + (div n 2) };
 374         a*a+ b*b, b*(a + c)
 375         end
 376       else
 377         begin
 378         assert { 2 * (div n 2) + 1 = (div n 2) + (div n 2) + 1 };
 379         b*(a + c), c*c + b*b
 380         end
 381     end
 382
 383   (* by induction, we easily prove that
 384
 385      (1 1)^n = (F(n+1) F(n)  )
 386      (1 0)     (F(n)   F(n-1))
 387
 388     thus, we can compute F(n) in O(log(n)) using funtion logfib above
 389   *)
 390
 391   let rec lemma fib_m (n: int)
 392     requires { n >= 0 }
 393     variant { n }
 394     ensures { let p = power m1110 n in fib (n+1) = p.a11 /\ fib n = p.a21 }
 395   = if n = 0 then () else fib_m (n-1)
 396
 397   let fibo n requires { n >= 0 } ensures { result = fib n } =
 398     let _, b = logfib n in b
 399
 400
 401   let test0 () = fibo 0
 402   let test1 () = fibo 1
 403   let test7 () = fibo 7
 404   let test42 () = fibo 42
 405   let test2014 () = fibo 2014
 406
 407   exception BenchFailure
 408
 409   let bench () raises { BenchFailure } =
 410     if test42 () <> 267914296 then raise BenchFailure;
 411     if test2014 () <> 3561413997540486142674781564382874188700994538849211456995042891654110985470076818421080236961243875711537543388676277339875963824466334432403730750376906026741819889036464401788232213002522934897299928844192803507157647764542466327613134605502785287441134627457615461304177503249289874066244145666889138852687147544158443155204157950294129177785119464446668374163746700969372438526182906768143740891051274219441912520127
 412     then raise BenchFailure
 413
 414 end