fibonacci.mlw 8.78 KB
Newer Older
MARCHE Claude's avatar
MARCHE Claude committed
1

2 3
theory FibonacciTest

4
  use import int.Fibonacci
5

6
  lemma isfib_2_1 : fib 2 = 1
7 8 9
  lemma isfib_6_8 : fib 6 = 8

  lemma not_isfib_2_2 : fib 2 <> 2
MARCHE Claude's avatar
MARCHE Claude committed
10 11 12 13 14

end

module FibonacciLinear

15
  use import int.Fibonacci
MARCHE Claude's avatar
MARCHE Claude committed
16
  use import int.Int
17
  use import ref.Ref
MARCHE Claude's avatar
MARCHE Claude committed
18

19 20 21 22
  let fib (n:int) : int
    requires { n >= 0 }
    ensures { fib n = result}
  = let y = ref 0 in
MARCHE Claude's avatar
MARCHE Claude committed
23
    let x = ref 1 in
24
    for i = 0 to n - 1 do
Andrei Paskevich's avatar
Andrei Paskevich committed
25
      invariant { 0 <= i <= n /\ fib (i+1) = !x /\ fib i = !y }
MARCHE Claude's avatar
MARCHE Claude committed
26
      let aux = !y in
27
      y := !x;
MARCHE Claude's avatar
MARCHE Claude committed
28 29 30 31 32
      x := !x + aux
    done;
    !y

end
33

34 35
module FibRecGhost "recursive version, using ghost code"

36
  use import int.Fibonacci
37 38 39
  use import int.Int

  let rec fib_aux (ghost n: int) (a b k: int) : int
40
    requires { k >= 0 }
41
    requires { 0 <= n && a = fib n && b = fib (n+1) }
42
    variant  { k }
43 44 45 46 47 48 49 50
    ensures  { result = fib (n+k) }
  = if k = 0 then a else fib_aux (n+1) b (a+b) (k-1)

  let fib (n: int) : int
    requires { 0 <= n }
    ensures  { result = fib n }
  = fib_aux 0 0 1 n

MARCHE Claude's avatar
MARCHE Claude committed
51 52 53 54 55 56 57
  let test42 () = fib 42

  exception BenchFailure

  let bench () raises { BenchFailure -> true } =
    if test42 () <> 267914296 then raise BenchFailure

58 59 60 61
end

module FibRecNoGhost "recursive version, without ghost code"

62
  use import int.Fibonacci
63 64 65
  use import int.Int

  let rec fib_aux (a b k: int) : int
66
    requires { k >= 0 }
67
    requires { exists n: int. 0 <= n && a = fib n && b = fib (n+1) }
68
    variant  { k }
69 70 71 72 73 74 75 76 77 78 79
    ensures  { forall n: int. 0 <= n && a = fib n && b = fib (n+1) ->
                              result = fib (n+k) }
  = if k = 0 then a else fib_aux b (a+b) (k-1)

  let fib (n: int) : int
    requires { 0 <= n }
    ensures  { result = fib n }
  = fib_aux 0 1 n

end

80 81 82 83 84 85 86 87 88
module SmallestFibAbove

  use import int.Fibonacci
  use import int.Int
  use import int.MinMax
  use import ref.Ref

  let smallest_fib_above (x: int) : int
    requires { 0 <= x }
89
    ensures  { exists k: int. 0 <= k /\ fib k <= x < fib (k+1) = result }
90 91 92 93 94
  =
    let a = ref 0 in
    let b = ref 1 in
    while !b <= x do
      invariant { exists k: int. 0 <= k /\ !a = fib k <= x /\ !b = fib (k+1) }
95
      invariant { 0 <= !a /\ 1 <= !b }
96 97 98 99 100 101 102 103 104
      variant   { 2*x - (!a + !b) }
      let f = !a + !b in
      a := !b;
      b := f
   done;
   !b

end

105 106 107 108 109 110 111 112 113
(**
Zeckendorf's theorem states that every positive integer can be
represented uniquely as the sum of one or more distinct Fibonacci
numbers in such a way that the sum does not include any two
consecutive Fibonacci numbers.

Cf https://en.wikipedia.org/wiki/Zeckendorf%27s_theorem
*)

