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
253
254
255
256
257
  (* 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
    | _ , _ -> let (a1,q1) = decomp 2 0 l1 l1 in
      let (a2,q2) = decomp 2 0 l2 l2 in
      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