realarith.ml

   1 (* ========================================================================= *)
   2 (* Framework for universal real decision procedures, and a simple instance.  *)
   3 (*                                                                           *)
   4 (*       John Harrison, University of Cambridge Computer Laboratory          *)
   5 (*                                                                           *)
   6 (*            (c) Copyright, University of Cambridge 1998                    *)
   7 (*              (c) Copyright, John Harrison 1998-2007                       *)
   8 (* ========================================================================= *)
   9
  10 needs "calc_int.ml";;
  11
  12 (* ------------------------------------------------------------------------- *)
  13 (* Some lemmas needed now just to drive the decision procedure.              *)
  14 (* ------------------------------------------------------------------------- *)
  15
  16 let REAL_LTE_TOTAL = prove
  17  (`!x y. x < y \/ y <= x`,
  18   REWRITE_TAC[real_lt] THEN CONV_TAC TAUT);;
  19
  20 let REAL_LET_TOTAL = prove
  21  (`!x y. x <= y \/ y < x`,
  22   REWRITE_TAC[real_lt] THEN CONV_TAC TAUT);;
  23
  24 let REAL_LT_IMP_LE = prove
  25  (`!x y. x < y ==> x <= y`,
  26   MESON_TAC[real_lt; REAL_LE_TOTAL]);;
  27
  28 let REAL_LTE_TRANS = prove
  29  (`!x y z. x < y /\ y <= z ==> x < z`,
  30   MESON_TAC[real_lt; REAL_LE_TRANS]);;
  31
  32 let REAL_LET_TRANS = prove
  33  (`!x y z. x <= y /\ y < z ==> x < z`,
  34   MESON_TAC[real_lt; REAL_LE_TRANS]);;
  35
  36 let REAL_LT_TRANS = prove
  37  (`!x y z. x < y /\ y < z ==> x < z`,
  38   MESON_TAC[REAL_LTE_TRANS; REAL_LT_IMP_LE]);;
  39
  40 let REAL_LE_ADD = prove
  41  (`!x y. &0 <= x /\ &0 <= y ==> &0 <= x + y`,
  42   MESON_TAC[REAL_LE_LADD_IMP; REAL_ADD_RID; REAL_LE_TRANS]);;
  43
  44 let REAL_LTE_ANTISYM = prove
  45  (`!x y. ~(x < y /\ y <= x)`,
  46   MESON_TAC[real_lt]);;
  47
  48 let REAL_SUB_LE = prove
  49  (`!x y. &0 <= (x - y) <=> y <= x`,
  50   REWRITE_TAC[real_sub; GSYM REAL_LE_LNEG; REAL_LE_NEG2]);;
  51
  52 let REAL_NEG_SUB = prove
  53  (`!x y. --(x - y) = y - x`,
  54   REWRITE_TAC[real_sub; REAL_NEG_ADD; REAL_NEG_NEG] THEN
  55   REWRITE_TAC[REAL_ADD_AC]);;
  56
  57 let REAL_LE_LT = prove
  58  (`!x y. x <= y <=> x < y \/ (x = y)`,
  59   REWRITE_TAC[real_lt] THEN MESON_TAC[REAL_LE_ANTISYM; REAL_LE_TOTAL]);;
  60
  61 let REAL_SUB_LT = prove
  62  (`!x y. &0 < (x - y) <=> y < x`,
  63   REWRITE_TAC[real_lt] THEN ONCE_REWRITE_TAC[GSYM REAL_NEG_SUB] THEN
  64   REWRITE_TAC[REAL_LE_LNEG; REAL_ADD_RID; REAL_SUB_LE]);;
  65
  66 let REAL_NOT_LT = prove
  67  (`!x y. ~(x < y) <=> y <= x`,
  68   REWRITE_TAC[real_lt]);;
  69
  70 let REAL_SUB_0 = prove
  71  (`!x y. (x - y = &0) <=> (x = y)`,
  72   REPEAT GEN_TAC THEN REWRITE_TAC[GSYM REAL_LE_ANTISYM] THEN
  73   GEN_REWRITE_TAC (LAND_CONV o LAND_CONV) [GSYM REAL_NOT_LT] THEN
  74   REWRITE_TAC[REAL_SUB_LE; REAL_SUB_LT] THEN REWRITE_TAC[REAL_NOT_LT]);;
  75
  76 let REAL_LT_LE = prove
  77  (`!x y. x < y <=> x <= y /\ ~(x = y)`,
  78   MESON_TAC[real_lt; REAL_LE_TOTAL; REAL_LE_ANTISYM]);;
  79
  80 let REAL_LT_REFL = prove
  81  (`!x. ~(x < x)`,
  82   REWRITE_TAC[real_lt; REAL_LE_REFL]);;
  83
  84 let REAL_LTE_ADD = prove
  85  (`!x y. &0 < x /\ &0 <= y ==> &0 < x + y`,
  86   MESON_TAC[REAL_LE_LADD_IMP; REAL_ADD_RID; REAL_LTE_TRANS]);;
  87
  88 let REAL_LET_ADD = prove
  89  (`!x y. &0 <= x /\ &0 < y ==> &0 < x + y`,
  90   MESON_TAC[REAL_LTE_ADD; REAL_ADD_SYM]);;
  91
  92 let REAL_LT_ADD = prove
  93  (`!x y. &0 < x /\ &0 < y ==> &0 < x + y`,
  94   MESON_TAC[REAL_LT_IMP_LE; REAL_LTE_ADD]);;
  95
  96 let REAL_ENTIRE = prove
  97  (`!x y. (x * y = &0) <=> (x = &0) \/ (y = &0)`,
  98   REPEAT GEN_TAC THEN EQ_TAC THEN STRIP_TAC THEN
  99   ASM_REWRITE_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO] THEN
 100   ASM_CASES_TAC `x = &0` THEN ASM_REWRITE_TAC[] THEN
 101   FIRST_ASSUM(MP_TAC o AP_TERM `(*) (inv x)`) THEN
 102   REWRITE_TAC[REAL_MUL_ASSOC] THEN
 103   FIRST_ASSUM(fun th -> REWRITE_TAC[MATCH_MP REAL_MUL_LINV th]) THEN
 104   REWRITE_TAC[REAL_MUL_LID; REAL_MUL_RZERO]);;
 105
 106 let REAL_LE_NEGTOTAL = prove
 107  (`!x. &0 <= x \/ &0 <= --x`,
 108   REWRITE_TAC[REAL_LE_RNEG; REAL_ADD_LID; REAL_LE_TOTAL]);;
 109
 110 let REAL_LE_SQUARE = prove
 111  (`!x. &0 <= x * x`,
 112   GEN_TAC THEN DISJ_CASES_TAC(SPEC `x:real` REAL_LE_NEGTOTAL) THEN
 113   POP_ASSUM(fun th -> MP_TAC(MATCH_MP REAL_LE_MUL (CONJ th th))) THEN
 114   REWRITE_TAC[REAL_MUL_LNEG; REAL_MUL_RNEG; REAL_NEG_NEG]);;
 115
 116 let REAL_MUL_RID = prove
 117  (`!x. x * &1 = x`,
 118   MESON_TAC[REAL_MUL_LID; REAL_MUL_SYM]);;
 119
 120 let REAL_POW_2 = prove
 121  (`!x. x pow 2 = x * x`,
 122   REWRITE_TAC[num_CONV `2`; num_CONV `1`] THEN
 123   REWRITE_TAC[real_pow; REAL_MUL_RID]);;
 124
 125 let REAL_POLY_CLAUSES = prove
 126  (`(!x y z. x + (y + z) = (x + y) + z) /\
 127    (!x y. x + y = y + x) /\
 128    (!x. &0 + x = x) /\
 129    (!x y z. x * (y * z) = (x * y) * z) /\
 130    (!x y. x * y = y * x) /\
 131    (!x. &1 * x = x) /\
 132    (!x. &0 * x = &0) /\
 133    (!x y z. x * (y + z) = x * y + x * z) /\
 134    (!x. x pow 0 = &1) /\
 135    (!x n. x pow (SUC n) = x * x pow n)`,
 136   REWRITE_TAC[real_pow; REAL_ADD_LDISTRIB; REAL_MUL_LZERO] THEN
 137   REWRITE_TAC[REAL_MUL_ASSOC; REAL_ADD_LID; REAL_MUL_LID] THEN
 138   REWRITE_TAC[REAL_ADD_AC] THEN REWRITE_TAC[REAL_MUL_SYM]);;
 139
 140 let REAL_POLY_NEG_CLAUSES = prove
 141  (`(!x. --x = --(&1) * x) /\
 142    (!x y. x - y = x + --(&1) * y)`,
 143   REWRITE_TAC[REAL_MUL_LNEG; real_sub; REAL_MUL_LID]);;
 144
 145 let REAL_POS = prove
 146  (`!n. &0 <= &n`,
 147   REWRITE_TAC[REAL_OF_NUM_LE; LE_0]);;
 148
 149 (* ------------------------------------------------------------------------- *)
 150 (* Data structure for Positivstellensatz refutations.                        *)
 151 (* ------------------------------------------------------------------------- *)
 152
 153 type positivstellensatz =
 154    Axiom_eq of int
 155  | Axiom_le of int
 156  | Axiom_lt of int
 157  | Rational_eq of num
 158  | Rational_le of num
 159  | Rational_lt of num
 160  | Square of term
 161  | Eqmul of term * positivstellensatz
 162  | Sum of positivstellensatz * positivstellensatz
 163  | Product of positivstellensatz * positivstellensatz;;
 164
 165 (* ------------------------------------------------------------------------- *)
 166 (* Parametrized reals decision procedure.                                    *)
 167 (*                                                                           *)
 168 (* This is a bootstrapping version, and subsequently gets overwritten twice  *)
 169 (* with more specialized versions, once here and finally in "calc_rat.ml".   *)
 170 (* ------------------------------------------------------------------------- *)
 171
 172 let GEN_REAL_ARITH =
 173   let pth = prove
 174    (`(x < y <=> y - x > &0) /\
 175      (x <= y <=> y - x >= &0) /\
 176      (x > y <=> x - y > &0) /\
 177      (x >= y <=> x - y >= &0) /\
 178      ((x = y) <=> (x - y = &0)) /\
 179      (~(x < y) <=> x - y >= &0) /\
 180      (~(x <= y) <=> x - y > &0) /\
 181      (~(x > y) <=> y - x >= &0) /\
 182      (~(x >= y) <=> y - x > &0) /\
 183      (~(x = y) <=> x - y > &0 \/ --(x - y) > &0)`,
 184     REWRITE_TAC[real_gt; real_ge; REAL_SUB_LT; REAL_SUB_LE; REAL_NEG_SUB] THEN
 185     REWRITE_TAC[REAL_SUB_0; real_lt] THEN MESON_TAC[REAL_LE_ANTISYM])
 186   and pth_final = TAUT `(~p ==> F) ==> p`
 187   and pth_add = prove
 188    (`((x = &0) /\ (y = &0) ==> (x + y = &0)) /\
 189      ((x = &0) /\ y >= &0 ==> x + y >= &0) /\
 190      ((x = &0) /\ y > &0 ==> x + y > &0) /\
 191      (x >= &0 /\ (y = &0) ==> x + y >= &0) /\
 192      (x >= &0 /\ y >= &0 ==> x + y >= &0) /\
 193      (x >= &0 /\ y > &0 ==> x + y > &0) /\
 194      (x > &0 /\ (y = &0) ==> x + y > &0) /\
 195      (x > &0 /\ y >= &0 ==> x + y > &0) /\
 196      (x > &0 /\ y > &0 ==> x + y > &0)`,
 197     SIMP_TAC[REAL_ADD_LID; REAL_ADD_RID; real_ge; real_gt] THEN
 198     REWRITE_TAC[REAL_LE_LT] THEN
 199     MESON_TAC[REAL_ADD_LID; REAL_ADD_RID; REAL_LT_ADD])
 200   and pth_mul = prove
 201    (`((x = &0) /\ (y = &0) ==> (x * y = &0)) /\
 202      ((x = &0) /\ y >= &0 ==> (x * y = &0)) /\
 203      ((x = &0) /\ y > &0 ==> (x * y = &0)) /\
 204      (x >= &0 /\ (y = &0) ==> (x * y = &0)) /\
 205      (x >= &0 /\ y >= &0 ==> x * y >= &0) /\
 206      (x >= &0 /\ y > &0 ==> x * y >= &0) /\
 207      (x > &0 /\ (y = &0) ==> (x * y = &0)) /\
 208      (x > &0 /\ y >= &0 ==> x * y >= &0) /\
 209      (x > &0 /\ y > &0 ==> x * y > &0)`,
 210     SIMP_TAC[REAL_MUL_LZERO; REAL_MUL_RZERO; real_ge; real_gt] THEN
 211     SIMP_TAC[REAL_LT_LE; REAL_LE_MUL] THEN MESON_TAC[REAL_ENTIRE])
 212   and pth_emul = prove
 213    (`(y = &0) ==> !x. x * y = &0`,
 214     SIMP_TAC[REAL_MUL_RZERO])
 215   and pth_square = prove
 216    (`!x. x * x >= &0`,
 217     REWRITE_TAC[real_ge; REAL_POW_2; REAL_LE_SQUARE])
 218   and MATCH_MP_RULE th =
 219     let net = itlist
 220      (fun th -> net_of_conv (lhand(concl th)) (PART_MATCH lhand th))
 221      (CONJUNCTS th) empty_net in
 222     fun th -> MP (REWRITES_CONV net (concl th)) th
 223   and x_tm = `x:real` and y_tm = `y:real`
 224   and neg_tm = `(--):real->real`
 225   and gt_tm = `(>):real->real->bool`
 226   and ge_tm = `(>=):real->real->bool`
 227   and eq_tm = `(=):real->real->bool`
 228   and p_tm = `p:bool`
 229   and or_tm = `(\/)`
 230   and false_tm = `F`
 231   and z_tm = `&0 :real`
 232   and xy_lt = `(x:real) < y`
 233   and xy_nlt = `~((x:real) < y)`
 234   and xy_le = `(x:real) <= y`
 235   and xy_nle = `~((x:real) <= y)`
 236   and xy_gt = `(x:real) > y`
 237   and xy_ngt = `~((x:real) > y)`
 238   and xy_ge = `(x:real) >= y`
 239   and xy_nge = `~((x:real) >= y)`
 240   and xy_eq = `x:real = y`
 241   and xy_ne = `~(x:real = y)` in
 242   let is_ge = is_binop ge_tm
 243   and is_gt = is_binop gt_tm
 244   and is_req = is_binop eq_tm in
 245   fun (mk_numeric,
 246        NUMERIC_EQ_CONV,NUMERIC_GE_CONV,NUMERIC_GT_CONV,
 247        POLY_CONV,POLY_NEG_CONV,POLY_ADD_CONV,POLY_MUL_CONV,
 248        absconv1,absconv2,prover) ->
 249   let REAL_INEQ_CONV pth tm =
 250     let lop,r = dest_comb tm in
 251     let th = INST [rand lop,x_tm; r,y_tm] pth in
 252     TRANS th (LAND_CONV POLY_CONV (rand(concl th))) in
 253   let [REAL_LT_CONV; REAL_LE_CONV; REAL_GT_CONV; REAL_GE_CONV; REAL_EQ_CONV;
 254        REAL_NOT_LT_CONV; REAL_NOT_LE_CONV; REAL_NOT_GT_CONV;
 255        REAL_NOT_GE_CONV; _] =
 256        map REAL_INEQ_CONV (CONJUNCTS pth)
 257   and REAL_NOT_EQ_CONV =
 258      let pth = last(CONJUNCTS pth) in
 259      fun tm ->
 260       let l,r = dest_eq tm in
 261       let th = INST [l,x_tm; r,y_tm] pth in
 262       let th_p = POLY_CONV(lhand(lhand(rand(concl th)))) in
 263       let th_x = AP_TERM neg_tm th_p in
 264       let th_n = CONV_RULE (RAND_CONV POLY_NEG_CONV) th_x in
 265       let th' = MK_DISJ (AP_THM (AP_TERM gt_tm th_p) z_tm)
 266                         (AP_THM (AP_TERM gt_tm th_n) z_tm) in
 267       TRANS th th' in
 268   let net_single = itlist (uncurry net_of_conv)
 269    [xy_lt,REAL_LT_CONV;
 270     xy_nlt,(fun t -> REAL_NOT_LT_CONV(rand t));
 271     xy_le,REAL_LE_CONV;
 272     xy_nle,(fun t -> REAL_NOT_LE_CONV(rand t));
 273     xy_gt,REAL_GT_CONV;
 274     xy_ngt,(fun t -> REAL_NOT_GT_CONV(rand t));
 275     xy_ge,REAL_GE_CONV;
 276     xy_nge,(fun t -> REAL_NOT_GE_CONV(rand t));
 277     xy_eq,REAL_EQ_CONV;
 278     xy_ne,(fun t -> REAL_NOT_EQ_CONV(rand t))]
 279    empty_net
 280   and net_double = itlist (uncurry net_of_conv)
 281    [xy_lt,(fun t -> REAL_LT_CONV t,REAL_NOT_LT_CONV t);
 282     xy_le,(fun t -> REAL_LE_CONV t,REAL_NOT_LE_CONV t);
 283     xy_gt,(fun t -> REAL_GT_CONV t,REAL_NOT_GT_CONV t);
 284     xy_ge,(fun t -> REAL_GE_CONV t,REAL_NOT_GE_CONV t);
 285     xy_eq,(fun t -> REAL_EQ_CONV t,REAL_NOT_EQ_CONV t)]
 286    empty_net in
 287   let REAL_INEQ_NORM_CONV = REWRITES_CONV net_single
 288   and REAL_INEQ_NORM_DCONV = REWRITES_CONV net_double in
 289   let NNF_NORM_CONV =
 290     GEN_NNF_CONV false (REAL_INEQ_NORM_CONV,REAL_INEQ_NORM_DCONV) in
 291   let MUL_RULE =
 292     let rules = MATCH_MP_RULE pth_mul in
 293     fun th -> CONV_RULE(LAND_CONV POLY_MUL_CONV) (rules th)
 294   and ADD_RULE =
 295     let rules = MATCH_MP_RULE pth_add in
 296     fun th -> CONV_RULE(LAND_CONV POLY_ADD_CONV) (rules th)
 297   and EMUL_RULE =
 298     let rule = MATCH_MP pth_emul in
 299     fun tm th -> CONV_RULE (LAND_CONV POLY_MUL_CONV)
 300                            (SPEC tm (rule th))
 301   and SQUARE_RULE t =
 302     CONV_RULE (LAND_CONV POLY_MUL_CONV) (SPEC t pth_square) in
 303   let hol_of_positivstellensatz(eqs,les,lts) =
 304     let rec translate prf =
 305       match prf with
 306         Axiom_eq n -> el n eqs
 307       | Axiom_le n -> el n les
 308       | Axiom_lt n -> el n lts
 309       | Rational_eq x ->
 310           EQT_ELIM(NUMERIC_EQ_CONV(mk_comb(mk_comb(eq_tm,mk_numeric x),z_tm)))
 311       | Rational_le x ->
 312           EQT_ELIM(NUMERIC_GE_CONV(mk_comb(mk_comb(ge_tm,mk_numeric x),z_tm)))
 313       | Rational_lt x ->
 314           EQT_ELIM(NUMERIC_GT_CONV(mk_comb(mk_comb(gt_tm,mk_numeric x),z_tm)))
 315       | Square t -> SQUARE_RULE t
 316       | Eqmul(t,p) -> EMUL_RULE t (translate p)
 317       | Sum(p1,p2) -> ADD_RULE (CONJ (translate p1) (translate p2))
 318       | Product(p1,p2) -> MUL_RULE (CONJ (translate p1) (translate p2)) in
 319     fun prf ->
 320       CONV_RULE(FIRST_CONV[NUMERIC_GE_CONV; NUMERIC_GT_CONV; NUMERIC_EQ_CONV])
 321                (translate prf) in
 322   let init_conv =
 323     TOP_DEPTH_CONV BETA_CONV THENC
 324     PRESIMP_CONV THENC
 325     NNF_CONV THENC DEPTH_BINOP_CONV or_tm CONDS_ELIM_CONV THENC
 326     NNF_NORM_CONV THENC
 327     SKOLEM_CONV THENC
 328     PRENEX_CONV THENC
 329     WEAK_DNF_CONV in
 330   let rec overall dun ths =
 331     match ths with
 332       [] ->
 333         let eq,ne = partition (is_req o concl) dun in
 334         let le,nl = partition (is_ge o concl) ne in
 335         let lt = filter (is_gt o concl) nl in
 336         prover hol_of_positivstellensatz (eq,le,lt)
 337     | th::oths ->
 338         let tm = concl th in
 339         if is_conj tm then
 340           let th1,th2 = CONJ_PAIR th in
 341           overall dun (th1::th2::oths)
 342         else if is_disj tm then
 343           let th1 = overall dun (ASSUME (lhand tm)::oths)
 344           and th2 = overall dun (ASSUME (rand tm)::oths) in
 345           DISJ_CASES th th1 th2
 346         else overall (th::dun) oths in
 347   fun tm ->
 348     let NNF_NORM_CONV' =
 349       GEN_NNF_CONV false
 350         (CACHE_CONV REAL_INEQ_NORM_CONV,fun t -> failwith "") in
 351     let rec absremover t =
 352      (TOP_DEPTH_CONV(absconv1 THENC BINOP_CONV (LAND_CONV POLY_CONV)) THENC
 353       TRY_CONV(absconv2 THENC NNF_NORM_CONV' THENC BINOP_CONV absremover)) t in
 354     let th0 = init_conv(mk_neg tm) in
 355     let tm0 = rand(concl th0) in
 356     let th =
 357       if tm0 = false_tm then fst(EQ_IMP_RULE th0) else
 358       let evs,bod = strip_exists tm0 in
 359       let avs,ibod = strip_forall bod in
 360       let th1 = itlist MK_FORALL avs (DEPTH_BINOP_CONV or_tm absremover ibod) in
 361       let th2 = overall [] [SPECL avs (ASSUME(rand(concl th1)))] in
 362       let th3 =
 363         itlist SIMPLE_CHOOSE evs (PROVE_HYP (EQ_MP th1 (ASSUME bod)) th2) in
 364       DISCH_ALL(PROVE_HYP (EQ_MP th0 (ASSUME (mk_neg tm))) th3) in
 365     MP (INST [tm,p_tm] pth_final) th;;
 366
 367 (* ------------------------------------------------------------------------- *)
 368 (* Linear prover. This works over the rationals in general, but is designed  *)
 369 (* to be OK on integers provided the input contains only integers.           *)
 370 (* ------------------------------------------------------------------------- *)
 371
 372 let REAL_LINEAR_PROVER =
 373   let linear_add = combine (+/) (fun z -> z =/ num_0)
 374   and linear_cmul c = mapf (fun x -> c */ x)
 375   and one_tm = `&1` in
 376   let contradictory p (e,_) =
 377     (is_undefined e & not(p num_0)) or
 378     (dom e = [one_tm] & not(p(apply e one_tm))) in
 379   let rec linear_ineqs vars (les,lts) =
 380     try find (contradictory (fun x -> x >/ num_0)) lts with Failure _ ->
 381     try find (contradictory (fun x -> x >=/ num_0)) les with Failure _ ->
 382     if vars = [] then failwith "linear_ineqs: no contradiction" else
 383     let ineqs = les @ lts in
 384     let blowup v =
 385       length(filter (fun (e,_) -> tryapplyd e v num_0 >/ num_0) ineqs) *
 386       length(filter (fun (e,_) -> tryapplyd e v num_0 </ num_0) ineqs) in
 387     let v =
 388      fst(hd(sort (fun (_,i) (_,j) -> i < j)
 389                  (map (fun v -> v,blowup v) vars))) in
 390     let addup (e1,p1) (e2,p2) acc =
 391       let c1 = tryapplyd e1 v num_0 and c2 = tryapplyd e2 v num_0 in
 392       if c1 */ c2 >=/ num_0 then acc else
 393       let e1' = linear_cmul (abs_num c2) e1
 394       and e2' = linear_cmul (abs_num c1) e2
 395       and p1' = Product(Rational_lt(abs_num c2),p1)
 396       and p2' = Product(Rational_lt(abs_num c1),p2) in
 397       (linear_add e1' e2',Sum(p1',p2'))::acc in
 398     let les0,les1 = partition (fun (e,_) -> tryapplyd e v num_0 =/ num_0) les
 399     and lts0,lts1 = partition (fun (e,_) -> tryapplyd e v num_0 =/ num_0) lts in
 400     let lesp,lesn = partition (fun (e,_) -> tryapplyd e v num_0 >/ num_0) les1
 401     and ltsp,ltsn = partition
 402      (fun (e,_) -> tryapplyd e v num_0 >/ num_0) lts1 in
 403     let les' = itlist (fun ep1 -> itlist (addup ep1) lesp) lesn les0
 404     and lts' = itlist (fun ep1 -> itlist (addup ep1) (lesp@ltsp)) ltsn
 405                       (itlist (fun ep1 -> itlist (addup ep1) (lesn@ltsn)) ltsp
 406                               lts0) in
 407     linear_ineqs (subtract vars [v]) (les',lts') in
 408   let rec linear_eqs(eqs,les,lts) =
 409     try find (contradictory (fun x -> x =/ num_0)) eqs with Failure _ ->
 410     match eqs with
 411       [] -> let vars = subtract
 412              (itlist (union o dom o fst) (les@lts) []) [one_tm] in
 413             linear_ineqs vars (les,lts)
 414     | (e,p)::es -> if is_undefined e then linear_eqs(es,les,lts) else
 415                    let x,c = choose (undefine one_tm e) in
 416                    let xform(t,q as inp) =
 417                      let d = tryapplyd t x num_0 in
 418                      if d =/ num_0 then inp else
 419                      let k = minus_num d */ abs_num c // c in
 420                      let e' = linear_cmul k e
 421                      and t' = linear_cmul (abs_num c) t
 422                      and p' = Eqmul(term_of_rat k,p)
 423                      and q' = Product(Rational_lt(abs_num c),q) in
 424                      linear_add e' t',Sum(p',q') in
 425                    linear_eqs(map xform es,map xform les,map xform lts) in
 426   let linear_prover =
 427    fun (eq,le,lt) ->
 428     let eqs = map2 (fun p n -> p,Axiom_eq n) eq (0--(length eq-1))
 429     and les = map2 (fun p n -> p,Axiom_le n) le (0--(length le-1))
 430     and lts = map2 (fun p n -> p,Axiom_lt n) lt (0--(length lt-1)) in
 431     linear_eqs(eqs,les,lts) in
 432   let lin_of_hol =
 433     let one_tm = `&1`
 434     and zero_tm = `&0`
 435     and add_tm = `(+):real->real->real`
 436     and mul_tm = `(*):real->real->real` in
 437     let rec lin_of_hol tm =
 438       if tm = zero_tm then undefined
 439       else if not (is_comb tm) then (tm |=> Int 1)
 440       else if is_ratconst tm then (one_tm |=> rat_of_term tm) else
 441       let lop,r = dest_comb tm in
 442       if not (is_comb lop) then (tm |=> Int 1) else
 443       let op,l = dest_comb lop in
 444       if op = add_tm then linear_add (lin_of_hol l) (lin_of_hol r)
 445       else if op = mul_tm & is_ratconst l then (r |=> rat_of_term l)
 446       else (tm |=> Int 1) in
 447     lin_of_hol in
 448   let is_alien tm =
 449     match tm with
 450       Comb(Const("real_of_num",_),n) when not(is_numeral n) -> true
 451     | _ -> false in
 452   let n_tm = `n:num` in
 453   let pth = REWRITE_RULE[GSYM real_ge] (SPEC n_tm REAL_POS) in
 454   fun translator (eq,le,lt) ->
 455     let eq_pols = map (lin_of_hol o lhand o concl) eq
 456     and le_pols = map (lin_of_hol o lhand o concl) le
 457     and lt_pols = map (lin_of_hol o lhand o concl) lt in
 458     let aliens =  filter is_alien
 459       (itlist (union o dom) (eq_pols @ le_pols @ lt_pols) []) in
 460     let le_pols' = le_pols @ map (fun v -> (v |=> Int 1)) aliens in
 461     let _,proof = linear_prover(eq_pols,le_pols',lt_pols) in
 462     let le' = le @ map (fun a -> INST [rand a,n_tm] pth) aliens in
 463     translator (eq,le',lt) proof;;
 464
 465 (* ------------------------------------------------------------------------- *)
 466 (* Bootstrapping REAL_ARITH: trivial abs-elim and only integer constants.    *)
 467 (* ------------------------------------------------------------------------- *)
 468
 469 let REAL_ARITH =
 470   let REAL_POLY_NEG_CONV,REAL_POLY_ADD_CONV,REAL_POLY_SUB_CONV,
 471     REAL_POLY_MUL_CONV,REAL_POLY_POW_CONV,REAL_POLY_CONV =
 472   SEMIRING_NORMALIZERS_CONV REAL_POLY_CLAUSES REAL_POLY_NEG_CLAUSES
 473    (is_realintconst,
 474     REAL_INT_ADD_CONV,REAL_INT_MUL_CONV,REAL_INT_POW_CONV)
 475    (<) in
 476   let rule =
 477    GEN_REAL_ARITH
 478    (mk_realintconst,
 479     REAL_INT_EQ_CONV,REAL_INT_GE_CONV,REAL_INT_GT_CONV,
 480     REAL_POLY_CONV,REAL_POLY_NEG_CONV,REAL_POLY_ADD_CONV,REAL_POLY_MUL_CONV,
 481     NO_CONV,NO_CONV,REAL_LINEAR_PROVER)
 482   and deabs_conv = REWRITE_CONV[real_abs; real_max; real_min] in
 483   fun tm ->
 484     let th1 = deabs_conv tm in
 485     EQ_MP (SYM th1) (rule(rand(concl th1)));;
 486
 487 (* ------------------------------------------------------------------------- *)
 488 (* Slightly less parametrized GEN_REAL_ARITH with more intelligent           *)
 489 (* elimination of abs, max and min hardwired in.                             *)
 490 (* ------------------------------------------------------------------------- *)
 491
 492 let GEN_REAL_ARITH =
 493   let ABSMAXMIN_ELIM_CONV1 =
 494     GEN_REWRITE_CONV I [time REAL_ARITH
 495      `(--(&1) * abs(x) >= r <=>
 496        --(&1) * x >= r /\ &1 * x >= r) /\
 497       (--(&1) * abs(x) + a >= r <=>
 498        a + --(&1) * x >= r /\ a + &1 * x >= r) /\
 499       (a + --(&1) * abs(x) >= r <=>
 500        a + --(&1) * x >= r /\ a + &1 * x >= r) /\
 501       (a + --(&1) * abs(x) + b >= r <=>
 502        a + --(&1) * x + b >= r /\ a + &1 * x + b >= r) /\
 503       (a + b + --(&1) * abs(x) >= r <=>
 504        a + b + --(&1) * x >= r /\ a + b + &1 * x >= r) /\
 505       (a + b + --(&1) * abs(x) + c >= r <=>
 506        a + b + --(&1) * x + c >= r /\ a + b + &1 * x + c >= r) /\
 507       (--(&1) * max x y >= r <=>
 508        --(&1) * x >= r /\ --(&1) * y >= r) /\
 509       (--(&1) * max x y + a >= r <=>
 510        a + --(&1) * x >= r /\ a + --(&1) * y >= r) /\
 511       (a + --(&1) * max x y >= r <=>
 512        a + --(&1) * x >= r /\ a + --(&1) * y >= r) /\
 513       (a + --(&1) * max x y + b >= r <=>
 514        a + --(&1) * x + b >= r /\ a + --(&1) * y  + b >= r) /\
 515       (a + b + --(&1) * max x y >= r <=>
 516        a + b + --(&1) * x >= r /\ a + b + --(&1) * y >= r) /\
 517       (a + b + --(&1) * max x y + c >= r <=>
 518        a + b + --(&1) * x + c >= r /\ a + b + --(&1) * y  + c >= r) /\
 519       (&1 * min x y >= r <=>
 520        &1 * x >= r /\ &1 * y >= r) /\
 521       (&1 * min x y + a >= r <=>
 522        a + &1 * x >= r /\ a + &1 * y >= r) /\
 523       (a + &1 * min x y >= r <=>
 524        a + &1 * x >= r /\ a + &1 * y >= r) /\
 525       (a + &1 * min x y + b >= r <=>
 526        a + &1 * x + b >= r /\ a + &1 * y  + b >= r) /\
 527       (a + b + &1 * min x y >= r <=>
 528        a + b + &1 * x >= r /\ a + b + &1 * y >= r) /\
 529       (a + b + &1 * min x y + c >= r <=>
 530        a + b + &1 * x + c >= r /\ a + b + &1 * y  + c >= r) /\
 531       (min x y >= r <=>
 532         x >= r /\  y >= r) /\
 533       (min x y + a >= r <=>
 534        a + x >= r /\ a + y >= r) /\
 535       (a + min x y >= r <=>
 536        a + x >= r /\ a + y >= r) /\
 537       (a + min x y + b >= r <=>
 538        a + x + b >= r /\ a + y  + b >= r) /\
 539       (a + b + min x y >= r <=>
 540        a + b + x >= r /\ a + b + y >= r) /\
 541       (a + b + min x y + c >= r <=>
 542        a + b + x + c >= r /\ a + b + y + c >= r) /\
 543       (--(&1) * abs(x) > r <=>
 544        --(&1) * x > r /\ &1 * x > r) /\
 545       (--(&1) * abs(x) + a > r <=>
 546        a + --(&1) * x > r /\ a + &1 * x > r) /\
 547       (a + --(&1) * abs(x) > r <=>
 548        a + --(&1) * x > r /\ a + &1 * x > r) /\
 549       (a + --(&1) * abs(x) + b > r <=>
 550        a + --(&1) * x + b > r /\ a + &1 * x + b > r) /\
 551       (a + b + --(&1) * abs(x) > r <=>
 552        a + b + --(&1) * x > r /\ a + b + &1 * x > r) /\
 553       (a + b + --(&1) * abs(x) + c > r <=>
 554        a + b + --(&1) * x + c > r /\ a + b + &1 * x + c > r) /\
 555       (--(&1) * max x y > r <=>
 556        --(&1) * x > r /\ --(&1) * y > r) /\
 557       (--(&1) * max x y + a > r <=>
 558        a + --(&1) * x > r /\ a + --(&1) * y > r) /\
 559       (a + --(&1) * max x y > r <=>
 560        a + --(&1) * x > r /\ a + --(&1) * y > r) /\
 561       (a + --(&1) * max x y + b > r <=>
 562        a + --(&1) * x + b > r /\ a + --(&1) * y  + b > r) /\
 563       (a + b + --(&1) * max x y > r <=>
 564        a + b + --(&1) * x > r /\ a + b + --(&1) * y > r) /\
 565       (a + b + --(&1) * max x y + c > r <=>
 566        a + b + --(&1) * x + c > r /\ a + b + --(&1) * y  + c > r) /\
 567       (min x y > r <=>
 568         x > r /\  y > r) /\
 569       (min x y + a > r <=>
 570        a + x > r /\ a + y > r) /\
 571       (a + min x y > r <=>
 572        a + x > r /\ a + y > r) /\
 573       (a + min x y + b > r <=>
 574        a + x + b > r /\ a + y  + b > r) /\
 575       (a + b + min x y > r <=>
 576        a + b + x > r /\ a + b + y > r) /\
 577       (a + b + min x y + c > r <=>
 578        a + b + x + c > r /\ a + b + y + c > r)`]
 579   and ABSMAXMIN_ELIM_CONV2 =
 580     let pth_abs = prove
 581      (`P(abs x) <=> (x >= &0 /\ P x) \/ (&0 > x /\ P (--x))`,
 582       REWRITE_TAC[real_abs; real_gt; real_ge] THEN COND_CASES_TAC THEN
 583       ASM_REWRITE_TAC[real_lt])
 584     and pth_max = prove
 585      (`P(max x y) <=> (y >= x /\ P y) \/ (x > y /\ P x)`,
 586       REWRITE_TAC[real_max; real_gt; real_ge] THEN
 587       COND_CASES_TAC THEN ASM_REWRITE_TAC[real_lt])
 588     and pth_min = prove
 589     (`P(min x y) <=> (y >= x /\ P x) \/ (x > y /\ P y)`,
 590       REWRITE_TAC[real_min; real_gt; real_ge] THEN
 591       COND_CASES_TAC THEN ASM_REWRITE_TAC[real_lt])
 592     and abs_tm = `real_abs`
 593     and p_tm = `P:real->bool`
 594     and x_tm = `x:real`
 595     and y_tm = `y:real` in
 596     let is_max = is_binop `real_max`
 597     and is_min = is_binop `real_min`
 598     and is_abs t = is_comb t & rator t = abs_tm in
 599     let eliminate_construct p c tm =
 600       let t = find_term (fun t -> p t & free_in t tm) tm in
 601       let v = genvar(type_of t) in
 602       let th0 = SYM(BETA_CONV(mk_comb(mk_abs(v,subst[v,t] tm),t))) in
 603       let p,ax = dest_comb(rand(concl th0)) in
 604       CONV_RULE(RAND_CONV(BINOP_CONV(RAND_CONV BETA_CONV)))
 605                (TRANS th0 (c p ax)) in
 606     let elim_abs =
 607       eliminate_construct is_abs
 608         (fun p ax -> INST [p,p_tm; rand ax,x_tm] pth_abs)
 609     and elim_max =
 610       eliminate_construct is_max
 611         (fun p ax -> let ax,y = dest_comb ax in
 612                      INST [p,p_tm; rand ax,x_tm; y,y_tm] pth_max)
 613     and elim_min =
 614       eliminate_construct is_min
 615         (fun p ax -> let ax,y = dest_comb ax in
 616                      INST [p,p_tm; rand ax,x_tm; y,y_tm] pth_min) in
 617     FIRST_CONV [elim_abs; elim_max; elim_min] in
 618   fun (mkconst,EQ,GE,GT,NORM,NEG,ADD,MUL,PROVER) ->
 619         GEN_REAL_ARITH(mkconst,EQ,GE,GT,NORM,NEG,ADD,MUL,
 620                        ABSMAXMIN_ELIM_CONV1,ABSMAXMIN_ELIM_CONV2,PROVER);;
 621
 622 (* ------------------------------------------------------------------------- *)
 623 (* Incorporate that. This gets overwritten again in "calc_rat.ml".           *)
 624 (* ------------------------------------------------------------------------- *)
 625
 626 let REAL_ARITH =
 627   let REAL_POLY_NEG_CONV,REAL_POLY_ADD_CONV,REAL_POLY_SUB_CONV,
 628     REAL_POLY_MUL_CONV,REAL_POLY_POW_CONV,REAL_POLY_CONV =
 629   SEMIRING_NORMALIZERS_CONV REAL_POLY_CLAUSES REAL_POLY_NEG_CLAUSES
 630    (is_realintconst,
 631     REAL_INT_ADD_CONV,REAL_INT_MUL_CONV,REAL_INT_POW_CONV)
 632    (<) in
 633   GEN_REAL_ARITH
 634    (mk_realintconst,
 635     REAL_INT_EQ_CONV,REAL_INT_GE_CONV,REAL_INT_GT_CONV,
 636     REAL_POLY_CONV,REAL_POLY_NEG_CONV,REAL_POLY_ADD_CONV,REAL_POLY_MUL_CONV,
 637     REAL_LINEAR_PROVER);;