114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
module Zeckendorf

  use import int.Fibonacci
  use import int.Int
  use import int.MinMax
  use import ref.Ref
  use import list.List
  use import SmallestFibAbove

  function sum (l: list int) : int = match l with
  | Nil -> 0
  | Cons k r -> fib k + sum r
  end

  (* sorted in increasing order, above min, and non consecutive *)
  predicate wf (min: int) (l: list int) = match l with
  | Nil -> true
  | Cons k r -> min <= k /\ wf (k + 2) r
  end

  let rec lemma fib_nonneg (n: int) : unit
    requires { 0 <= n }
    ensures  { 0 <= fib n }
    variant  { n }
  = if n > 1 then begin fib_nonneg (n-2); fib_nonneg (n-1) end

  let rec lemma fib_increasing (k1 k2: int) : unit
    requires { 0 <= k1 <= k2 }
    ensures  { fib k1 <= fib k2 }
    variant  { k2 - k1 }
  = if k1 < k2 then fib_increasing (k1+1) k2

  let greatest_fib (x: int) : (int, int)
    requires { 0 < x }
    ensures  { let (k, fk) = result in
               2 <= k /\ 1 <= fk = fib k <= x < fib (k + 1) }
  =
    let a = ref 1 in
    let b = ref 1 in
    let k = ref 1 in
    while !b <= x do
      invariant { 1 <= !k /\ !a = fib !k <= x /\ !b = fib (!k + 1) }
      invariant { 1 <= !a /\ 1 <= !b }
      variant   { 2*x - (!a + !b) }
      let f = !a + !b in
      a := !b;
      b := f;
      k := !k + 1
   done;
   (!k, !a)

  let zeckendorf (n: int) : list int
    requires { 0 <= n }
    ensures  { wf 2 result }
    ensures  { n = sum result }
  =
    let x = ref n in
    let l = ref Nil in
    while !x > 0 do
      invariant { 0 <= !x <= n }
      invariant { wf 2 !l }
      invariant { !x + sum !l = n }
      invariant { match !l with Nil -> true | Cons k _ -> !x < fib (k-1) end }
      variant   { !x }
      let (k, fk) = greatest_fib !x in
      x := !x - fk;
      l := Cons k !l
    done;
    !l

184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
  (* a more efficient, linear implementation *)

  let zeckendorf_fast (n: int) : list int
    requires { 0 <= n }
    ensures  { wf 2 result }
    ensures  { n = sum result }
  =
    if n = 0 then Nil else
    let a = ref 1 in
    let b = ref 1 in
    let k = ref 1 in
    while !b <= n do
      invariant { 1 <= !k /\ !a = fib !k <= n /\ !b = fib (!k + 1) }
      invariant { 1 <= !a /\ 1 <= !b }
      variant   { 2*n - (!a + !b) }
      let f = !a + !b in
      a := !b;
      b := f;
      k := !k + 1
    done;
    assert { 2 <= !k /\ 1 <= !a = fib !k <= n < fib (!k + 1) = !b };
    let l = ref (Cons !k Nil) in
    let x = ref (n - !a) in
    while !x > 0 do
      invariant { 1 <= !k /\ !a = fib !k <= n /\ !x < !b = fib (!k + 1) }
      invariant { 1 <= !a /\ 1 <= !b }
      invariant { 0 <= !x <= n }
      invariant { wf 2 !l }
      invariant { !x + sum !l = n }
      invariant { match !l with Nil -> true | Cons k _ -> !x < fib (k-1) end }
      variant   { !k }
      if !a <= !x then begin
        x := !x - !a;
        l := Cons !k !l
      end;
      k := !k - 1;
      let tmp = !b - !a in
      b := !a;
      a := tmp
   done;
   !l
225

226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
  (* unicity proof *)

  function snoc (l:list int) (x:int) : list int = match l with
    | Nil -> Cons x Nil
    | Cons y q -> Cons y (snoc q x)
    end

  let rec lemma zeckendorf_unique (l1 l2:list int) : unit
    requires { wf 2 l1 /\ wf 2 l2 }
    requires { sum l1 = sum l2 }
    ensures { l1 = l2 }
    variant { sum l1 }
  = let rec decomp (k acc:int) (lc lb:list int) : (int,list int)
      requires { wf k lc }
      requires { k >= 2 /\ lc <> Nil }
      requires { 0 <= acc = sum lb - sum lc < fib (k-1) }
      returns { (x,p) -> fib x <= sum lb = acc + fib x + sum p < fib (x+1) }
      returns { (x,p) -> wf k p /\ x >= k /\ lc = snoc p x }
      variant { lc }
    = match lc with
      | Nil -> absurd
      | Cons x Nil -> (x,Nil)
      | Cons x q -> let (y,p) = decomp (x+2) (acc+fib x) q lb in (y,Cons x p)
      end in
    match l1 , l2 with
    | Nil , Nil -> ()
    | Nil , l | l , Nil -> let _ = decomp 2 0 l l in absurd
