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

2
theory Fibonacci "Fibonacci numbers"
MARCHE Claude's avatar
MARCHE Claude committed
3

4
  use export int.Int
MARCHE Claude's avatar
MARCHE Claude committed
5

Andrei Paskevich's avatar
Andrei Paskevich committed
6
  function fib int : int
MARCHE Claude's avatar
MARCHE Claude committed
7

8
9
  axiom fib0: fib 0 = 0
  axiom fib1: fib 1 = 1
10
  axiom fibn: forall n:int. n >= 2 -> fib n = fib (n-1) + fib (n-2)
11
12
13
14
15
16
17

end

theory FibonacciTest

  use import Fibonacci

18
  lemma isfib_2_1 : fib 2 = 1
19
20
21
  lemma isfib_6_8 : fib 6 = 8

  lemma not_isfib_2_2 : fib 2 <> 2
MARCHE Claude's avatar
MARCHE Claude committed
22
23
24
25
26
27
28

end

module FibonacciLinear

  use import Fibonacci
  use import int.Int
29
  use import ref.Ref
MARCHE Claude's avatar
MARCHE Claude committed
30

31
32
33
34
  let fib (n:int) : int
    requires { n >= 0 }
    ensures { fib n = result}
  = let y = ref 0 in
MARCHE Claude's avatar
MARCHE Claude committed
35
    let x = ref 1 in
36
    for i = 0 to n - 1 do
Andrei Paskevich's avatar
Andrei Paskevich committed
37
      invariant { 0 <= i <= n /\ fib (i+1) = !x /\ fib i = !y }
MARCHE Claude's avatar
MARCHE Claude committed
38
      let aux = !y in
39
      y := !x;
MARCHE Claude's avatar
MARCHE Claude committed
40
41
42
43
44
      x := !x + aux
    done;
    !y

end
45

46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
module FibRecGhost "recursive version, using ghost code"

  use import Fibonacci
  use import int.Int

  let rec fib_aux (ghost n: int) (a b k: int) : int
    requires { 0 <= n && a = fib n && b = fib (n+1) }
    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

end

module FibRecNoGhost "recursive version, without ghost code"

  use import Fibonacci
  use import int.Int

  let rec fib_aux (a b k: int) : int
    requires { exists n: int. 0 <= n && a = fib n && b = fib (n+1) }
    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

81
82
83
84
theory Mat22 "2x2 integer matrices"

  use import int.Int

85
  type t = { a11: int; a12: int; a21: int; a22: int }
86

87
  constant id : t = { a11 = 1; a12 = 0; a21 = 0; a22 = 1 }
88

Andrei Paskevich's avatar
Andrei Paskevich committed
89
  function mult (x: t) (y: t) : t =
90
    {
91
    a11 = x.a11 * y.a11 + x.a12 * y.a21; a12 = x.a11 * y.a12 + x.a12 * y.a22;
92
    a21 = x.a21 * y.a11 + x.a22 * y.a21; a22 = x.a21 * y.a12 + x.a22 * y.a22;
93
    }
94

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

98
  clone export
Andrei Paskevich's avatar
Andrei Paskevich committed
99
    int.Exponentiation with type t = t, function one = id, function (*) = mult
100
101
102
103
104
105
106

end

module FibonacciLogarithmic

  use import Fibonacci
  use import int.EuclideanDivision
107
  use import Mat22
108

109
110
  constant m1110 : t = { a11 = 1; a12 = 1;
                         a21 = 1; a22 = 0 }
111

112
113
  (* computes ((1 1) (1 0))^n in O(log(n)) time

114
     since it is a matrix of the shape ((a+b b) (b a)),
115
116
     we only return the pair (a, b) *)

117
118
119
120
121
  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
122
      (1, 0)
123
124
125
126
    else begin
      let a, b = logfib (div n 2) in
      let c = a + b in
      if mod n 2 = 0 then
127
        (a*a + b*b, b*(a + c))
128
      else
129
        (b*(a + c), c*c + b*b)
130
131
    end

132
133
134
135
136
137
138
139
  (* 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
  *)

140
  lemma fib_m :
141
    forall n: int. n >= 0 ->
Andrei Paskevich's avatar
Andrei Paskevich committed
142
    let p = power m1110 n in fib (n+1) = p.a11 /\ fib n = p.a21
143

144
  let fibo n requires { n >= 0 } ensures { result = fib n } =
145
146
    let _, b = logfib n in b

147
148
149
150
151
152

  let test0 () = fibo 0
  let test1 () = fibo 1
  let test7 () = fibo 7
  let test42 () = fibo 42

153
end