Toccata

(* Greatest common divisor, using the Euclidean algorithm *)

module EuclideanAlgorithm

  use mach.int.Int
  use number.Gcd

  let rec euclid (u v: int) : int
    variant  { v }
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
  =
    if v = 0 then
      u
    else
      euclid v (u % v)

end

module EuclideanAlgorithmIterative

  use mach.int.Int
  use ref.Ref
  use number.Gcd

  let euclid (u0 v0: int) : int
    requires { u0 >= 0 /\ v0 >= 0 }
    ensures  { result = gcd u0 v0 }
  =
    let ref u = u0 in
    let ref v = 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

  use mach.int.Int
  use number.Parity
  use number.Gcd

  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

  use number.Coprime

  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
    variant  { v, u }
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
  =
    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

  use mach.int.Int31
  use number.Gcd

  let rec euclid (u v: int31) : int31
    variant  { to_int v }
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
  =
    if v = 0 then
      u
    else
      euclid v (u % v)

end
*)

module EuclideanAlgorithm63

  use mach.int.Int63
  use number.Gcd

  let rec euclid (u v: int63) : int63
    variant  { to_int v }
    requires { u >= 0 /\ v >= 0 }
    ensures  { result = gcd u v }
  =
    if v = 0 then
      u
    else
      euclid v (u % v)

end