diff --git a/approxCompErrScript.sml b/approxCompErrScript.sml index 0e2e095f0902a5d1216bab22d77f90e0f9ee4758..c2292de586962854b63e5e7efd073ef28bbde818 100644 --- a/approxCompErrScript.sml +++ b/approxCompErrScript.sml @@ -1,7 +1,7 @@ open IntervalArithTheory ErrorValidationTheory sqrtApproxTheory; open mcLaurinApproxTheory realPolyTheory realPolyProofsTheory transcLangTheory approxPolyTheory transcIntvSemTheory; -open preamble; +open preambleDandelion; val _ = new_theory "approxCompErr"; diff --git a/approxPolyScript.sml b/approxPolyScript.sml index 770e80eb8efe558099249ac2061817ad8f9c9426..f0149e76971aa913af78d9bc253bdb586e91f091 100644 --- a/approxPolyScript.sml +++ b/approxPolyScript.sml @@ -1,7 +1,7 @@ (** **) open realTheory realLib RealArith transcTheory; open realPolyTheory realPolyProofsTheory mcLaurinApproxTheory transcLangTheory; -open preamble; +open preambleDandelion; val _ = new_theory "approxPoly"; @@ -66,6 +66,8 @@ Definition approxPoly_def: else NONE | Sin => SOME (sinPoly approxSteps, sinErr iv approxSteps) | Cos => SOME (cosPoly approxSteps, cosErr iv approxSteps) + | Tan => NONE + | Sqrt => NONE End (** Simple properties of polynomials used for proofs later **) @@ -486,10 +488,7 @@ Proof >> qexists_tac ‘abs (x pow (2 * m'))’ >> gs[ABS_LE, GSYM POW_ABS] >> irule POW_LE >> gs[ABS_POS] >> irule RealSimpsTheory.maxAbs >> gs[]) - (* Tan function *) - >- ( cheat ) - (* Sqrt function *) - >> cheat + (* Tan and Sqrt function missing here as they are unimplemented *) QED val _ = export_theory(); diff --git a/bitArithScript.sml b/bitArithScript.sml index 4123036a6594eeaba3a8e7b171e5ffa79b8c6f15..a64d96ade74d6989acc2489962c31409af741c73 100644 --- a/bitArithScript.sml +++ b/bitArithScript.sml @@ -1,5 +1,5 @@ open HolKernel Parse boolLib bossLib arithmeticTheory; -open preamble; +open preambleDandelion; (** Check reduceLib.num_compset() **) val _ = new_theory "bitArith"; diff --git a/checkerDefsScript.sml b/checkerDefsScript.sml index ecf12bbcc3ca9f9210d230fb3f1c3ddf6ad676a3..07ef42787cbb16ac2f0ee20d39db8d872efe5911 100644 --- a/checkerDefsScript.sml +++ b/checkerDefsScript.sml @@ -1,6 +1,6 @@ open realTheory realLib RealArith stringTheory; open renameTheory realPolyTheory transcLangTheory hintTheory; -open preamble; +open preambleDandelion; val _ = new_theory"checkerDefs"; diff --git a/checkerScript.sml b/checkerScript.sml index 4d7c69b4b5f14be85fb063b12142431ba0b665b7..67ed56954f74dbcea46d97b5c4cc0758100ee7ab 100644 --- a/checkerScript.sml +++ b/checkerScript.sml @@ -6,7 +6,7 @@ open realTheory realLib RealArith stringTheory polyTheory transcTheory; open renameTheory realPolyTheory transcLangTheory sturmComputeTheory sturmTheory drangTheory checkerDefsTheory pointCheckerTheory mcLaurinApproxTheory hintTheory realPolyProofsTheory; -open preamble; +open preambleDandelion; val _ = new_theory "checker"; diff --git a/drangScript.sml b/drangScript.sml index e5a64f2bcf0bdfd1ecb1e3f0423c02a55edfaf00..ee615caea4959684e20d51e5591b53c411af1e4d 100644 --- a/drangScript.sml +++ b/drangScript.sml @@ -1,6 +1,6 @@ open bossLib RealArith realTheory polyTheory limTheory; open renameTheory; -open preamble; +open preambleDandelion; val _ = new_theory "drang"; diff --git a/euclidDivScript.sml b/euclidDivScript.sml index 4d80b8ee1ae7e549dfd08515fca88960a3b05e70..c1d906791da25ef92762094f5cfe6b9e67cde00c 100644 --- a/euclidDivScript.sml +++ b/euclidDivScript.sml @@ -2,7 +2,7 @@ open pred_setTheory listTheory bossLib RealArith realTheory polyTheory; open realPolyTheory sturmTheory realPolyProofsTheory; open renameTheory; -open preamble; +open preambleDandelion; val _ = new_theory "euclidDiv"; diff --git a/floverConnScript.sml b/floverConnScript.sml index c01959a4761c8840c08946341538c2256d9e66d9..92540476c97d642b9a648865abe5b8636928cff4 100644 --- a/floverConnScript.sml +++ b/floverConnScript.sml @@ -3,7 +3,7 @@ **) open ExpressionsTheory ExpressionSemanticsTheory realPolyTheory -open preamble; +open preambleDandelion; val _ = new_theory "floverConn"; diff --git a/hintScript.sml b/hintScript.sml index c4c8192b11d2d73d05e8c65af22dee370c35a163..ad3f5a38551c437d3b742f0b33b8c8fc6567546c 100644 --- a/hintScript.sml +++ b/hintScript.sml @@ -1,5 +1,5 @@ open realTheory realLib RealArith stringTheory polyTheory transcTheory; -open preamble; +open preambleDandelion; val _ = new_theory "hint"; diff --git a/mcLaurinApproxScript.sml b/mcLaurinApproxScript.sml index 9d7b5454ac2471d523229948866bfa6570e42fcd..6ce55b3bdaff5c3e28052ebf06130d8ea2568154 100644 --- a/mcLaurinApproxScript.sml +++ b/mcLaurinApproxScript.sml @@ -2,7 +2,7 @@ in the transcLang file **) open moreRealTheory realPolyTheory realPolyProofsTheory; -open preamble; +open preambleDandelion; val _ = new_theory "mcLaurinApprox"; diff --git a/moreRealScript.sml b/moreRealScript.sml index 27e83dc19fe4c1bda040f0e215e0729c57f0d13a..7a8c25baa7e040622ea80a214a4e81965d2a0725 100644 --- a/moreRealScript.sml +++ b/moreRealScript.sml @@ -1,4 +1,4 @@ -open preamble; +open preambleDandelion; val _ = new_theory "moreReal"; diff --git a/pointCheckerProofsScript.sml b/pointCheckerProofsScript.sml index ca13cf43b006f5fa5753d427ee71b4c68f92215e..4b2e06dca06eaaab12cc4ec2634e44f5d7a357ed 100644 --- a/pointCheckerProofsScript.sml +++ b/pointCheckerProofsScript.sml @@ -4,7 +4,7 @@ open realTheory realLib RealArith stringTheory; open renameTheory realPolyTheory checkerDefsTheory pointCheckerTheory; open realPolyProofsTheory; -open preamble; +open preambleDandelion; val _ = new_theory "pointCheckerProofs"; diff --git a/preamble.sml b/preambleDandelion.sml similarity index 98% rename from preamble.sml rename to preambleDandelion.sml index 5d098b58aaf71fd7f22ad7799c462add312f7ad4..33ca1a3acce837cc0969d83839092b6e1c61bc18 100644 --- a/preamble.sml +++ b/preambleDandelion.sml @@ -2,7 +2,7 @@ Proof tools (e.g. tactics) used throughout the development. Copied from CakeML (https://github.com/CakeML/cakeml) *) -structure preamble = +structure preambleDandelion = struct local open intLib in end; open set_relationTheory; (* comes first so relationTheory takes precedence *) diff --git a/realPolyProofsScript.sml b/realPolyProofsScript.sml index ace1d932a1e75743f06561c7ad854b5604907034..737986e4c4c9c4f2259437898ee37cc0e63b0301 100644 --- a/realPolyProofsScript.sml +++ b/realPolyProofsScript.sml @@ -3,7 +3,7 @@ **) open realTheory realLib RealArith renameTheory polyTheory; open realPolyTheory; -open preamble; +open preambleDandelion; val _ = new_theory "realPolyProofs"; diff --git a/realPolyScript.sml b/realPolyScript.sml index b939000cd63e6fa61bea948f99638c64944f05ab..7adecfebc927cb2ba6e0aa276775e24d922cb05b 100644 --- a/realPolyScript.sml +++ b/realPolyScript.sml @@ -8,7 +8,7 @@ **) open realTheory realLib RealArith bossLib polyTheory; open renameTheory; -open preamble; +open preambleDandelion; val _ = new_theory "realPoly"; diff --git a/realZeroLib.sml b/realZeroLib.sml index e1485fea4780e10e269a2434146d57294eee8439..71f77d214135d9a4577295c90fdcc2a638dfd75a 100644 --- a/realZeroLib.sml +++ b/realZeroLib.sml @@ -4,12 +4,13 @@ struct open realPolyTheory realPolyProofsTheory checkerTheory sturmComputeTheory transcLangTheory transcIntvSemTheory approxPolyTheory transcApproxSemTheory transcReflectTheory euclidDivTheory; - open bossLib preamble; + open bitArithLib bossLib preamble; - exception ZeroLibErr of string; + exception ZeroLibErr of string; val _ = computeLib.add_thms [REAL_INV_1OVER] computeLib.the_compset; val _ = computeLib.add_funs [polyTheory.poly_diff_def, polyTheory.poly_diff_aux_def, polyTheory.poly_def] + val _ = bitArithLib.use_karatsuba(); fun appOrErr (P:term -> bool) (f:term -> 'a) (t:term) (errMsg:string) :'a = if P t then f t else raise ZeroLibErr errMsg; @@ -357,3 +358,172 @@ Q.prove (‘(! x. xlo <= x /\ x <= xhi ==> (P x /\ Q x)) <=> (* val t = REAL_ZERO_CONV “! x. 90 <= x /\ x <= 100 ==> evalPoly [1; 2/100; 3] x <= 100:real†*) end; + +(** UNUSED TESTS + + local + fun rhs_fun_conv c = RIGHT_CONV_RULE c + + fun poly_mul_rec_conv tm = let + val p1 = rand $ rator tm + val p2 = rand tm + in + if listSyntax.is_nil p1 then EVAL tm + else if listSyntax.is_nil p2 then SPEC p1 (GEN_ALL mul_0_right) + else + let + val (hdTm, tlTm) = listSyntax.dest_cons p2 + val tmpTh = + REWR_CONV (hd $ tl (CONJ_LIST 2 poly_mul_def)) tm + |> timed "mul_cst" (rhs_fun_conv (RATOR_CONV $ RAND_CONV EVAL)) + |> timed "cond" (rhs_fun_conv (RAND_CONV $ RATOR_CONV $ RATOR_CONV $ RAND_CONV EVAL)) + |> rhs_fun_conv (RAND_CONV COND_CONV) + val rhsTm = tmpTh |> rhs o concl |> rand + in + if listSyntax.is_nil rhsTm then tmpTh (* |> timed "add" (rhs_fun_conv EVAL) *) + else + tmpTh |> rhs_fun_conv (RAND_CONV $ RAND_CONV poly_mul_rec_conv) + (* |> timed "add" (rhs_fun_conv EVAL) *) + end + end; + + fun poly_add_0_cons_conv tm = + (REWR_CONV (GSYM reduce_add_r) tm + |> timed "push_reduce" rhs_fun_conv (RAND_CONV rec_reduce_push_conv) + |> timed "prep" rhs_fun_conv (REWR_CONV poly_add_assoc) + |> timed "add_single" rhs_fun_conv (RATOR_CONV $ RAND_CONV EVAL) + |> rhs_fun_conv poly_add_0_cons_conv) + (* |> rhs_fun_conv EVAL) *) + handle HOL_ERR _ => + (print ("Error found for \n"^term_to_string tm); + EVAL tm); + + fun poly_mul_conv tm = + poly_mul_rec_conv tm + |> timed "add" (rhs_fun_conv $ poly_add_0_cons_conv) + +Theorem poly_add_cnst: + ! c p1 p2. reduce (c::(p1 +p p2)) = (0::p1 +p c::p2) +Proof + Induct_on ‘p1’ >> rw[poly_add_def, poly_add_aux_def] + >- gs[reduce_def, reduce_idempotent] + >> Cases_on ‘p2’ >> gs[] + >> gs[reduce_def, reduce_idempotent] +QED + +Theorem reduce_reduce: + reduce (c::cs) = reduce (c::reduce cs) +Proof + Induct_on ‘cs’ >> rw[reduce_def] >> gs[reduce_def, reduce_idempotent] +QED + +fun iterconv n c = if n <= 0 then REFL else c THENC (iterconv (n-1) c) + +(** Assumes term of shape reduce (x :: y :: ...) **) +fun rec_reduce_push_conv tm = + REWR_CONV poly_add_cnst tm + handle HOL_ERR _ => + tm |> REWR_CONV reduce_reduce + |> rhs_fun_conv (REWR_CONV reduce_idempotent ORELSEC REFL) + |> rhs_fun_conv $ RAND_CONV $ RAND_CONV rec_reduce_push_conv + |> rhs_fun_conv $ REWR_CONV poly_add_cnst + +val tm = + “[124544987/128396; 192832957/ 1984792857; 345637978/ 18293756379; 348965798567/1090768937; + 124544987/128396; 192832957/ 1984792857; 345637978/ 18293756379; 348965798567/1090768937; + 124544987/128396; 192832957/ 1984792857; 345637978/ 18293756379; 348965798567/1090768937; + 124544987/128396; 192832957/ 1984792857; 345637978/ 18293756379; 348965798567/1090768937; + 124544987/128396; 192832957/ 1984792857; 345637978/ 18293756379; 348965798567/1090768937] *p + [8934769857/1289768935; 890754897/8975641; 128978967389/857923; 77865897298/2089765489; + 8934769857/1289768935; 890754897/8975641; 128978967389/857923; 77865897298/2089765489; + 8934769857/1289768935; 890754897/8975641; 128978967389/857923; 77865897298/2089765489; + 8934769857/1289768935; 890754897/8975641; 128978967389/857923; 77865897298/2089765489; + 8934769857/1289768935; 890754897/8975641; 128978967389/857923; 77865897298/2089765489]†+ +poly_mul_conv tm + val tm = poly_mul_rec_conv tm |> rhs o concl + val tm2 = poly_add_0_cons_conv tm + val tm1 = rhs o concl $ + (REWR_CONV (GSYM reduce_add_r) tm + |> rhs_fun_conv (RAND_CONV $ iterconv 0 $ REWR_CONV reduce_reduce) + |> rhs_fun_conv (RAND_CONV $ REPEATC $ REWRITE_CONV [poly_add_cnst]) + |> timed "prep" rhs_fun_conv (REWR_CONV poly_add_assoc) + |> timed "add_single" rhs_fun_conv (RATOR_CONV $ RAND_CONV EVAL)) + val tm2 = + rhs o concl $ + (REWR_CONV (GSYM reduce_add_r) tm1 + |> rhs_fun_conv (RAND_CONV $ iterconv 1 $ REWR_CONV reduce_reduce) + |> rhs_fun_conv (RAND_CONV $ REPEATC $ REWRITE_CONV [poly_add_cnst]) + |> timed "prep" rhs_fun_conv (REWR_CONV poly_add_assoc) + |> timed "add_single" rhs_fun_conv (RATOR_CONV $ RAND_CONV EVAL)) + val tm3 = + rhs o concl $ + (REWR_CONV (GSYM reduce_add_r) tm2 + |> rhs_fun_conv (RAND_CONV $ iterconv 2 $ REWR_CONV reduce_reduce) + |> rhs_fun_conv (RAND_CONV $ REPEATC $ REWRITE_CONV [poly_add_cnst]) + |> timed "prep" rhs_fun_conv (REWR_CONV poly_add_assoc) + |> timed "add_single" rhs_fun_conv (RATOR_CONV $ RAND_CONV EVAL)) +p + +EVAL tm + +[1/2] *p [1/4]†+ +TOR_CONV $ RATOR_CONV $ RAND_CONV EVAL))) + +raise ZeroLibErr "Foo" end; + + + fun karatsuba_rec_conv tm = let + val bs1 = rand $ rator $ rand tm + val bs2 = rand $ rand tm + val len_bs1 = numSyntax.dest_numeral $ rhs o concl $ EVAL “LENGTH ^bs1†+ val len_bs2 = numSyntax.dest_numeral $ rhs o concl $ EVAL “LENGTH ^bs2†+ in + if Arbnum.<(len_bs1,!karatsuba_lim) orelse Arbnum.<(len_bs2,!karatsuba_lim) + then (RAND_CONV EVAL) tm + else let + val _ = print "." + val th = SPECL [bs1, bs2] karatsuba_bit + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let d = ... *) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let x1 = ... *) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let x0 = ... *) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let y1 = ... *) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let y0 = ... *) + (** Ugly: Get rid of lets for mults **) + val z0 = th |> rhs o concl |> rand |> rand + val th = th |> rhs_fun_conv subst_let_conv (* inline z0 *) + val z2 = th |> rhs o concl |> rand |> rand + val th = th |> rhs_fun_conv subst_let_conv (* inline z2 *) + (** Continue eval **) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let z1a = ... *) + val th = th |> rhs_fun_conv (eval_let_conv EVAL) (* let z1b = ... *) + (** Ugly: More inlining **) + val z1Mul = th |> rhs o concl |> rand |> rand + val th = th |> rhs_fun_conv subst_let_conv (* inline z1Mul *) + val th = th |> rhs_fun_conv subst_let_conv (* inline z1Sub *) + (** Now evaluate the terms we inlined **) + val z0Eval = karatsuba_rec_conv “bleval ^z0†(* z0 *) + val z2Eval = karatsuba_rec_conv “bleval ^z2†(* z2 *) + val z1MulEval = karatsuba_rec_conv “bleval ^z1Mul†(* z1Mul **) + (** Now first get down to ‘bleval (mul _ _)’ terms **) + (** TODO: Prove a single rewriting theorem for this part in the needed shape and apply it with REWR_CONV **) + val th2 = th |> REWRITE_RULE [add_thm, mulpow2_thm, sub_thm] + val th3 = th2 |> REWRITE_RULE [z0Eval, z2Eval, z1MulEval] + (** TODO: Use a single use theorem here too *) + val th4 = th3 |> REWRITE_RULE [GSYM add_thm, GSYM mulpow2_thm, GSYM sub_thm] + in + th4 |> rhs_fun_conv EVAL + end end; + in + fun karatsuba_conv tm = + let val (arg1, arg2) = numSyntax.dest_mult tm + val arg1_bval = ONCE_REWRITE_CONV [GSYM tobl_correct] arg1 |> RIGHT_CONV_RULE $ RAND_CONV EVAL + val arg2_bval = ONCE_REWRITE_CONV [GSYM tobl_correct] arg2 |> RIGHT_CONV_RULE $ RAND_CONV EVAL + val th = REWRITE_CONV [arg1_bval, arg2_bval] tm |> REWRITE_RULE [GSYM mul_thm] + val karat = th |> RIGHT_CONV_RULE $ karatsuba_rec_conv + in + karat |> REWRITE_RULE [GSYM fromBL_correct] |> RIGHT_CONV_RULE EVAL + end + end; +**) diff --git a/renameScript.sml b/renameScript.sml index 3e26fa587886ea4f267d6194cf12a3459514bb69..c60f7ead9b49dec877d09113195f8f714860911b 100644 --- a/renameScript.sml +++ b/renameScript.sml @@ -2,7 +2,7 @@ renaming theory to unify naming of theorems **) -open preamble; +open preambleDandelion; val _ = new_theory "rename"; diff --git a/sturmComputeScript.sml b/sturmComputeScript.sml index 3ad681e7ef85e783a099e95a6d33e83a0b220312..0932c53b421d7a05b169275c6cd3986e17294fe4 100644 --- a/sturmComputeScript.sml +++ b/sturmComputeScript.sml @@ -5,10 +5,138 @@ open pred_setTheory listTheory bossLib RealArith realTheory polyTheory; open realPolyTheory sturmTheory realPolyProofsTheory euclidDivTheory; open renameTheory; -open preamble; +open preambleDandelion; val _ = new_theory "sturmCompute"; +Definition sturm_seq_step_def: + sturm_seq_step p q = + if zerop q ∨ LENGTH (reduce q) < 2 then NONE else + (let g = --p (rm p (inv (coeff q (deg q)) *c q)) in + SOME g) +End + +fun gcd (a:Arbint.int) b = + if Arbintcore.< (a, b) then gcd b a + else if (Arbintcore.compare (b, Arbintcore.zero) = EQUAL) then a + else gcd b (Arbintcore.mod (a, b)) + +fun scale p = + if not (type_of p = mk_type ("list", [mk_type ("real", [])])) then p + else + let + val (tms, ty) = listSyntax.dest_list p + val tms_mapped = List.map (fn t => if (realSyntax.is_div t) then snd (realSyntax.dest_div t) else t) tms + val tms_filtered = List.filter (fn t => realSyntax.is_real_literal t andalso (not (realSyntax.int_of_term t = Arbint.zero))) tms_mapped + val tms_arbints = List.map realSyntax.int_of_term tms_filtered + val gcd_list = List.foldl (fn (a,b) => + let val res = gcd a b in + if (Arbintcore.compare (res,Arbintcore.zero) = EQUAL) + then b + else (print (Arbintcore.toString res); res) end) (hd tms_arbints) (tl tms_arbints) + in + Parse.Term ‘^(realSyntax.term_of_int gcd_list) *c ^p’ |> EVAL |> rhs o concl + end; + (* +val p = “[-6494149570753833533576226955119636841118335723877 / + 12500000000000000000000000000000000000000000000000; + 13843969975321583631977517825362156145274639129639 / + 12500000000000000000000000000000000000000000000000; + -12955340407972480536169523901435240986756980419159 / + 15625000000000000000000000000000000000000000000000; + 349184292465257220192120790613898861920461058616637 / + 1500000000000000000000000000000000000000000000000000; 0; + -1 / 120; 0; 1 / 5040; 0; -1 / 362880; 0; 1 / 39916800; 0; + -1 / 6227020800; 0; 1 / 1307674368000; 0; -1 / 355687428096000]:real list†+ +val q = “[-18671393139647857184471035907336045056581497192383 / + 100000000000000000000000000000000000000000000000000; + 110400542702815170070795858237033826299011707305909 / + 225000000000000000000000000000000000000000000000000; + -13843969975321583631977517825362156145274639129639 / + 28125000000000000000000000000000000000000000000000; + 1439482267552497837352169322381693442972997824351 / + 6250000000000000000000000000000000000000000000000; + -2444290047256800541344845534297292033443227410316459 / + 54000000000000000000000000000000000000000000000000000; 0; + 1 / 1080; 0; -1 / 72576; 0; 1 / 8164800; 0; -1 / 1437004800; + 0; 1 / 392302310400; 0; -1 / 188305108992000]:real list†+val deg = EVAL “deg ^p - 1†+ +val scaleFactor = “inv (10 pow 10)†+val pScale = EVAL “^scaleFactor *c (12500000000000000000000000000000000000000000000000:real *c ^p)†|> rhs o concl +val qScale = EVAL “^scaleFactor *c (12500000000000000000000000000000000000000000000000:real *c ^q)†|> rhs o concl + +val ss1 = EVAL “sturm_seq_step ^p ^q†(* 3.31 s *) +val ss1_term = ss1 |> rhs o concl |> optionSyntax.dest_some +val ss1Scaled = EVAL “sturm_seq_step ^pScale ^qScale†(* 2.65 s *) +val ss1Scaled_term = ss1Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^scaleFactor *c (inv 12500000000000000000000000000000000000000000000000 *c ^t)†|> rhs o concl) + (* 0.56 s *) + +val ss2 = EVAL “sturm_seq_step ^q ^ss1_term†(* 4.18 s *) +val ss2_term = ss2 |> rhs o concl |> optionSyntax.dest_some +val qScale2 = EVAL “(100000000000000000000000000000000000000000000000000:real *c ^q)†|> rhs o concl +val ss1Scale = EVAL “(100000000000000000000000000000000000000000000000000:real *c ^ss1Scaled_term)†|> rhs o concl +val ss2Scaled = EVAL “sturm_seq_step ^qScale2 ^ss1Scale†(* 3.01 s *) +val ss2Scaled_term = ss2Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv 100000000000000000000000000000000000000000000000000 *c ^t†|> rhs o concl) (* 0.30 s *) + +val ss3 = EVAL “sturm_seq_step ^ss1_term ^ss2_term†(* 5.20 s *) +val ss3_term = ss3 |> rhs o concl |> optionSyntax.dest_some +val tF =“125000000000000000000000000000000000000000000000000000000000:realâ€; +val ss1Scale = EVAL “(^tF *c ^ss1Scaled_term)†|> rhs o concl +val ss2Scale = EVAL “(^tF *c ^ss2Scaled_term)†|> rhs o concl +val ss3Scaled = EVAL “sturm_seq_step ^ss1Scale ^ss2Scale†(* 4.13 s *) +val ss3Scaled_term = ss3Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) (* 0.45 s *) + +val ss4 = EVAL “sturm_seq_step ^ss2_term ^ss3_term†(* 6.57 s *) +val ss4_term = ss4 |> rhs o concl |> optionSyntax.dest_some +val tF =“120000000000000000000000000000000000000000000000000:realâ€; +val ss2Scale = EVAL “(^tF *c ^ss2Scaled_term)†|> rhs o concl +val ss3Scale = EVAL “(^tF *c ^ss3Scaled_term)†|> rhs o concl +val ss4Scaled = EVAL “sturm_seq_step ^ss2Scale ^ss3Scale†(* 3.97 s *) +val ss4Scaled_term = ss4Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) (* 0.37 s *) + +val ss5 = EVAL “sturm_seq_step ^ss3_term ^ss4_term†(* 8.40 s *) +val ss5_term = ss5 |> rhs o concl |> optionSyntax.dest_some +val tF =“125000000000000000000000000000000000000000000000000000000000:realâ€; +val ss3Scale = EVAL “(^tF *c ^ss3Scaled_term)†|> rhs o concl +val ss4Scale = EVAL “(^tF *c ^ss4Scaled_term)†|> rhs o concl +val ss5Scaled = EVAL “sturm_seq_step ^ss3Scale ^ss4Scale†(* 5.64 s *) +val ss5Scaled_term = ss5Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) + +val ss6 = EVAL “sturm_seq_step ^ss4_term ^ss5_term†(* 10.19 s *) +val ss6_term = ss6 |> rhs o concl |> optionSyntax.dest_some +val tF =“100000000000000000000000000000000000000000000000000:realâ€; +val ss4Scale = EVAL “(^tF *c ^ss4Scaled_term)†|> rhs o concl +val ss5Scale = EVAL “(^tF *c ^ss5Scaled_term)†|> rhs o concl +val ss6Scaled = EVAL “sturm_seq_step ^ss4Scale ^ss5Scale†(* 5.83 s *) +val ss6Scaled_term = ss6Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) + +val ss7 = EVAL “sturm_seq_step ^ss5_term ^ss6_term†(* 51.43 s *) +val ss7_term = ss7 |> rhs o concl |> optionSyntax.dest_some +val tF =“1875817730304000000000000000000000000000000000000000000000000000:realâ€; +val ss5Scale = EVAL “(^tF *c ^ss5Scaled_term)†|> rhs o concl +val ss6Scale = EVAL “(^tF *c ^ss6Scaled_term)†|> rhs o concl +val ss7Scaled = EVAL “sturm_seq_step ^ss5Scale ^ss6Scale†(* 37.62 s *) +val ss7Scaled_term = ss7Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) + +val ss8g = EVAL “sturm_seq_step ^ss6_term ^ss7_term†(* 325.40 s *) +val ss8g_term = ss7 |> rhs o concl |> optionSyntax.dest_some +val tF =“52359382622366919003859256612525309952421121811441798729857138480169343485333963069940040914451717670615216237357336816000000000000000000000000000000000000000000000000000:realâ€; +val ss6Scale = EVAL “(^tF *c ^ss6Scaled_term)†|> rhs o concl +val ss7Scale = EVAL “(^tF *c ^ss7Scaled_term)†|> rhs o concl +val ss8gScaled = EVAL “sturm_seq_step ^ss6Scale ^ss7Scale†(* 278.34 s *) +val ss8gScaled_term = ss7Scaled |> rhs o concl |> optionSyntax.dest_some + |> (fn t => EVAL “inv ^tF *c ^t†|> rhs o concl) +*) + Definition sturm_seq_aux_def: sturm_seq_aux (0:num) (p:poly) (q:poly) = (if (rm p (inv (coeff q (deg q)) *c q) = [] ∧ ~ zerop q) @@ -36,80 +164,6 @@ Definition sturm_seq_def: | _ => NONE End -Theorem poly_neg_evals: - ∀p x. poly (--p p) x = - poly p x -Proof - gs[GSYM poly_compat, eval_simps] -QED - -Theorem poly_nill_left_add: - ∀ p x. poly ([] +p p) x = poly p x -Proof - gs[poly_add_rid] >> gs[GSYM poly_compat, reduce_preserving] -QED - -Theorem poly_reduce_evals: - ∀ p x. poly (reduce p) x = poly p x -Proof - gs[GSYM poly_compat, reduce_preserving] -QED - -Theorem poly_neg_neg_evals: - ∀ t. poly (--p (--p t)) = poly t -Proof - rpt strip_tac >> rewrite_tac[FUN_EQ_THM, GSYM poly_compat, eval_simps] - >> real_tac -QED - -Theorem poly_add_aux_evals: - ∀p t x. poly (poly_add_aux p t) x = poly p x + poly t x -Proof - Induct_on ‘p’ - >- gs[poly_add_aux_def, poly_def] - >> rpt strip_tac >> gs[poly_add_aux_def] - >> Cases_on ‘t’ - >- gs[poly_add_aux_lid, poly_def] - >> gs[poly_def] - >> ‘h + x * poly p x + (h' + x * poly t' x) = h + h' + - x * (poly p x + poly t' x)’ by REAL_ARITH_TAC - >> pop_assum $ rewrite_tac o single -QED - -Theorem poly_sub_evals: - ∀ p q x. poly (p -p q) x = poly p x - poly q x -Proof - gs[GSYM poly_compat, eval_simps] -QED - -Theorem poly_cst_evals: - ∀ p c x. c * poly p x = poly (c *c p) x -Proof - gs[GSYM poly_compat, eval_simps] -QED - -Theorem poly_shift_eval: - ∀p q. q ≠[] ⇒ p ≠[] ⇒ deg q < deg p ⇒ - poly (monom (deg p − deg q) [LAST p / LAST q]) x * poly q x = - poly (LAST p / LAST q *c monom (deg p − deg q) q) x -Proof - rpt strip_tac - >> assume_tac LESS_ADD - >> pop_assum $ qspecl_then [‘deg p’, ‘deg q’] assume_tac - >> res_tac - >> ‘deg p - deg q = p'’ by gs[] - >> pop_assum $ rewrite_tac o single - >> pop_assum kall_tac - >> pop_assum kall_tac - >> pop_assum kall_tac - >> Induct_on ‘p'’ - >- ( - gs[monom_def, poly_def] - >> gs[poly_cst_evals, real_div] - ) - >> gs[monom_def, poly_def, GSYM poly_cst_evals] - >> gs[poly_def, GSYM REAL_MUL_ASSOC] >> REAL_ARITH_TAC -QED - (* a / b = c where c * b + r = a*) (* b divides (a + k * r) ⇔ ∃ c. (a + k * r) = b * c *) (* We say p2 divides p1 if ∃ q. p1 * q = p2 *) diff --git a/sturmScript.sml b/sturmScript.sml index a85dfdfb02dadb6077737f88e20662ca88c1faab..9563d510075bbf7f2d5144a99f30dde39ca24338 100644 --- a/sturmScript.sml +++ b/sturmScript.sml @@ -1,6 +1,6 @@ open pred_setTheory listTheory bossLib RealArith realTheory polyTheory; open renameTheory; -open preamble; +open preambleDandelion; val _ = new_theory "sturm"; diff --git a/transcApproxSemScript.sml b/transcApproxSemScript.sml index 34a2bc6a988ccf1ec47385ec8efee62af4f0b047..ba8affe43b3f80932ed89d02b2032162099424bd 100644 --- a/transcApproxSemScript.sml +++ b/transcApproxSemScript.sml @@ -2,7 +2,7 @@ open realTheory realLib RealArith transcTheory; open IntervalArithTheory ErrorValidationTheory sqrtApproxTheory; open realPolyTheory transcLangTheory approxPolyTheory transcIntvSemTheory approxCompErrTheory; -open preamble; +open preambleDandelion; val _ = new_theory "transcApproxSem"; @@ -41,12 +41,16 @@ End Definition errorPropSinCos_def: errorPropSinCos cfg iv err = - do - assert (approxPolySideCond cfg.steps); - (polyCos, errCos) <- approxPoly Cos (0,err) cfg.steps; - (polySin, errSin) <- approxPoly Sin (0,err) cfg.steps; - return (abs ((evalPoly polyCos err - errCos) - 1) + evalPoly polySin err + errSin); - od + if approxPolySideCond cfg.steps + then + case approxPoly Cos (0,err) cfg.steps of + | NONE => NONE + | SOME (polyCos, errCos) => + case approxPoly Sin (0,err) cfg.steps of + | NONE => NONE + | SOME (polySin, errSin) => + SOME (abs ((evalPoly polyCos err - errCos) - 1) + evalPoly polySin err + errSin) + else NONE End Definition errorPropFun_def: @@ -57,16 +61,14 @@ Definition errorPropFun_def: * (2 * err)) (* propagated error from f *) ∧ errorPropFun Sin cfg (iv:real#real) (err:real) (pWiden:poly) (errPWiden:real) = - do - propErr <- errorPropSinCos cfg iv err; - return (errPWiden + propErr) - od + (case errorPropSinCos cfg iv err of + | NONE => NONE + | SOME propErr => SOME (errPWiden + propErr)) ∧ errorPropFun Cos cfg (iv:real#real) (err:real) (pWiden:poly) (errPWiden:real) = - do - propErr <- errorPropSinCos cfg iv err; - return (errPWiden + propErr) - od + (case errorPropSinCos cfg iv err of + | NONE => NONE + | SOME propErr => SOME (errPWiden + propErr)) ∧ errorPropFun _ _ _ _ _ _ = NONE (** TODO **) End @@ -200,8 +202,7 @@ Proof >> rpt $ disch_then drule >> gs[]) (* sin *) >- ( - gs[CaseEq"option", errorPropSinCos_def] - >> ntac 2 (Cases_on ‘x'’ >> gs[CaseEq"option"]) + gs[CaseEq"option", CaseEq"prod", errorPropSinCos_def] >> rpt VAR_EQ_TAC >> irule MCLAURIN_SIN_COMPOSITE_ERR >> gs[PULL_EXISTS] >> qexists_tac ‘cfg.steps’ >> gs[] @@ -210,8 +211,7 @@ Proof >> rpt $ disch_then drule >> gs[]) (* cos *) >- ( - gs[CaseEq"option", errorPropSinCos_def] - >> ntac 2 (Cases_on ‘x'’ >> gs[CaseEq"option"]) + gs[CaseEq"option", CaseEq"prod", errorPropSinCos_def] >> rpt VAR_EQ_TAC >> irule MCLAURIN_COS_COMPOSITE_ERR >> gs[PULL_EXISTS] >> qexists_tac ‘cfg.steps’ >> gs[] diff --git a/transcIntvSemScript.sml b/transcIntvSemScript.sml index 0498d01755cab847e2716a34c7be9d12102fb754..4a61ab05304f196414628f31b86c878faa7c7d63 100644 --- a/transcIntvSemScript.sml +++ b/transcIntvSemScript.sml @@ -4,7 +4,7 @@ open realTheory realLib RealArith transcTheory; open IntervalArithTheory sqrtApproxTheory IntervalValidationTheory; open realPolyTheory transcLangTheory approxPolyTheory; -open preamble; +open preambleDandelion; val _ = new_theory "transcIntvSem"; diff --git a/transcReflectScript.sml b/transcReflectScript.sml index 1e0d143b5686837efdffc02d2be0806dd2d97f38..50384ba3412b9f7b411b8c73f1807f4a5a7cbc93 100644 --- a/transcReflectScript.sml +++ b/transcReflectScript.sml @@ -2,7 +2,7 @@ TODO **) open realPolyTheory realPolyProofsTheory transcLangTheory; -open preamble; +open preambleDandelion; val _ = new_theory"transcReflect"; diff --git a/translateScript.sml b/translateScript.sml new file mode 100644 index 0000000000000000000000000000000000000000..7d3439697c540dcd0db3e3949b39a983c34aa812 --- /dev/null +++ b/translateScript.sml @@ -0,0 +1,433 @@ +open RealArith realTheory realLib realSyntax polyTheory; +open realPolyTheory realPolyProofsTheory checkerTheory sturmComputeTheory + transcLangTheory transcIntvSemTheory approxPolyTheory + transcApproxSemTheory transcReflectTheory euclidDivTheory; +open AbbrevsTheory RealSimpsTheory sqrtApproxTheory IntervalArithTheory; +open RatProgTheory basis; +open bossLib preambleDandelion; + +val _ = new_theory "translate"; + +val _ = translation_extends "RatProg"; + +val _ = translate oEL_def; + +(** Copied from FloVer **) +(* translation of real_div *) +val _ = (next_ml_names := ["real_div"]); +val res = translate realTheory.real_div; +val real_div_side = prove( + let val tm = hyp res |> hd |> rand val v = rand tm + in mk_eq(tm,``^v <> 0``) end,EVAL_TAC) + |> update_precondition; +(* / translation of real_div *) + + +(** Translate basic functions **) +val _ = map translate [FACT, reduce_def, zerop_def, deg_def, coeff_def, monom_def, + poly_add_aux_def, poly_add_def, + poly_mul_cst_aux_def, poly_mul_cst_def, poly_mul_def, + poly_neg_def, poly_sub_def, + evalPoly_def, + divmod_aux_def, divmod_def, rm_def]; + +val expPoly_res = translate expPoly_def; +val exppoly_side_def = fetch "-" "exppoly_side_def"; + +Theorem FACT_NOT_ZERO: + ∀ n. FACT n ≠0 +Proof + Induct_on ‘n’ >> gs[FACT] +QED + +Theorem FACT_NOT_ZERO_REAL: + ∀ n. &FACT n ≠0 +Proof + gs[FACT_NOT_ZERO] +QED + +Theorem exppoly_side_true: + ∀ n. exppoly_side n = T +Proof + Induct_on ‘n’ >> simp[Once exppoly_side_def] + >> rpt strip_tac >> ‘FACT n ≠0’ suffices_by gs[] + >> gs[FACT_NOT_ZERO] +QED + +val _ = exppoly_side_true |> update_precondition; + +Definition sturm_seq_aux_alt_def: + sturm_seq_aux_alt n p q = + case n of + | 0 => + if coeff q (deg q) = 0 then NONE else + if (rm p ((coeff q (deg q))â»Â¹ *c q) = []) ∧ ¬zerop q then SOME [] + else NONE + | SUC n => + if zerop q ∨ (LENGTH (reduce q) < 2) ∨ (coeff q (deg q) = 0) then NONE + else + (let + g = --p (rm p (inv (coeff q (deg q)) *c q)) + in + case sturm_seq_aux_alt n q g of + NONE => NONE + | SOME ss => SOME (g::ss)) +End + +(** sturm_seq_aux gets a precond, but we discharge it later **) +val sturm_seq_aux_res = translate sturm_seq_aux_alt_def; +val sturm_seq_aux_side_def = fetch "-" "sturm_seq_aux_alt_side_def" + +Theorem sturm_seq_aux_side_true: + ∀ n p1 p2. + sturm_seq_aux_alt_side n p1 p2 = T +Proof + Induct_on ‘n’ >> simp[Once sturm_seq_aux_side_def] +QED + +val _ = update_precondition sturm_seq_aux_side_true; + +Definition sturm_seq_alt_def: + sturm_seq_alt p q = + if zerop q ∨ deg p ≤ 1 then NONE + else + case sturm_seq_aux_alt (deg p - 1) p q of + NONE => NONE + | SOME sseq => case LAST sseq of + [] => NONE + | [x] => if x ≠0 then SOME sseq else NONE + | x::xs => NONE +End + +val res = translate sturm_seq_alt_def +val sturm_seq_alt_side_def = fetch "-" "sturm_seq_alt_side_def" + +Theorem sturm_seq_alt_side_true: + ∀ p q. sturm_seq_alt_side p q = T +Proof + gs[sturm_seq_alt_side_def] + >> rpt strip_tac + >> ‘sturm_seq_aux (deg p - 1) p q = SOME []’ by cheat + >> imp_res_tac sturm_seq_aux_length >> gs[] +QED + +val _ = update_precondition sturm_seq_alt_side_true + +val _ = map translate [pow] + +val expErrSmall_res = translate expErrSmall_def; +val _ = fetch "-" "experrsmall_side_def" |> REWRITE_RULE [FACT_NOT_ZERO] |> update_precondition + +val expErrBig_res = translate expErrBig_def; +val _ = fetch "-" "experrbig_side_def" |> REWRITE_RULE [FACT_NOT_ZERO] |> update_precondition + +(** TODO: This uses LEAST, i.e. a loop to find the ceiling, maybe replace with a +faster version using division? **) +val NUM_CEILING_res = translate NUM_CEILING_def; +val num_ceiling_side_def = fetch "-" "num_ceiling_side_def" + +Theorem NUM_CEILING_side_true: + ∀ x. num_ceiling_side x = T +Proof + gs[num_ceiling_side_def, FUN_EQ_THM] >> rpt strip_tac + >> qexists_tac ‘clg (x)’ + >> gs[LE_NUM_CEILING] +QED + +val _ = update_precondition NUM_CEILING_side_true + +val _ = map translate [min4_def, max4_def, IVlo_def, IVhi_def] + +(** Taken from FloVer once more **) +fun LET_CONV var_name body = + (UNBETA_CONV body THENC + RATOR_CONV (ALPHA_CONV (mk_var(var_name, type_of body))) THENC + REWR_CONV (GSYM LET_THM)); + +val absIntvUpd_eq = + IntervalArithTheory.absIntvUpd_def + |> SPEC_ALL + |> CONV_RULE (RAND_CONV (LET_CONV "lo1" ``IVlo iv1`` THENC + LET_CONV "lo2" ``IVlo iv2`` THENC + LET_CONV "hi1" ``IVhi iv1`` THENC + LET_CONV "hi2" ``IVhi iv2``)); + +val addInterval_eq = + IntervalArithTheory.addInterval_def + |> REWRITE_RULE [absIntvUpd_eq] +val _ = translate addInterval_eq + +val multInterval_eq = + IntervalArithTheory.multInterval_def + |> REWRITE_RULE [absIntvUpd_eq] +val _ = translate multInterval_eq +(** End of copy **) + +(** newton and getFunIv have a side condition for the second argument being + unequal to 0, this is discharged later when dealing with interpIntv **) +val _ = map translate [getAnn_def, abs, newton_def, + validate_newton_down_def, validate_newton_up_def, + getFunIv_def + |> SIMP_RULE std_ss [internalSteps_def, newtonIters_def], + negateInterval_def, invertInterval_def, + subtractInterval_def, divideInterval_def, + widenInterval_def, + intvBop_def, intvUop_def] + +val evalpolyintv_res = translate evalPolyIntv_def; +val intvbop_side_def = fetch "-" "intvbop_side_def" +val evalpolyintv_side_def = fetch "-" "evalpolyintv_side_def" |> SIMP_RULE std_ss [intvbop_side_def] + +Theorem evalpolyintv_side_true: + ∀ p iv. evalpolyintv_side p iv = T +Proof + Induct_on ‘p’ >> simp[Once evalpolyintv_side_def] +QED + +val _ = update_precondition evalpolyintv_side_true; + +Definition interpIntv_alt_def: + interpIntv_alt (Var x) (env:(string#(real#real)) list) = + (case FIND (λ (y, iv). y = x) env of + | NONE => NONE + | SOME xv => SOME (VarAnn (SND xv) x)) + ∧ + interpIntv_alt (Cst c) env = SOME (CstAnn (c,c) c) ∧ + interpIntv_alt (Uop uop t) env = + (case interpIntv_alt t env of + | NONE => NONE + | SOME r => + if (~ (uop = Inv)) ∨ (SND (getAnn r) < 0 ∨ 0 < FST (getAnn r)) + then SOME (UopAnn (intvUop uop (getAnn r)) uop r) + else NONE) + ∧ + interpIntv_alt (Bop bop t1 t2) env = + (case interpIntv_alt t1 env of + | NONE => NONE + | SOME r1 => + case interpIntv_alt t2 env of + | NONE => NONE + | SOME r2 => + let iv1 = getAnn r1; iv2 = getAnn r2; in + if bop = Div ⇒ (SND iv2 < 0 ∨ 0 < FST iv2) + then SOME (BopAnn (intvBop bop iv1 iv2) bop r1 r2) + else NONE) + ∧ + interpIntv_alt (Fun s t) env = + (case interpIntv_alt t env of + | NONE => NONE + | SOME r => + let iv = getAnn r in + (* Sqrt defined for positive values only *) + (* Tan cannot be done at 0 because we approximate it with sin x/cos x *) + if ((~ (s = Sqrt)) ∨ (0 < FST iv)) ∧ + ((~ (s = Tan)) ∨ (SND iv < 0 ∨ 0 < FST iv)) + then + case getFunIv s iv of + | NONE => NONE + | SOME ivRes => SOME (FunAnn ivRes s r) + else NONE) + ∧ + interpIntv_alt (Poly p t) env = + case interpIntv_alt t env of + | NONE => NONE + | SOME r => + let iv = getAnn r in + SOME (PolyAnn (evalPolyIntv p iv) p r) +End + +Theorem interpIntv_alt_sound: + interpIntv_alt p env = SOME ann ⇒ + interpIntv p env = SOME ann +Proof + cheat (** TODO **) +QED + +val _ = translate interpIntv_alt_def +val interpintv_alt_side_def = fetch "-" "interpintv_alt_side_def" +val getfuniv_side_def = fetch "-" "getfuniv_side_def" +val newton_side_def = fetch "-" "newton_side_def" +val intvuop_side_def = fetch "-" "intvuop_side_def" +val divideinterval_side_def = fetch "-" "divideinterval_side_def" +val invertinterval_side_def = fetch "-" "invertinterval_side_def" + +Theorem newton_side_inductive: + ∀ n x y. + 0 ≤ x ∧ 0 ≤ y ∧ x ≠0 ∧ y ≠0 ⇒ + newton_side n x y +Proof + Induct_on ‘n’ >> simp[Once newton_side_def] + >> rpt strip_tac + >> ‘0 ≤ x / y ’ by gs[real_div, moreRealTheory.REAL_ZERO_LE_MUL1] + >> first_x_assum irule >> rpt conj_tac + >- gs[] + >- ( + ‘x/y ≠0’ by gs[real_div] + >> ‘(y + x / y) ≠0’ by (gs[real_div] >> real_tac) + >> gs[real_div]) + >> gs[] + >> real_tac +QED + +Theorem interpintv_alt_side_true: + ∀ p env. interpintv_alt_side p env = T +Proof + Induct_on ‘p’ >> simp[Once interpintv_alt_side_def] + >- ( + rpt strip_tac >> gs[getfuniv_side_def] + >> ‘FST (getAnn x1) ≤ SND (getAnn x1)’ by cheat + >> rpt strip_tac >> res_tac >> gs[] + >> irule newton_side_inductive >> gs[] + >> rpt conj_tac >> real_tac) + >- ( + rpt strip_tac + >> gs[intvbop_side_def] >> rpt strip_tac >> res_tac + >> gs[divideinterval_side_def, invertinterval_side_def] + >> ‘FST (getAnn x2) ≤ SND (getAnn x2)’ by cheat + >> conj_tac >> real_tac) + >> rpt strip_tac + >> gs[intvuop_side_def] >> rpt strip_tac >> res_tac + >> gs[invertinterval_side_def] + >> ‘FST (getAnn x4) ≤ SND (getAnn x4)’ by cheat + >> conj_tac >> real_tac +QED + +val _ = update_precondition interpintv_alt_side_true + +Definition approxTransc_alt_def: + approxTransc_alt cfg (VarAnn iv s) = SOME (VarAnn 0 s) ∧ + approxTransc_alt cfg (CstAnn iv r) = SOME (CstAnn 0 r) ∧ + approxTransc_alt cfg (BopAnn iv b e1 e2) = + (case approxTransc_alt cfg e1 of + | NONE => NONE + | SOME e1Appr => + case approxTransc_alt cfg e2 of + | NONE => NONE + | SOME e2Appr => + if (b = Div ⇒ + SND (widenInterval (getAnn e2) (getAnn e2Appr)) < 0 ∨ + 0 < FST (widenInterval (getAnn e2) (getAnn e2Appr))) + then + let propError = errorPropBop b (getAnn e1) (getAnn e1Appr) (getAnn e2) (getAnn e2Appr) + in SOME (BopAnn propError b e1Appr e2Appr) + else NONE) + ∧ + approxTransc_alt cfg (UopAnn iv u e) = + (case approxTransc_alt cfg e of + | NONE => NONE + | SOME eAppr => + if (u = Inv ⇒ + SND (widenInterval (getAnn e) (getAnn eAppr)) < 0 ∨ + 0 < FST (widenInterval (getAnn e) (getAnn eAppr))) + then + let propError = errorPropUop u (getAnn e) (getAnn eAppr) + in SOME (UopAnn propError u eAppr) + else NONE) + ∧ + approxTransc_alt cfg (FunAnn iv f e) = + (case approxTransc_alt cfg e of (* recursive call *) + | NONE => NONE + | SOME eAppr => + (* approximate polynomial on widened interval *) + case approxPoly f (widenInterval (getAnn e) (getAnn eAppr) ) cfg.steps of + | NONE => NONE + | SOME (pWiden,errWiden) => + if (approxPolySideCond cfg.steps ∧ getAnn eAppr ≤ inv 2 ∧ 0 ≤ getAnn eAppr) + then + case errorPropFun f cfg (getAnn e) (getAnn eAppr) pWiden errWiden of + | NONE => NONE + | SOME fullError => SOME (PolyAnn fullError pWiden eAppr) + else NONE) + ∧ + (* We do not support partial approximations for now *) + approxTransc_alt cfg (PolyAnn iv p e) = NONE +End + +Theorem approxTransc_alt_sound: + approxTransc_alt cfg p = SOME ann ⇒ + approxTransc cfg p = SOME ann +Proof + cheat +QED + +val _ = map translate [sinPoly_def, cosPoly_def, sinErr_def, cosErr_def]; + +val sinpoly_side_def = fetch "-" "sinpoly_side_def" |> SIMP_RULE std_ss [] +val cospoly_side_def = fetch "-" "cospoly_side_def" |> SIMP_RULE std_ss [] +val sinerr_side_def = fetch "-" "sinerr_side_def" |> SIMP_RULE std_ss [FACT_NOT_ZERO] |> update_precondition +val coserr_side_def = fetch "-" "coserr_side_def" |> SIMP_RULE std_ss [FACT_NOT_ZERO] |> update_precondition + +Theorem sinpoly_side_true: + ∀ n. sinpoly_side n = T +Proof + Induct_on ‘n’ >> simp[Once sinpoly_side_def, FACT_NOT_ZERO] +QED + +val _ = update_precondition sinpoly_side_true + +Theorem cospoly_side_true: + ∀ n. cospoly_side n = T +Proof + Induct_on ‘n’ >> simp[Once cospoly_side_def, FACT_NOT_ZERO] +QED + +val _ = update_precondition cospoly_side_true + +val _ = map translate [minAbsFun_def, maxAbs_def, approxPolySideCond_def, + approxPoly_def] + +val approxpoly_side_def = fetch "-" "approxpoly_side_def" + +val errorpropsincos_side_def = fetch "-" "errorpropsincos_side_def" +val errorpropfun_side_def = fetch "-" "errorpropfun_side_def" +val approxtransc_alt_side_def = fetch "-" "approxtransc_alt_side_def" + +errorPropUop_def, errorPropBop_def, + errorPropSinCos_def, errorPropFun_def, + approxTransc_alt_def] + + + + fun eval t = Parse.Term t |> EVAL + fun getSome t = if optionSyntax.is_some t then optionSyntax.dest_some t else raise ZeroLibErr "Found NONE instead of SOME" + (* extract components from certificate *) + val (transc, poly, eps, iv) = destCert (defTh |> concl |> rhs) + val ivTm = eval ‘if (LENGTH ^iv = 1) then SOME (SND (HD ^iv)) else NONE’ + |> rhs o concl |> getSome + val var = eval ‘if (LENGTH ^iv = 1) then SOME (FST (HD ^iv)) else NONE’ + |> rhs o concl |> getSome + val iv_FST = EVAL “FST ^ivTm†+ val iv_SND = EVAL “SND ^ivTm†+ val approxSideThm = eval ‘approxPolySideCond ^numApprox’ |> SIMP_RULE std_ss [EQ_CLAUSES] + val ivAnnotThm = eval ‘interpIntv ^transc ^iv’ + val ivAnnotTm = ivAnnotThm |> rhs o concl |> getSome + val ivSoundThm = MATCH_MP interpIntv_sound ivAnnotThm + val approxThm = eval ‘approxTransc (^approxCfg with steps := ^numApprox) ^ivAnnotTm’ + val approxTm = approxThm |> rhs o concl |> getSome + val length1Thm = eval ‘LENGTH ^iv = 1’ |> REWRITE_RULE [EQ_CLAUSES] + val approxSoundThm = + MATCH_MP + (MATCH_MP + (MATCH_MP + (REWRITE_RULE [GSYM AND_IMP_INTRO] approxTransc_sound_single) + length1Thm) + ivSoundThm) + approxThm + |> SIMP_RULE std_ss [erase_def, getAnn_def] + val transpThm = eval ‘reflectToPoly (erase (^approxTm)) ^var’ + val reflectOkThm = MATCH_MP reflectSemEquiv transpThm |> REWRITE_RULE [erase_def] + val varEqThm = EVAL “FST (HD ^iv)†+ val ivEqThm = EVAL “SND (HD ^iv)†+ val approxSoundPolyThm = REWRITE_RULE [varEqThm, ivEqThm, reflectOkThm, optionGet_SOME, AND_IMP_INTRO] approxSoundThm + val transpTm = transpThm |> rhs o concl |> getSome + val transpGetThm = Q.ISPEC ‘^(transpThm |> lhs o concl)’ optionGet_def + |> SIMP_RULE std_ss [SimpR “$=â€, transpThm] + val err = Parse.Term ‘getAnn ^approxTm’ |> EVAL |> concl |> rhs + val errorpThm = Parse.Term ‘^transpTm -p ^poly’ |> EVAL + val errorp = errorpThm |> rhs o concl + val polyErrThm = (testSollya(); + REAL_ZERO_CONV (Parse.Term ‘! x. FST (^ivTm) <= x /\ x <= SND (^ivTm) ==> evalPoly ^errorp x <= ^eps - ^err’)) + val polyErrThm_simped = REWRITE_RULE [GSYM errorpThm, eval_simps, GSYM poly_compat, Once $ GSYM transpGetThm, Once $ GSYM iv_FST, transpThm, optionGet_SOME] polyErrThm + val final_thm = MATCH_MP (MATCH_MP REAL_ABS_TRIANGLE_PRE approxSoundPolyThm) polyErrThm_simped + +val _ = export_theory();