253 254
    | _ , _ -> let (_,q1) = decomp 2 0 l1 l1 in
      let (_,q2) = decomp 2 0 l2 l2 in
255 256 257
      zeckendorf_unique q1 q2
    end

258 259
end

260 261 262 263
theory Mat22 "2x2 integer matrices"

  use import int.Int

264
  type t = { a11: int; a12: int; a21: int; a22: int }
265

266
  constant id : t = { a11 = 1; a12 = 0; a21 = 0; a22 = 1 }
267

Andrei Paskevich's avatar
Andrei Paskevich committed
268
  function mult (x: t) (y: t) : t =
269
    {
270
    a11 = x.a11 * y.a11 + x.a12 * y.a21; a12 = x.a11 * y.a12 + x.a12 * y.a22;
271
    a21 = x.a21 * y.a11 + x.a22 * y.a21; a22 = x.a21 * y.a12 + x.a22 * y.a22;
272
    }
273

274
  (* holds, but not useful *)
Andrei Paskevich's avatar
Andrei Paskevich committed
275
  (* clone algebra.Assoc with type t = t, function op = mult, lemma Assoc *)
276

277
  clone export
Andrei Paskevich's avatar
Andrei Paskevich committed
278
    int.Exponentiation with type t = t, function one = id, function (*) = mult
279 280 281 282 283

end

module FibonacciLogarithmic

284
  use import int.Int
285
  use import int.Fibonacci
286
  use import int.EuclideanDivision
287
  use import Mat22
288

289 290
  constant m1110 : t = { a11 = 1; a12 = 1;
                         a21 = 1; a22 = 0 }
291

292 293
  (* computes ((1 1) (1 0))^n in O(log(n)) time

294
     since it is a matrix of the shape ((a+b b) (b a)),
295 296
     we only return the pair (a, b) *)

297 298 299 300 301
  let rec logfib n variant { n }
    requires { n >= 0 }
    ensures  { let a, b = result in
      power m1110 n = { a11 = a+b; a12 = b; a21 = b; a22 = a } }
  = if n = 0 then
302
      (1, 0)
303 304 305 306
    else begin
      let a, b = logfib (div n 2) in
      let c = a + b in
      if mod n 2 = 0 then
307
        (a*a + b*b, b*(a + c))
308
      else
309
        (b*(a + c), c*c + b*b)
310 311
    end

312 313 314 315 316 317 318 319
  (* by induction, we easily prove that

     (1 1)^n = (F(n+1) F(n)  )
     (1 0)     (F(n)   F(n-1))

    thus, we can compute F(n) in O(log(n)) using funtion logfib above
  *)

320
  lemma fib_m :
321
    forall n: int. n >= 0 ->
Andrei Paskevich's avatar
Andrei Paskevich committed
322
    let p = power m1110 n in fib (n+1) = p.a11 /\ fib n = p.a21
323

324
  let fibo n requires { n >= 0 } ensures { result = fib n } =
325 326
    let _, b = logfib n in b

327 328 329 330 331

  let test0 () = fibo 0
  let test1 () = fibo 1
  let test7 () = fibo 7
  let test42 () = fibo 42
332
  let test2014 () = fibo 2014
333

MARCHE Claude's avatar
MARCHE Claude committed
334 335 336 337
  exception BenchFailure

  let bench () raises { BenchFailure -> true } =
    if test42 () <> 267914296 then raise BenchFailure;
MARCHE Claude's avatar
MARCHE Claude committed
338
    if test2014 () <> 3561413997540486142674781564382874188700994538849211456995042891654110985470076818421080236961243875711537543388676277339875963824466334432403730750376906026741819889036464401788232213002522934897299928844192803507157647764542466327613134605502785287441134627457615461304177503249289874066244145666889138852687147544158443155204157950294129177785119464446668374163746700969372438526182906768143740891051274219441912520127
MARCHE Claude's avatar
MARCHE Claude committed
339
    then raise BenchFailure
MARCHE Claude's avatar
MARCHE Claude committed
340

341
end