Euclid's algorithm

gcd.ml

open Core.Std

(* r = t * a + s * b *)
type step = { r : int; t : int; s : int }

let euler a b =
  let rec loop cur prev history =
    (* prev.r = q * cur.r + next.r *)
    let q = prev.r / cur.r in
    (* next.r = prev.r - q * cur.r
     *        = (prev.t * a + prev.s * b) - q * (cur.t * a + cur.s * b)
     *        = (prev.t - q * cur.t) * a + (prev.s - q * cur.s) * b      *)
    let next =
      { r = prev.r - q * cur.r
      ; t = prev.t - q * cur.t
      ; s = prev.s - q * cur.s
      }
    in
    assert (next.r = prev.r mod cur.r);
    assert (next.t * a + next.s * b = next.r);

    if next.r = 0
    then List.rev (cur::prev::history)
    else loop next cur (prev::history)
  in
  loop
    { r = b; t = 0; s = 1 }
    { r = a; t = 1; s = 0 }
    []

let rec print_steps = function
  | (i0::i1::i2::rest) ->
    printf "%5i = q * %-5i + %5i\t\t" i0.r i1.r i2.r;
    printf "%5i = %5i * a + %5i * b\n" i2.r i2.t i2.s;
    print_steps (i1::i2::rest)
  | [_;_] | [_] | [] -> ()

let print_euler a b =
  print_steps (euler a b)

let command =
  Command.basic
    ~summary:"euler"
    (let open Command.Spec in
     empty +> anon ("a" %: int) +> anon ("b" %: int)
    )
    (fun a b () -> print_euler a b)

let () = Command.run command