Jordan/num_ext_gcd.ml

   1 (*
   2         Author: Thomas C. Hales, 2003
   3
   4         GCD_CONV takes two HOL-light terms (NUMERALs) a and b and
   5         produces a theorem of the form
   6                 |- GCD a b = g
   7
   8         (In particular, the arguments cannot be negative.)
   9
  10 *)
  11
  12
  13 prioritize_num();;
  14
  15 let DIVIDE = new_definition(`DIVIDE a b = ?m. (b = m*a )`);;
  16
  17 parse_as_infix("||",(16,"right"));;
  18
  19 override_interface("||",`DIVIDE:num->num->bool`);;
  20
  21 (* Now prove the lemmas *)
  22
  23 let DIV_TAC t =   EVERY[ REP_GEN_TAC;
  24    REWRITE_TAC[DIVIDE];
  25    DISCH_ALL_TAC;
  26    REPEAT (FIRST_X_ASSUM CHOOSE_TAC);
  27    TRY (EXISTS_TAC t)];;
  28
  29
  30 let DIVIDE_DIVIDE = prove_by_refinement(
  31   `!a b c. (((a || b) /\ (b || c)) ==> (a || c))`,
  32    [
  33    DIV_TAC `m'*m`;
  34    ASM_REWRITE_TAC[MULT_ASSOC]
  35    ]);;
  36
  37 let DIVIDE_EQ = prove_by_refinement(
  38    `! a b. (((a || b) /\ (b || a)) ==> (a = b))`,
  39   [
  40   DIV_TAC `1`;
  41   FIRST_X_ASSUM (fun th -> (POP_ASSUM MP_TAC) THEN REWRITE_TAC[th]);
  42   ASM_CASES_TAC `b=0`;
  43   ASM_REWRITE_TAC[];
  44   ARITH_TAC;
  45   REWRITE_TAC[ARITH_RULE `(b = m*m'*b) = (1*b = m*m'*b)`];
  46   ASM_REWRITE_TAC[MULT_ASSOC;EQ_MULT_RCANCEL];
  47   DISCH_THEN (fun th -> MP_TAC (REWRITE_RULE[MULT_EQ_1] (GSYM th)) );
  48   DISCH_THEN (fun th -> REWRITE_TAC[CONJUNCT2 th] THEN ARITH_TAC);
  49   ]);;
  50
  51 let DIVIDE_SUM = prove_by_refinement(
  52   `!a b h. (((h || a) /\ (h||b)) ==> (h || (a+b)))`,
  53   [
  54   DIV_TAC `m+m'`;
  55   ASM_REWRITE_TAC[ARITH;RIGHT_ADD_DISTRIB];
  56   ]);;
  57
  58 let DIVIDE_SUMMAND = prove_by_refinement(
  59   `!a b h. (((h|| b) /\ (h || (a+b))) ==> (h|| a))`,
  60    [
  61    DIV_TAC `m'-m`;
  62    REWRITE_TAC[RIGHT_SUB_DISTRIB];
  63    REPEAT (FIRST_X_ASSUM  (fun th -> REWRITE_TAC[GSYM th]));
  64    ARITH_TAC;
  65    ]);;
  66
  67 let DIVIDE_PROD = prove_by_refinement(
  68    `!a b h. (((h|| a) ==> (h || (b*a))))`,
  69    [
  70    DIV_TAC `b*m`;
  71    ASM_REWRITE_TAC[MULT_ASSOC];
  72    ]);;
  73
  74 let DIVIDE_PROD2 = prove_by_refinement(
  75    `!a b h. (((h|| a) ==> (h || (a*b))))`,
  76    [
  77    DIV_TAC `b*m`;
  78    ASM_REWRITE_TAC[MULT_AC]
  79    ]);;
  80
  81 let GCD = new_definition(`GCD a b = @g.
  82         ((g || a) /\ (g || b) /\
  83         (!h. (((h || a) /\ (h || b)) ==> (h || g))))`);;
  84
  85 let gcd_certificate = prove(`!a b g. ((? r s r' s' a' b'.
  86         ((a = a'*g) /\ (b = b'*g) /\ (g +r'*a+s'*b= r*a + s*b)))
  87         ==> (GCD a b = g))`,
  88         let tac1 = (
  89         (REPEAT GEN_TAC)
  90         THEN (DISCH_TAC)
  91         THEN (REPEAT (POP_ASSUM CHOOSE_TAC))
  92         THEN (REWRITE_TAC[GCD])
  93         THEN (MATCH_MP_TAC SELECT_UNIQUE)
  94         THEN BETA_TAC
  95         THEN GEN_TAC
  96         THEN EQ_TAC) and
  97
  98         ygbranch = (
  99         DISCH_TAC
 100         THEN (MATCH_MP_TAC DIVIDE_EQ)
 101         THEN CONJ_TAC) and
 102
 103         ydivg_branch = (
 104         (SUBGOAL_TAC (` (y || (r*a + s*b))/\ (y || (r'*a +s'*b))`))
 105         THENL [((ASM MESON_TAC)[DIVIDE_SUM;DIVIDE_PROD]);
 106         ((ASM MESON_TAC)[DIVIDE_SUMMAND])]
 107         ) and
 108
 109         gdivy_branch = (
 110         (UNDISCH_TAC
 111           (`(y||a) /\ (y ||b) /\ (!h. (((h||a)/\(h||b))==> (h||y)))`))
 112         THEN (TAUT_TAC (` (A ==> B) ==> ((C /\ D/\ A)==> B)`))
 113         THEN (DISCH_TAC)
 114         THEN (POP_ASSUM MATCH_MP_TAC)
 115         THEN (REWRITE_TAC[DIVIDE])
 116         THEN (CONJ_TAC)
 117         THEN ((ASM MESON_TAC)[])
 118                 ) and
 119
 120         yghyp_branch = (
 121         (DISCH_TAC)
 122         THEN (let x t = REWRITE_TAC[t] in (POP_ASSUM x))
 123         THEN (CONJ_TAC)
 124         THENL [((ASM MESON_TAC)[DIVIDE]);ALL_TAC]
 125         THEN (CONJ_TAC)
 126         THENL [((ASM MESON_TAC)[DIVIDE]);ALL_TAC]
 127         THEN GEN_TAC
 128         THEN DISCH_TAC
 129         THEN (SUBGOAL_TAC (` (h || (r*a + s*b))/\ (h || (r'*a+s'*b))`))
 130         THENL [((ASM MESON_TAC)[DIVIDE_SUM;DIVIDE_PROD]);
 131                 ((ASM MESON_TAC)[DIVIDE_SUMMAND])]
 132                 ) in
 133         tac1 THENL [ygbranch THENL [ydivg_branch;gdivy_branch];yghyp_branch]);;
 134
 135 (* Now compute gcd with CAML num calculations,
 136    then check the answer in HOL-light *)
 137 let gcd_num x1 x2 =
 138         let rec gcd_data (a1,b1,x1,a2,b2,x2) =
 139         if (x1 < (Int 0)) then
 140                 gcd_data(minus_num a1,minus_num b1,minus_num x1,a2,b2,x2)
 141         else if (x2 < (Int 0)) then gcd_data(a1,b1,x1,minus_num a2,minus_num
 142         b2,minus_num x2)
 143         else if (x1 = (Int 0)) then (a2,b2,x2)
 144         else if (x1>x2) then gcd_data (a2,b2,x2,a1,b1,x1)
 145         else (
 146                 let r = (quo_num x2 x1) in
 147                 gcd_data (a1,b1,x1,a2 -/ r*/ a1,b2 -/ r*/ b1, x2 -/ r*/ x1)
 148              ) in
 149         gcd_data ((Int 1),(Int 0),x1,(Int 0),(Int 1),x2);;
 150
 151 let gcd_num x1 x2 =
 152         let rec gcd_data (a1,b1,x1,a2,b2,x2) =
 153         if (x1 < (Int 0)) then
 154                 gcd_data(minus_num a1,minus_num b1,minus_num x1,a2,b2,x2)
 155         else if (x2 < (Int 0)) then gcd_data(a1,b1,x1,minus_num a2,minus_num
 156         b2,minus_num x2)
 157         else if (x1 = (Int 0)) then (a2,b2,x2)
 158         else if (x1>x2) then gcd_data (a2,b2,x2,a1,b1,x1)
 159         else (
 160                 let r = (quo_num x2 x1) in
 161                 gcd_data (a1,b1,x1,a2 -/ r*/ a1,b2 -/ r*/ b1, x2 -/ r*/ x1)
 162              ) in
 163         gcd_data ((Int 1),(Int 0),x1,(Int 0),(Int 1),x2);;
 164
 165         (* g = gcd, (a',b') = (a,b)/g, g +r1'*a+s1'*b = r1*a+s1*b *)
 166 let gcd_numdata a b =
 167         let a = abs_num a in
 168         let b = abs_num b in
 169         let Z = Int 0 in
 170         let (r,s,g) = gcd_num a b in
 171         let a' = if (g=Z) then Z else round_num(a//g) in
 172         let b' = if (g=Z) then Z else round_num(b//g) in
 173         let _ = if not(a=a'*/g) then failwith "GCD_CONV a" else 0 in
 174         let _ = if not(b=b'*/g) then failwith "GCD_CONV b" else 0 in
 175         let _ = if not(g=r*/a+/s*/b) then failwith "GCD_CONV g" else 0 in
 176         let (r1,r1') = if (r >/ Z) then (r,Z) else (Z,minus_num r) in
 177         let (s1,s1') = if (s >/ Z) then (s,Z) else (Z,minus_num s) in
 178         (g,a,b,a',b',r1',s1',r1,s1);;
 179
 180 (* Here is the conversion.
 181         Example:
 182                 GCD_CONV (`66`) (`144`)
 183
 184 *)
 185 let GCD_CONV at bt =
 186         let a = dest_numeral at in
 187         let b = dest_numeral bt in
 188         let (g,a,b,a',b',r1',s1',r1,s1) = gcd_numdata a b in
 189         prove(parse_term("GCD "^(string_of_num a)^" "^(string_of_num b)^" = "^
 190                 (string_of_num g)),
 191                 (MATCH_MP_TAC gcd_certificate)
 192                 THEN (EXISTS_TAC (mk_numeral r1))
 193                 THEN (EXISTS_TAC (mk_numeral s1))
 194                 THEN (EXISTS_TAC (mk_numeral r1'))
 195                 THEN (EXISTS_TAC (mk_numeral s1'))
 196                 THEN (EXISTS_TAC (mk_numeral a'))
 197                 THEN (EXISTS_TAC (mk_numeral b'))
 198                 THEN (ARITH_TAC));;
 199
 200 (* Example:
 201         hol_gcd 66 144
 202
 203    This version can overflow on CAML integers before it reaches hol-light.
 204    Example:
 205         hol_gcd 1000000000000000000 10000000000000000000000
 206         - : thm = |- GCD 660865024 843055104 = 262144
 207 *)
 208
 209 let hol_gcd a b = GCD_CONV (mk_small_numeral a) (mk_small_numeral b);;
 210
 211 remove_interface ("||");;
 212 pop_priority();;
 213
 214
 215 (* test code *)
 216
 217 exception Test_suite_num_ext_gcd of string;;
 218
 219 (* For the tests we use integers a and b.  These can overflow if
 220    a and b are too large, so that we should confine ourselves to
 221    tests that are not too large.
 222 *)
 223
 224 let test_num_ext_gcd (a, b) =
 225   let a1 = string_of_int (abs a) in
 226   let b1 = string_of_int (abs b) in
 227   let c = gcd a b in
 228   let c1 = string_of_int (abs c) in
 229   let th = GCD_CONV (mk_small_numeral a) (mk_small_numeral b) in
 230   if (not (hyp th = ([]:term list))) then raise
 231     (failwith ("num_ext_gcd test suite failure "^a1^" "^b1))
 232   else if (not (concl th = (parse_term ("GCD "^a1^" "^b1^"="^c1))))
 233     then raise (failwith ("num_ext_gcd test suite failure "^a1^" "^b1))
 234   else ();;
 235
 236
 237 let test_suite_num_ext_gcd  =
 238   let _ =
 239     map test_num_ext_gcd
 240       [(0,0);(0,1);(1,0);(-0,-0);
 241        (2,3);(4,6);
 242        (0,2);(2,0);
 243        (10,100);(100,10);(17,100);(100,17)] in
 244    print_string "num_ext_gcd loaded\n";;
 245
 246 let divide = DIVIDE and
 247     gcd = GCD and
 248     gcd_conv = GCD_CONV;;
 249