100/gcd.ml

   1 (* ========================================================================= *)
   2 (* Euclidean GCD algorithm.                                                  *)
   3 (* ========================================================================= *)
   4
   5 needs "Library/prime.ml";;
   6
   7 let egcd = define
   8  `egcd(m,n) = if m = 0 then n
   9               else if n = 0 then m
  10               else if m <= n then egcd(m,n - m)
  11               else egcd(m - n,n)`;;
  12
  13 (* ------------------------------------------------------------------------- *)
  14 (* Main theorems.                                                            *)
  15 (* ------------------------------------------------------------------------- *)
  16
  17 let EGCD_INVARIANT = prove
  18  (`!m n d. d divides egcd(m,n) <=> d divides m /\ d divides n`,
  19   GEN_TAC THEN GEN_TAC THEN WF_INDUCT_TAC `m + n` THEN
  20   GEN_TAC THEN ONCE_REWRITE_TAC[egcd] THEN
  21   ASM_CASES_TAC `m = 0` THEN ASM_REWRITE_TAC[DIVIDES_0] THEN
  22   ASM_CASES_TAC `n = 0` THEN ASM_REWRITE_TAC[DIVIDES_0] THEN
  23   COND_CASES_TAC THEN
  24   (W(fun (asl,w) -> FIRST_X_ASSUM(fun th ->
  25     MP_TAC(PART_MATCH (lhs o snd o dest_forall o rand) th (lhand w)))) THEN
  26    ANTS_TAC THENL [ASM_ARITH_TAC; ALL_TAC]) THEN
  27   ASM_MESON_TAC[DIVIDES_SUB; DIVIDES_ADD; SUB_ADD; LE_CASES]);;
  28
  29 (* ------------------------------------------------------------------------- *)
  30 (* Hence we get the proper behaviour, and it's equal to the real GCD.        *)
  31 (* ------------------------------------------------------------------------- *)
  32
  33 let EGCD_GCD = prove
  34  (`!m n. egcd(m,n) = gcd(m,n)`,
  35   ONCE_REWRITE_TAC[GSYM GCD_UNIQUE] THEN
  36   MESON_TAC[EGCD_INVARIANT; DIVIDES_REFL]);;
  37
  38 let EGCD = prove
  39  (`!a b. (egcd (a,b) divides a /\ egcd (a,b) divides b) /\
  40          (!e. e divides a /\ e divides b ==> e divides egcd (a,b))`,
  41   REWRITE_TAC[EGCD_GCD; GCD]);;