gcd_vc_sp.mlw 2.96 KB
Newer Older
MARCHE Claude's avatar
MARCHE Claude committed
1

2
(** Greatest common divisor, using the Euclidean algorithm *)
MARCHE Claude's avatar
MARCHE Claude committed
3 4 5

module EuclideanAlgorithm

6 7
  use mach.int.Int
  use number.Gcd
MARCHE Claude's avatar
MARCHE Claude committed
8 9 10 11 12

  let rec euclid (u v: int) : int
    variant  { v }
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
13
  = [@vc:sp]
MARCHE Claude's avatar
MARCHE Claude committed
14 15 16 17 18 19 20 21 22
    if v = 0 then
      u
    else
      euclid v (u % v)

end

module EuclideanAlgorithmIterative

23 24 25
  use mach.int.Int
  use ref.Ref
  use number.Gcd
MARCHE Claude's avatar
MARCHE Claude committed
26 27 28 29

  let euclid (u0 v0: int) : int
    requires { u0 >= 0 /\ v0 >= 0 }
    ensures  { result = gcd u0 v0 }
30
  = [@vc:sp]
MARCHE Claude's avatar
MARCHE Claude committed
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
    let u = ref u0 in
    let v = ref v0 in
    while !v <> 0 do
      invariant { !u >= 0 /\ !v >= 0 }
      invariant { gcd !u !v = gcd u0 v0 }
      variant   { !v }
      let tmp = !v in
      v := !u % !v;
      u := tmp
    done;
    !u

end

module BinaryGcd

47 48 49
  use mach.int.Int
  use number.Parity
  use number.Gcd
MARCHE Claude's avatar
MARCHE Claude committed
50 51 52 53 54

  lemma even1: forall n: int. 0 <= n -> even n <-> n = 2 * div n 2
  lemma odd1: forall n: int. 0 <= n -> not (even n) <-> n = 2 * div n 2 + 1
  lemma div_nonneg: forall n: int. 0 <= n -> 0 <= div n 2

55
  use number.Coprime
MARCHE Claude's avatar
MARCHE Claude committed
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75

  lemma gcd_even_even: forall u v: int. 0 <= v -> 0 <= u ->
    gcd (2 * u) (2 * v) = 2 * gcd u v
  lemma gcd_even_odd: forall u v: int. 0 <= v -> 0 <= u ->
    gcd (2 * u) (2 * v + 1) = gcd u (2 * v + 1)
  lemma gcd_even_odd2: forall u v: int. 0 <= v -> 0 <= u ->
    even u -> odd v -> gcd u v = gcd (div u 2) v
  lemma odd_odd_div2: forall u v: int. 0 <= v -> 0 <= u ->
    div ((2*u+1) - (2*v+1)) 2 = u - v

  let lemma gcd_odd_odd (u v: int)
    requires { 0 <= v <= u }
    ensures { gcd (2 * u + 1) (2 * v + 1) = gcd (u - v) (2 * v + 1) }
  = assert { gcd (2 * u + 1) (2 * v + 1) =
             gcd ((2*u+1) - 1*(2*v+1)) (2 * v + 1) }

  lemma gcd_odd_odd2: forall u v: int. 0 <= v <= u ->
    odd u -> odd v -> gcd u v = gcd (div (u - v) 2) v

  let rec binary_gcd (u v: int) : int
76
    variant  { v, u }
MARCHE Claude's avatar
MARCHE Claude committed
77 78
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
79
  = [@vc:sp]
MARCHE Claude's avatar
MARCHE Claude committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
    if v > u then binary_gcd v u else
    if v = 0 then u else
    if mod u 2 = 0 then
      if mod v 2 = 0 then 2 * binary_gcd (u / 2) (v / 2)
                     else binary_gcd (u / 2) v
      else
      if mod v 2 = 0 then binary_gcd u (v / 2)
                     else binary_gcd ((u - v) / 2) v

end

(** With machine integers.
    Note that we assume parameters u, v to be nonnegative.
    Otherwise, for u = v = min_int, the gcd could not be represented. *)

(* does not work with extraction driver ocaml64
module EuclideanAlgorithm31

98 99
  use mach.int.Int31
  use number.Gcd
MARCHE Claude's avatar
MARCHE Claude committed
100 101 102

  let rec euclid (u v: int31) : int31
    variant  { to_int v }
103 104
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
MARCHE Claude's avatar
MARCHE Claude committed
105
  =
106
    if v = 0 then
MARCHE Claude's avatar
MARCHE Claude committed
107 108 109 110 111 112 113 114 115
      u
    else
      euclid v (u % v)

end
*)

module EuclideanAlgorithm63

116 117
  use mach.int.Int63
  use number.Gcd
MARCHE Claude's avatar
MARCHE Claude committed
118 119 120

  let rec euclid (u v: int63) : int63
    variant  { to_int v }
121 122
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
123
  = [@vc:sp]
124
    if v = 0 then
MARCHE Claude's avatar
MARCHE Claude committed
125 126 127 128 129
      u
    else
      euclid v (u % v)

end