examples/binomial.mlw

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