1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
(** Binomial coefficients
The binomial coefficient C(n,k) is equal to
n*(n-1)*(n-2)*...*(n-k+1)
-------------------------
k*(k-1)*(k-2)*...*1
This can be readily computed with three lines of C:
int c = 1;
for (int i = 1; i <= k ; i++)
c = c * (n - i + 1) / i;
In the code above, it is not obvious that each division is exact.
Below is a proof.
Author: Jean-Christophe Filliâtre (CNRS)
*)
use int.Int
use int.EuclideanDivision
use int.MinMax
let function (/) (x: int) (y: int)
requires { [@expl:check division by zero] y <> 0 }
= div x y
let rec function comb (n k: int) : int
requires { 0 <= k <= n }
variant { n }
ensures { result >= 1 }
= if k = 0 || k = n then 1 else comb (n-1) k + comb (n-1) (k-1)
let rec lemma prop1 (n k: int)
requires { 0 <= k <= n }
ensures { comb n (n - k) = comb n k }
variant { n }
= if 0 < k < n then (prop1 (n-1) k; prop1 (n-1) (k-1))
let rec lemma prop2 (n k: int)
requires { 1 <= k <= n }
ensures { k * comb n k = comb n (k - 1) * (n - k + 1) }
variant { n }
= if k < n then prop2 (n-1) k;
if 1 < k then prop2 (n-1) (k-1)
let compute (n k: int) : (r: int)
requires { 0 <= k <= n }
ensures { r = comb n k }
= let k = min k (n - k) in
let ref r = 1 in
for i = 1 to k do
invariant { 1 <= i <= k + 1 }
invariant { r = comb n (i - 1) }
r <- r * (n - i + 1) / i;
prop2 n i;
done;
r
|
