text_formalization/jordan/num_ext_gcd.hl

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