fibonacci.mlw 4.19 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
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
52
    requires { k >= 0 }
53
    requires { 0 <= n && a = fib n && b = fib (n+1) }
54
    variant  { k }
55
56
57
58
59
60
61
62
    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
63
64
65
66
67
68
69
  let test42 () = fib 42

  exception BenchFailure

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

70
71
72
73
74
75
76
77
end

module FibRecNoGhost "recursive version, without ghost code"

  use import Fibonacci
  use import int.Int

  let rec fib_aux (a b k: int) : int
78
    requires { k >= 0 }
79
    requires { exists n: int. 0 <= n && a = fib n && b = fib (n+1) }
80
    variant  { k }
81
82
83
84
85
86
87
88
89
90
91
    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

92
93
94
95
theory Mat22 "2x2 integer matrices"

  use import int.Int

96
  type t = { a11: int; a12: int; a21: int; a22: int }
97

98
  constant id : t = { a11 = 1; a12 = 0; a21 = 0; a22 = 1 }
99

Andrei Paskevich's avatar
Andrei Paskevich committed
100
  function mult (x: t) (y: t) : t =
101
    {
102
    a11 = x.a11 * y.a11 + x.a12 * y.a21; a12 = x.a11 * y.a12 + x.a12 * y.a22;
103
    a21 = x.a21 * y.a11 + x.a22 * y.a21; a22 = x.a21 * y.a12 + x.a22 * y.a22;
104
    }
105

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

109
  clone export
Andrei Paskevich's avatar
Andrei Paskevich committed
110
    int.Exponentiation with type t = t, function one = id, function (*) = mult
111
112
113
114
115
116
117

end

module FibonacciLogarithmic

  use import Fibonacci
  use import int.EuclideanDivision
118
  use import Mat22
119

120
121
  constant m1110 : t = { a11 = 1; a12 = 1;
                         a21 = 1; a22 = 0 }
122

123
124
  (* computes ((1 1) (1 0))^n in O(log(n)) time

125
     since it is a matrix of the shape ((a+b b) (b a)),
126
127
     we only return the pair (a, b) *)

128
129
130
131
132
  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
133
      (1, 0)
134
135
136
137
    else begin
      let a, b = logfib (div n 2) in
      let c = a + b in
      if mod n 2 = 0 then
138
        (a*a + b*b, b*(a + c))
139
      else
140
        (b*(a + c), c*c + b*b)
141
142
    end

143
144
145
146
147
148
149
150
  (* 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
  *)

151
  lemma fib_m :
152
    forall n: int. n >= 0 ->
Andrei Paskevich's avatar
Andrei Paskevich committed
153
    let p = power m1110 n in fib (n+1) = p.a11 /\ fib n = p.a21
154

155
  let fibo n requires { n >= 0 } ensures { result = fib n } =
156
157
    let _, b = logfib n in b

158
159
160
161
162

  let test0 () = fibo 0
  let test1 () = fibo 1
  let test7 () = fibo 7
  let test42 () = fibo 42
163
  let test2014 () = fibo 2014
164

MARCHE Claude's avatar
MARCHE Claude committed
165
166
167
168
  exception BenchFailure

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

172
end