Rqe/poly_ext.ml

   1 let poly_tm = `poly`;;
   2
   3 let dest_poly tm =
   4   let poly,[l;var] = strip_ncomb 2 tm in
   5   if not (poly = poly_tm) then failwith "dest_poly: not a poly"
   6   else l,var;;
   7
   8 let is_poly tm = fst (strip_comb tm) = `poly`;;
   9
  10 (* ------------------------------------------------------------------------- *)
  11 (* Get the lead variable in polynomial; &1 if a constant.                    *)
  12 (* ------------------------------------------------------------------------- *)
  13
  14 let polyvar =
  15   let dummy_tm = `&1` in
  16   fun tm -> if is_ratconst tm then dummy_tm else lhand(rand tm);;
  17
  18
  19 (*
  20 let k00 = `&3 * x * y pow 2 + &2 * x pow 2 * y * z + z * x + &3 * y * z`
  21 let k0 = `(&0 + y * (&0 + z * &3)) + x * (((&0 + z * &1) + y * (&0 + y * &3)) +  x * (&0 + y * (&0 + z * &2)))`;;
  22 # polyvar k0;;
  23 val it : Term.term = `x`
  24 *)
  25
  26 (* ---------------------------------------------------------------------- *)
  27 (*  Is a constant polynomial (wrt variable ordering)                      *)
  28 (* ---------------------------------------------------------------------- *)
  29
  30 let is_constant vars p =
  31   assert (not (vars = []));
  32   try
  33     let l,r = dest_plus p in
  34     let x,r2 = dest_mult r in
  35       if x = hd vars then false else true
  36   with _ ->
  37     if p = hd vars then false else true;;
  38
  39 (* ------------------------------------------------------------------------- *)
  40 (* We only use this as a handy way to do derivatives.                        *)
  41 (* ------------------------------------------------------------------------- *)
  42
  43 let POLY = prove
  44  (`(poly [] x = &0) /\
  45    (poly [__c__] x = __c__) /\
  46    (poly (CONS __h__ __t__) x = __h__ + x * poly __t__ x)`,
  47    REWRITE_TAC[poly] THEN REAL_ARITH_TAC);;
  48
  49 (* ------------------------------------------------------------------------- *)
  50 (* Convert in and out of list representations.                               *)
  51 (* ------------------------------------------------------------------------- *)
  52
  53 (* THIS IS BAD CODE!!!  It depends on the names of the variables in POLY *)
  54 let POLY_ENLIST_CONV vars =
  55   let lem = GEN rx POLY in
  56   let [cnv_0; cnv_1; cnv_2] =
  57     map (fun th -> GEN_REWRITE_CONV I [GSYM th]) (CONJUNCTS (ISPEC (hd vars) lem))
  58   and zero_tm = rzero in
  59   let rec conv tm =
  60     if polyvar tm = hd vars then
  61       (funpow 2 RAND_CONV conv THENC cnv_2) tm
  62     else if tm = zero_tm then cnv_0 tm
  63     else cnv_1 tm in
  64     conv;;
  65
  66
  67 (*
  68 map GSYM (CONJUNCTS (ISPEC (hd vars) lem))
  69
  70 POLY_ENLIST_CONV vars p in
  71
  72 let tm = `&0 + c * &1`
  73
  74
  75
  76 POLY_ENLIST_CONV vars tm
  77
  78 #trace conv
  79 POLY_ENLIST_CONV vars tm
  80 let vars = [ry;rx]
  81 let tm =  `&0 + y * (&0 + x * &1)`
  82
  83
  84 let k1 = rhs(concl (POLY_ENLIST_CONV [`x:real`;`y:real`;`z:real`] k0));;
  85 POLY_ENLIST_CONV [`x:real`;`y:real`;`z:real`] k0;;
  86 val it : Hol.thm =
  87   |- k0 =
  88      poly [&0 + y * (&0 + z * &3);
  89            &0 * z * &1 + y * (&0 + y * &3);
  90            &0 + y * (&0 + z * &2)] x
  91 *)
  92
  93 let POLY_DELIST_CONV =
  94   let [cnv_0; cnv_1; cnv_2] =
  95     map (fun th -> GEN_REWRITE_CONV I [th]) (CONJUNCTS POLY) in
  96   let rec conv tm =
  97     (cnv_0 ORELSEC cnv_1 ORELSEC (cnv_2 THENC funpow 2 RAND_CONV conv)) tm in
  98   conv;;
  99
 100 (*
 101 # POLY_DELIST_CONV `poly [&5; &6; &7] x`;;
 102 val it : Hol.thm = |- poly [&5; &6; &7] x = &5 + x * (&6 + x * &7)
 103 *)
 104
 105 (* ------------------------------------------------------------------------- *)
 106 (* Differentiation within list representation.                               *)
 107 (* ------------------------------------------------------------------------- *)
 108
 109 (* let poly_diff_aux = new_recursive_definition list_RECURSION *)
 110 (*   `(poly_diff_aux n [] = []) /\ *)
 111 (*    (poly_diff_aux n (CONS h t) = CONS (&n * h) (poly_diff_aux (SUC n) t))`;; *)
 112
 113 (* let poly_diff = new_definition *)
 114 (*   `poly_diff l = if l = [] then [] else poly_diff_aux 1 (TL l)`;; *)
 115
 116 let POLY_DIFF_CLAUSES = prove
 117  (`(poly_diff [] = []) /\
 118    (poly_diff [c] = []) /\
 119    (poly_diff (CONS h t) = poly_diff_aux 1 t)`,
 120   REWRITE_TAC[poly_diff; NOT_CONS_NIL; HD; TL; poly_diff_aux]);;
 121
 122 let POLY_DIFF_LEMMA = prove
 123  (`!l n x. ((\x. (x pow (SUC n)) * poly l x) diffl
 124                    ((x pow n) * poly (poly_diff_aux (SUC n) l) x))(x)`,
 125 (* {{{ Proof *)
 126
 127   LIST_INDUCT_TAC THEN
 128   REWRITE_TAC[poly; poly_diff_aux; REAL_MUL_RZERO; DIFF_CONST] THEN
 129   MAP_EVERY X_GEN_TAC [`n:num`; `x:real`] THEN
 130   REWRITE_TAC[REAL_LDISTRIB; REAL_MUL_ASSOC] THEN
 131   ONCE_REWRITE_TAC[GSYM(ONCE_REWRITE_RULE[REAL_MUL_SYM] (CONJUNCT2 pow))] THEN
 132   POP_ASSUM(MP_TAC o SPECL [`SUC n`; `x:real`]) THEN
 133   SUBGOAL_THEN `(((\x. (x pow (SUC n)) * h)) diffl
 134                         ((x pow n) * &(SUC n) * h))(x)`
 135   (fun th -> DISCH_THEN(MP_TAC o CONJ th)) THENL
 136    [REWRITE_TAC[REAL_MUL_ASSOC] THEN ONCE_REWRITE_TAC[REAL_MUL_SYM] THEN
 137     MP_TAC(SPEC `\x. x pow (SUC n)` DIFF_CMUL) THEN BETA_TAC THEN
 138     DISCH_THEN MATCH_MP_TAC THEN
 139     MP_TAC(SPEC `SUC n` DIFF_POW) THEN REWRITE_TAC[SUC_SUB1] THEN
 140     DISCH_THEN(MATCH_ACCEPT_TAC o ONCE_REWRITE_RULE[REAL_MUL_SYM]);
 141     DISCH_THEN(MP_TAC o MATCH_MP DIFF_ADD) THEN BETA_TAC THEN
 142     REWRITE_TAC[REAL_MUL_ASSOC]]);;
 143
 144 (* }}} *)
 145
 146 let POLY_DIFF = prove
 147  (`!l x. ((\x. poly l x) diffl (poly (poly_diff l) x))(x)`,
 148 (* {{{ Proof *)
 149
 150   LIST_INDUCT_TAC THEN REWRITE_TAC[POLY_DIFF_CLAUSES] THEN
 151   ONCE_REWRITE_TAC[SYM(ETA_CONV `\x. poly l x`)] THEN
 152   REWRITE_TAC[poly; DIFF_CONST] THEN
 153   MAP_EVERY X_GEN_TAC [`x:real`] THEN
 154   MP_TAC(SPECL [`t:(real)list`; `0`; `x:real`] POLY_DIFF_LEMMA) THEN
 155   REWRITE_TAC[SYM(num_CONV `1`)] THEN REWRITE_TAC[pow; REAL_MUL_LID] THEN
 156   REWRITE_TAC[POW_1] THEN
 157   DISCH_THEN(MP_TAC o CONJ (SPECL [`h:real`; `x:real`] DIFF_CONST)) THEN
 158   DISCH_THEN(MP_TAC o MATCH_MP DIFF_ADD) THEN BETA_TAC THEN
 159   REWRITE_TAC[REAL_ADD_LID]);;
 160
 161 (* }}} *)
 162
 163 let CANON_POLY_DIFF_CONV =
 164   let aux_conv0 = GEN_REWRITE_CONV I [CONJUNCT1 poly_diff_aux]
 165   and aux_conv1 = GEN_REWRITE_CONV I [CONJUNCT2 poly_diff_aux]
 166   and diff_conv0 = GEN_REWRITE_CONV I (butlast (CONJUNCTS POLY_DIFF_CLAUSES))
 167   and diff_conv1 = GEN_REWRITE_CONV I [last (CONJUNCTS POLY_DIFF_CLAUSES)] in
 168   let rec POLY_DIFF_AUX_CONV tm =
 169    (aux_conv0 ORELSEC
 170     (aux_conv1 THENC
 171      RAND_CONV (LAND_CONV NUM_SUC_CONV THENC POLY_DIFF_AUX_CONV))) tm in
 172   diff_conv0 ORELSEC
 173   (diff_conv1 THENC POLY_DIFF_AUX_CONV);;
 174
 175 (*
 176
 177 # POLY_DIFF_CONV (mk_comb(`poly_diff`,k2));;
 178 val it : Hol.thm =
 179   |- poly_diff k2 =
 180      [&1 * (&0 * z * &1 + y * (&0 + y * &3)); &2 * (&0 + y * (&0 + z * &2))]
 181 *)
 182
 183 (* ------------------------------------------------------------------------- *)
 184 (* Whether the first of two items comes earlier in the list.                 *)
 185 (* ------------------------------------------------------------------------- *)
 186
 187 let rec earlier l x y =
 188   match l with
 189     h::t -> if h = y then false
 190               else if h = x then true
 191               else earlier t x y
 192   | [] -> false;;
 193
 194 (* ------------------------------------------------------------------------- *)
 195 (* Add polynomials.                                                          *)
 196 (* ------------------------------------------------------------------------- *)
 197
 198 let POLY_ADD_CONV =
 199   let [cnv_r; cnv_l; cnv_2; cnv_0] = (map REWR_CONV o CONJUNCTS o REAL_ARITH)
 200     `(pol1 + (d + y * q) = (pol1 + d) + y * q) /\
 201      ((c + x * p) + pol2 = (c + pol2) + x * p) /\
 202      ((c + x * p) + (d + x * q) = (c + d) + x * (p + q)) /\
 203      (c + x * &0 = c)`
 204   and dest_add = dest_binop `(+)` in
 205   let rec POLY_ADD_CONV vars tm =
 206     let pol1,pol2 = dest_add tm in
 207     let x = polyvar pol1 and y = polyvar pol2 in
 208     if not(is_var x) & not(is_var y) then REAL_RAT_REDUCE_CONV tm else
 209     if not(is_var y) or earlier vars x y then
 210       (cnv_l THENC LAND_CONV (POLY_ADD_CONV vars)) tm
 211     else if not(is_var x) or earlier vars y x then
 212       (cnv_r THENC LAND_CONV (POLY_ADD_CONV vars)) tm
 213     else
 214       (cnv_2 THENC COMB_CONV(RAND_CONV(POLY_ADD_CONV vars)) THENC
 215        TRY_CONV cnv_0) tm in
 216   POLY_ADD_CONV;;
 217
 218 (*
 219 # POLY_ADD_CONV [`x:real`;`y:real`;`z:real`] (mk_binop `(+)` k0 k0) ;;
 220 val it : Hol.thm =
 221   |- ((&0 + y * (&0 + z * &3)) +
 222       x *
 223       (((&0 + z * &1) + y * (&0 + y * &3)) + x * (&0 + y * (&0 + z * &2)))) +
 224      (&0 + y * (&0 + z * &3)) +
 225      x * (((&0 + z * &1) + y * (&0 + y * &3)) + x * (&0 + y * (&0 + z * &2))) =
 226      (&0 + y * (&0 + z * &6)) +
 227      x * (((&0 + z * &2) + y * (&0 + y * &6)) + x * (&0 + y * (&0 + z * &4)))
 228 *)
 229
 230 (* ------------------------------------------------------------------------- *)
 231 (* Negate polynomials.                                                       *)
 232 (* ------------------------------------------------------------------------- *)
 233
 234 let POLY_NEG_CONV =
 235   let cnv = REWR_CONV(REAL_ARITH `--(c + x * p) = --c + x * --p`) in
 236   let rec POLY_NEG_CONV tm =
 237     if is_ratconst(rand tm) then REAL_RAT_NEG_CONV tm else
 238     (cnv THENC COMB_CONV(RAND_CONV POLY_NEG_CONV)) tm in
 239   POLY_NEG_CONV;;
 240
 241 (* ------------------------------------------------------------------------- *)
 242 (* Subtract polynomials.                                                     *)
 243 (* ------------------------------------------------------------------------- *)
 244
 245 let POLY_SUB_CONV =
 246   let cnv = REWR_CONV real_sub in
 247   fun vars -> cnv THENC RAND_CONV POLY_NEG_CONV THENC POLY_ADD_CONV vars;;
 248
 249 (* ------------------------------------------------------------------------- *)
 250 (* Multiply polynomials.                                                     *)
 251 (* ------------------------------------------------------------------------- *)
 252
 253 let POLY_MUL_CONV =
 254   let [cnv_l1; cnv_r1; cnv_2; cnv_l0; cnv_r0] =
 255     (map REWR_CONV o CONJUNCTS o REAL_ARITH)
 256     `(pol1 * (d + y * q) = (pol1 * d) + y * (pol1 * q)) /\
 257      ((c + x * p) * pol2 = (c * pol2) + x * (p * pol2)) /\
 258      (pol1 * (d + x * q) = pol1 * d + (&0 + x * pol1 * q)) /\
 259      (&0 * pol2 = &0) /\
 260      (pol1 * &0 = &0)`
 261   and dest_mul = dest_binop `( * )`
 262   and zero_tm = `&0` in
 263   let rec POLY_MUL_CONV vars tm =
 264     let pol1,pol2 = dest_mul tm in
 265     if pol1 = zero_tm then cnv_l0 tm
 266     else if pol2 = zero_tm then cnv_r0 tm
 267     else if is_ratconst pol1 & is_ratconst pol2 then REAL_RAT_MUL_CONV tm else
 268     let x = polyvar pol1 and y = polyvar pol2 in
 269     if not(is_var y) or earlier vars x y then
 270       (cnv_r1 THENC COMB_CONV(RAND_CONV(POLY_MUL_CONV vars))) tm
 271     else if not(is_var x) or earlier vars y x then
 272       (cnv_l1 THENC COMB_CONV(RAND_CONV(POLY_MUL_CONV vars))) tm
 273     else
 274       (cnv_2 THENC COMB2_CONV (RAND_CONV(POLY_MUL_CONV vars))
 275                               (funpow 2 RAND_CONV (POLY_MUL_CONV vars)) THENC
 276        POLY_ADD_CONV vars) tm in
 277   POLY_MUL_CONV;;
 278
 279 (*
 280
 281 # POLY_MUL_CONV [`x:real`;`y:real`;`z:real`] (mk_binop `( * )` k0 k0) ;;
 282 val it : Hol.thm =
 283   |- ((&0 + y * (&0 + z * &3)) +
 284       x *
 285       (((&0 + z * &1) + y * (&0 + y * &3)) + x * (&0 + y * (&0 + z * &2)))) *
 286      ((&0 + y * (&0 + z * &3)) +
 287       x *
 288       (((&0 + z * &1) + y * (&0 + y * &3)) + x * (&0 + y * (&0 + z * &2)))) =
 289      (&0 + y * (&0 + y * (&0 + z * (&0 + z * &9)))) +
 290      x *
 291      ((&0 + y * ((&0 + z * (&0 + z * &6)) + y * (&0 + y * (&0 + z * &18)))) +
 292       x *
 293       (((&0 + z * (&0 + z * &1)) +
 294         y * (&0 + y * ((&0 + z * (&6 + z * &12)) + y * (&0 + y * &9)))) +
 295        x *
 296        ((&0 + y * ((&0 + z * (&0 + z * &4)) + y * (&0 + y * (&0 + z * &12)))) +
 297         x * (&0 + y * (&0 + y * (&0 + z * (&0 + z * &4)))))))
 298 *)
 299
 300
 301 (* ------------------------------------------------------------------------- *)
 302 (* Exponentiate polynomials.                                                 *)
 303 (* ------------------------------------------------------------------------- *)
 304
 305 let POLY_POW_CONV =
 306   let [cnv_0; cnv_1] = map REWR_CONV (CONJUNCTS real_pow)
 307   and zero_tm = `0` in
 308   let rec POLY_POW_CONV vars tm =
 309     if rand tm = zero_tm then cnv_0 tm else
 310     (RAND_CONV num_CONV THENC cnv_1 THENC
 311      RAND_CONV (POLY_POW_CONV vars) THENC
 312      POLY_MUL_CONV vars) tm in
 313   POLY_POW_CONV;;
 314
 315 (*
 316 # POLY_POW_CONV [`x:real`;`y:real`;`z:real`] (mk_binop `(pow)` k0 `2`) ;;
 317 val it : Hol.thm =
 318   |- ((&0 + y * (&0 + z * &3)) +
 319       x *
 320       (((&0 + z * &1) + y * (&0 + y * &3)) + x * (&0 + y * (&0 + z * &2)))) pow
 321      2 =
 322      (&0 + y * (&0 + y * (&0 + z * (&0 + z * &9)))) +
 323      x *
 324      ((&0 + y * ((&0 + z * (&0 + z * &6)) + y * (&0 + y * (&0 + z * &18)))) +
 325       x *
 326       (((&0 + z * (&0 + z * &1)) +
 327         y * (&0 + y * ((&0 + z * (&6 + z * &12)) + y * (&0 + y * &9)))) +
 328        x *
 329        ((&0 + y * ((&0 + z * (&0 + z * &4)) + y * (&0 + y * (&0 + z * &12)))) +
 330         x * (&0 + y * (&0 + y * (&0 + z * (&0 + z * &4)))))))
 331 *)
 332
 333 (* ------------------------------------------------------------------------- *)
 334 (* Convert expression to canonical polynomials.                              *)
 335 (* ------------------------------------------------------------------------- *)
 336
 337 let POLYNATE_CONV =
 338   let cnv_var = REWR_CONV(REAL_ARITH `x = &0 + x * &1`)
 339   and cnv_div = REWR_CONV real_div
 340   and neg_tm = `(--)`
 341   and add_tm = `(+)`
 342   and sub_tm = `(-)`
 343   and mul_tm = `( * )`
 344   and pow_tm = `(pow)`
 345   and div_tm = `(/)` in
 346   let rec POLYNATE_CONV vars tm =
 347     if is_var tm then cnv_var tm
 348     else if is_ratconst tm then REFL tm else
 349     let lop,r = dest_comb tm in
 350     if lop = neg_tm
 351     then (RAND_CONV(POLYNATE_CONV vars) THENC POLY_NEG_CONV) tm else
 352     let op,l = dest_comb lop in
 353     if op = pow_tm then
 354       (LAND_CONV(POLYNATE_CONV vars) THENC POLY_POW_CONV vars) tm
 355     else if op = div_tm then
 356       (cnv_div THENC
 357        COMB2_CONV (RAND_CONV(POLYNATE_CONV vars)) REAL_RAT_REDUCE_CONV THENC
 358        POLY_MUL_CONV vars) tm else
 359     let cnv = if op = add_tm then POLY_ADD_CONV
 360               else if op = sub_tm then POLY_SUB_CONV
 361               else if op = mul_tm then POLY_MUL_CONV
 362               else failwith "POLYNATE_CONV: unknown operation" in
 363     (BINOP_CONV (POLYNATE_CONV vars) THENC cnv vars) tm in
 364   POLYNATE_CONV;;
 365
 366 (*
 367 POLYNATE_CONV [`x:real`;`y:real`] `x + y`;;
 368 POLYNATE_CONV [`x:real`;`y:real`] `x * y + &2 * y`;;
 369 *)
 370
 371 (* ------------------------------------------------------------------------- *)
 372 (* Pure term manipulation versions; will optimize eventually.                *)
 373 (* ------------------------------------------------------------------------- *)
 374
 375 let poly_add_ =
 376   let add_tm = `(+)` in
 377   fun vars p1 p2 ->
 378     rand(concl(POLY_ADD_CONV vars (mk_comb(mk_comb(add_tm,p1),p2))));;
 379
 380 let poly_sub_ =
 381   let sub_tm = `(-)` in
 382   fun vars p1 p2 ->
 383     rand(concl(POLY_SUB_CONV vars (mk_comb(mk_comb(sub_tm,p1),p2))));;
 384
 385 let poly_mul_ =
 386   let mul_tm = `( * )` in
 387   fun vars p1 p2 ->
 388     rand(concl(POLY_MUL_CONV vars (mk_comb(mk_comb(mul_tm,p1),p2))));;
 389
 390 let poly_neg_ =
 391   let neg_tm = `(--)` in
 392   fun p -> rand(concl(POLY_NEG_CONV(mk_comb(neg_tm,p))));;
 393
 394 let poly_pow_ =
 395   let pow_tm = `(pow)` in
 396   fun vars p k ->
 397     rand(concl(POLY_POW_CONV vars
 398       (mk_comb(mk_comb(pow_tm,p),mk_small_numeral k))));;
 399
 400 (* ------------------------------------------------------------------------- *)
 401 (* Get the degree of a polynomial.                                           *)
 402 (* ------------------------------------------------------------------------- *)
 403
 404 let rec degree_ vars tm =
 405   if polyvar tm = hd vars then 1 + degree_ vars (funpow 2 rand tm)
 406   else 0;;
 407
 408 (* ------------------------------------------------------------------------- *)
 409 (* Get the list of coefficients.                                             *)
 410 (* ------------------------------------------------------------------------- *)
 411
 412 let rec coefficients vars tm =
 413   if polyvar tm = hd vars then (lhand tm)::coefficients vars (funpow 2 rand tm)
 414   else [tm];;
 415
 416 (* ------------------------------------------------------------------------- *)
 417 (* Get the head constant.                                                    *)
 418 (* ------------------------------------------------------------------------- *)
 419
 420 let head vars p = last(coefficients vars p);;
 421
 422 (* ---------------------------------------------------------------------- *)
 423 (*  Remove the head coefficient                                           *)
 424 (* ---------------------------------------------------------------------- *)
 425
 426 let rec behead vars tm =
 427   try
 428     let c,r = dest_plus tm in
 429     let x,p = dest_mult r in
 430     if not (x = hd vars) then failwith "" else
 431     let p' = behead vars p in
 432       if p' = rzero then c
 433       else mk_plus c (mk_mult x p')
 434   with _ -> rzero;;
 435
 436 (*
 437 behead [`x:real`] `&1 + x * (&1 + x * (&0 + y * &1))`
 438 *)
 439
 440
 441 let BEHEAD =
 442   let lem = ARITH_RULE `a + b * &0 = a` in
 443   fun vars zthm tm ->
 444     let tm' = behead vars tm in
 445       (* note: pure rewrite is ok here, as tm is in canonical form *)
 446     let thm1 = PURE_REWRITE_CONV[zthm] tm in
 447     let thm2 = PURE_REWRITE_CONV[lem] (rhs(concl thm1)) in
 448     let thm3 = TRANS thm1 thm2 in
 449       thm3;;
 450
 451 let BEHEAD3 =
 452   let lem = ARITH_RULE `a + b * &0 = a` in
 453   fun vars zthm tm ->
 454     let tm' = behead vars tm in
 455     (* note slight hack here:
 456        BEHEAD was working fine if
 457        p = a + x * b where a <> b.  But
 458        when they were equal, dropping multiple levels
 459        broke the reconstruction.  Thus, we only do conversion
 460        on the right except when the head variable has been fully eliminated *)
 461     let conv =
 462       let l,r = dest_binop rp tm in
 463       let l1,r1 = dest_binop rm r in
 464         if l1 = hd vars then RAND_CONV(PURE_ONCE_REWRITE_CONV[zthm])
 465         else PURE_ONCE_REWRITE_CONV[zthm] in
 466     let thm1 = conv tm in
 467     let thm2 = PURE_REWRITE_CONV[lem] (rhs(concl thm1)) in
 468     let thm3 = TRANS thm1 thm2 in
 469       thm3;;
 470
 471 let BEHEAD = BEHEAD3;;
 472
 473
 474 (*
 475 let vars = [`z:real`;`x:real`]
 476 let zthm = (ASSUME `&0 + x * &1 = &0`)
 477 let tm = `(&0 + x * &1) + z * (&0 + x * &1)`
 478 behead vars tm
 479 BEHEAD vars zthm tm
 480 BEHEAD2 vars zthm tm
 481 BEHEAD3 vars zthm tm
 482
 483 let tm = `(&0 + x * &1)`
 484 BEHEAD3 vars zthm tm
 485
 486
 487
 488 let vars = [`x:real`]
 489 let tm = `&1 + x * (&1 + x * (&0 + y * &1))`
 490 let zthm = (ASSUME `&0 + y * &1 = &0`)
 491 BEHEAD vars zthm tm
 492 BEHEAD2 vars zthm tm
 493
 494
 495
 496 *)
 497
 498 (* ------------------------------------------------------------------------- *)
 499 (* Test whether a polynomial is a constant w.r.t. the head variable.         *)
 500 (* ------------------------------------------------------------------------- *)
 501
 502 let is_const_poly vars tm = polyvar tm <> hd vars;;
 503
 504 (* ------------------------------------------------------------------------- *)
 505 (* Get the constant multiple of the "maximal" monomial (implicit lex order)  *)
 506 (* ------------------------------------------------------------------------- *)
 507
 508 let rec headconst p =
 509   try rat_of_term p with Failure _ -> headconst(funpow 2 rand p);;
 510
 511 (* ------------------------------------------------------------------------- *)
 512 (* Monicize; return |- const * pol = monic-pol                               *)
 513 (* ------------------------------------------------------------------------- *)
 514
 515 let MONIC_CONV =
 516   let mul_tm = `( * ):real->real->real` in
 517   fun vars p ->
 518     let c = Int 1 // headconst p in
 519     POLY_MUL_CONV vars (mk_comb(mk_comb(mul_tm,term_of_rat c),p));;
 520
 521 (* ------------------------------------------------------------------------- *)
 522 (* Pseudo-division of s by p; head coefficient of p assumed nonzero.         *)
 523 (* Returns |- a^k s = p q + r for some q and r with deg(r) < deg(p).         *)
 524 (* Optimized only for the trivial case of equal head coefficients; no GCDs.  *)
 525 (* ------------------------------------------------------------------------- *)
 526
 527 let PDIVIDE =
 528   let zero_tm = `&0`
 529   and add0_tm = `(+) (&0)`
 530   and add_tm = `(+)`
 531   and mul_tm = `( * )`
 532   and pow_tm = `(pow)`
 533   and one_tm = `&1` in
 534   let mk_varpow vars k =
 535     let mulx_tm = mk_comb(mul_tm,hd vars) in
 536     funpow k (fun t -> mk_comb(add0_tm,mk_comb(mulx_tm,t))) one_tm in
 537   let rec pdivide_aux vars a n p s =
 538     if s = zero_tm then (0,zero_tm,s) else
 539     let b = head vars s and m = degree_ vars s in
 540     if m < n then (0,zero_tm,s) else
 541     let xp = mk_varpow vars (m - n) in
 542     let p' = poly_mul_ vars xp p in
 543     if a = b then
 544       let (k,q,r) = pdivide_aux vars a n p (poly_sub_ vars s p') in
 545       (k,poly_add_ vars q (poly_mul_ vars xp (poly_pow_ vars a k)),r)
 546     else
 547       let (k,q,r) = pdivide_aux vars a n p
 548         (poly_sub_ vars (poly_mul_ vars a s) (poly_mul_ vars b p')) in
 549       let q' = poly_add_ vars q (poly_mul_ vars b
 550                 (poly_mul_ vars (poly_pow_ vars a k) xp)) in
 551       (k+1,q',r) in
 552   fun vars s p ->
 553     let a = head vars p in
 554     let (k,q,r) = pdivide_aux vars a (degree_ vars p) p s in
 555     let th1 = POLY_MUL_CONV vars (mk_comb(mk_comb(mul_tm,q),p)) in
 556     let th2 = AP_THM (AP_TERM add_tm th1) r in
 557     let th3 = CONV_RULE(RAND_CONV(POLY_ADD_CONV vars)) th2 in
 558     let th4 = POLY_POW_CONV vars
 559      (mk_comb(mk_comb(pow_tm,a),mk_small_numeral k)) in
 560     let th5 = AP_THM (AP_TERM mul_tm th4) s in
 561     let th6 = CONV_RULE(RAND_CONV(POLY_MUL_CONV vars)) th5 in
 562     TRANS th6 (GSYM th3);;
 563
 564 (* ------------------------------------------------------------------------- *)
 565 (* Produce sign theorem for rational constant.                               *)
 566 (* ------------------------------------------------------------------------- *)
 567
 568 let SIGN_CONST =
 569   let zero = Int 0
 570   and zero_tm = `&0`
 571   and eq_tm = `(=):real->real->bool`
 572   and gt_tm = `(>):real->real->bool`
 573   and lt_tm = `(<):real->real->bool` in
 574   fun tm ->
 575     let x = rat_of_term tm in
 576     if x =/ zero then
 577       EQT_ELIM(REAL_RAT_EQ_CONV(mk_comb(mk_comb(eq_tm,tm),zero_tm)))
 578     else if x >/ zero then
 579       EQT_ELIM(REAL_RAT_GT_CONV(mk_comb(mk_comb(gt_tm,tm),zero_tm)))
 580     else
 581       EQT_ELIM(REAL_RAT_LT_CONV(mk_comb(mk_comb(lt_tm,tm),zero_tm)));;
 582
 583 (*
 584 SIGN_CONST `-- &5`;;
 585 val it : Hol.thm = |- &5 > &0
 586 *)
 587
 588
 589 (* ------------------------------------------------------------------------- *)
 590 (* Differentiation conversion in main representation.                        *)
 591 (* ------------------------------------------------------------------------- *)
 592
 593 let POLY_DERIV_CONV =
 594   let poly_diff_tm = `poly_diff`
 595   and pth = GEN_REWRITE_RULE I [SWAP_FORALL_THM] POLY_DIFF in
 596   fun vars tm ->
 597     let th1 = POLY_ENLIST_CONV vars tm in
 598     let th2 = SPECL [hd vars; lhand(rand(concl th1))] pth in
 599     CONV_RULE(RATOR_CONV
 600       (COMB2_CONV (RAND_CONV(ABS_CONV(POLY_DELIST_CONV)))
 601                   (LAND_CONV(CANON_POLY_DIFF_CONV THENC
 602                              LIST_CONV (POLY_MUL_CONV vars)) THENC
 603                    POLY_DELIST_CONV))) th2;;
 604
 605 (*
 606 let k0 = (rhs o concl) (POLYNATE_CONV [`x:real`] `x pow 2 * y`);;
 607 let vars = [`x:real`]
 608 let tm = k0
 609 let k1 = concl th2
 610 let k2 = rator k1
 611 let l,r = dest_comb k2
 612
 613 RATOR_CONV
 614 (RAND_CONV(ABS_CONV(POLY_DELIST_CONV))) l
 615 (LAND_CONV(POLY_DIFF_CONV THENC LIST_CONV (CANON_POLY_MUL_CONV vars)) THENC POLY_DELIST_CONV) r
 616 (LAND_CONV(POLY_DIFF_CONV THENC LIST_CONV (CANON_POLY_MUL_CONV vars))) r
 617 (LAND_CONV(POLY_DIFF_CONV)) r
 618
 619
 620 POLY_DERIV_CONV [`x:real`] (rhs(concl((POLYNATE_CONV [`x:real`] `x pow 2 * y`))));;
 621 val it : Hol.thm =
 622   |- ((\x. &0 + x * (&0 + x * (&0 + y * &1))) diffl &0 + x * (&0 + y * &2)) x
 623 *)
 624
 625 (* ---------------------------------------------------------------------- *)
 626 (*  POLYATOM_CONV                                                         *)
 627 (* ---------------------------------------------------------------------- *)
 628
 629 (*
 630   This is the AFN_CONV argument to the lifting function LIFT_QELIM_CONV
 631 *)
 632
 633 let lt_lem = prove_by_refinement(
 634   `!x y. x < y <=> x - y < &0`,
 635 (* {{{ Proof *)
 636 [
 637   REAL_ARITH_TAC;
 638 ]);;
 639 (* }}} *)
 640
 641 let le_lem = prove_by_refinement(
 642   `!x y. x <= y <=> x - y <= &0`,
 643 (* {{{ Proof *)
 644 [
 645   REAL_ARITH_TAC;
 646 ]);;
 647 (* }}} *)
 648
 649 let eq_lem = prove_by_refinement(
 650   `!x y. (x = y) <=> (x - y = &0)`,
 651 (* {{{ Proof *)
 652 [
 653   REAL_ARITH_TAC;
 654 ]);;
 655 (* }}} *)
 656
 657 let POLYATOM_CONV vars tm =
 658   let thm1 = ONCE_REWRITE_CONV[real_gt;real_ge;eq_lem] tm in
 659   let l,r = dest_eq (concl thm1) in
 660   let thm2 = ONCE_REWRITE_CONV[lt_lem;le_lem] r in
 661   let op,l',r' = get_binop (rhs (concl thm2)) in
 662   let thm3a = POLYNATE_CONV vars l' in
 663   let thm3b = AP_TERM op thm3a in
 664   let thm3 = AP_THM thm3b rzero in
 665     end_itlist TRANS [thm1;thm2;thm3];;
 666
 667 (*
 668
 669 let k0 = `x pow 2 + y * x - &5 > x + &10`
 670 let k0 = `x pow 2 + y * x - &5 >= x + &10`
 671 let k0 = `x pow 2 + y * x - &5 < x + &10`
 672 let k0 = `x pow 2 + y * x - &5 <= x + &10`
 673 let k0 = `x pow 2 + y * x - &5 = x + &10`
 674 let tm = k0;;
 675 let vars = [`x:real`;`y:real`]
 676 POLYATOM_CONV vars  k0
 677
 678 let vars = [`e:real`; `k:real`;`f:real`;`a:real`]
 679 prioritize_real()
 680 let tm = `k < e`
 681
 682 let liouville =
 683  `&6 * (w pow 2 + x pow 2 + y pow 2 + z pow 2) pow 2 =
 684    (((w + x) pow 4 + (w + y) pow 4 + (w + z) pow 4 +
 685      (x + y) pow 4 + (x + z) pow 4 + (y + z) pow 4) +
 686     ((w - x) pow 4 + (w - y) pow 4 + (w - z) pow 4 +
 687      (x - y) pow 4 + (x - z) pow 4 + (y - z) pow 4))`
 688
 689 let lvars = [`w:real`;`x:real`;`y:real`; `z:real`]
 690
 691 POLYATOM_CONV lvars liouville
 692
 693 *)
 694
 695 (* ---------------------------------------------------------------------- *)
 696 (*  Factoring                                                             *)
 697 (* ---------------------------------------------------------------------- *)
 698
 699 let weakfactor x pol =
 700   let rec weakfactor k x pol =
 701     try
 702       let ls,rs = dest_plus pol in
 703       if not (ls = rzero) then failwith "" else
 704       let lm,rm = dest_mult rs in
 705       if not (lm = x) then failwith "" else
 706         weakfactor (k + 1) x rm
 707     with Failure _ ->
 708       k,pol in
 709     weakfactor 0 x pol;;
 710
 711 let poly_var x = mk_plus rzero (mk_mult x rone);;
 712 (*
 713    poly_var rx
 714 *)
 715
 716 let POW_PROD_SUM = prove_by_refinement(
 717   `!x n m. (x pow n) * x pow m = x pow (n + m)`,
 718 (* {{{ Proof *)
 719 [
 720   STRIP_TAC THEN STRIP_TAC THEN INDUCT_TAC;
 721   REWRITE_TAC[real_pow];
 722   NUM_SIMP_TAC;
 723   REAL_SIMP_TAC;
 724   REWRITE_TAC[real_pow];
 725   REWRITE_TAC[ARITH_RULE `n + SUC m = SUC (n + m)`];
 726   REWRITE_TAC[real_pow];
 727   POP_ASSUM (SUBST1_TAC o GSYM);
 728   REAL_ARITH_TAC;
 729 ]);;
 730 (* }}} *)
 731
 732 let lem1 = REAL_ARITH `x * x = x pow 2`;;
 733 let lem2 = GSYM (CONJUNCT2 real_pow);;
 734 let lem3 = REAL_ARITH `!x. x = x pow 1`;;
 735
 736 let SIMP_POW_CONV tm =
 737   let thm1 = ((REWRITE_CONV [GSYM REAL_MUL_ASSOC;lem1;lem2;POW_PROD_SUM]) THENC (ARITH_SIMP_CONV[])) tm in
 738   let _,r = dest_eq (concl thm1) in
 739   if can dest_pow r then thm1 else
 740   let thm2 = ISPEC r lem3 in
 741     thm2;;
 742
 743 (*
 744   SIMP_POW_CONV `x * x * x * x * x`
 745   SIMP_POW_CONV `x * x * (x * x) * x`
 746   SIMP_POW_CONV `x * (x * (x * x)) *(x * x)`
 747   SIMP_POW_CONV `x:real`
 748
 749 *)
 750
 751
 752 let WEAKFACTOR_CONV x pol =
 753   let k,pol' = weakfactor x pol in
 754   let thm1 = ((itlist2 (fun x y z -> ((funpow y RAND_CONV) x) THENC z)
 755                (replicate (GEN_REWRITE_CONV I [REAL_ADD_LID]) k)
 756                (0--(k-1)) ALL_CONV) THENC
 757              (PURE_REWRITE_CONV[REAL_MUL_ASSOC])) pol in
 758   let thm2 = (CONV_RULE (RAND_CONV (LAND_CONV SIMP_POW_CONV))) thm1 in
 759     thm2;;
 760
 761
 762
 763 (*
 764   let pol = `&0 + x * (&0 + x * (&0 + y * &1))`
 765   let pol = `&0 + x * (&0 + x * (&0 + x * (&0 + x * (&0 + x * (&0 + x * (&0 + y * &1))))))`
 766   let pol = `&0 + x * (&0 + x * (&0 + x * (&0 + x * (&0 + x * (&1 + x * (&0 + y * &1))))))`
 767   let pol = `&1 + x * (&0 + x * (&0 + y * &1))`
 768   let pol = `&0 + x * (&1 + x * (&0 + y * &1))`
 769   WEAKFACTOR_CONV rx pol
 770   weakfactor rx pol
 771
 772 *)