diff --git a/.gitignore b/.gitignore
index 7add466e870aa718c8d279e246433b0817781229..2f16c11cbe005c852493d1d9943900f5a411ea8e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,6 @@
*.dat
*.sig
*Theory*
+examples/sollya_inputs/*_coeffs_results*
+examples/sollya_inputs/*proof*
diff --git a/checkerScript.sml b/checkerScript.sml
index 771469cd00fc953c6637cd1e26c3a56188207189..ca7649f57e78b5e3178e4c784a1f488bfa942afc 100644
--- a/checkerScript.sml
+++ b/checkerScript.sml
@@ -10,16 +10,43 @@ open preamble;
val _ = new_theory "checker";
-Definition exp_steps_def[simp]:
- exp_steps:num = 10
+Definition approx_steps_def[simp]:
+ approx_steps:num = 10
End
+Theorem exp_err_EVAL_THM = EVAL “inv (&FACT approx_steps * 2 pow (approx_steps - 1))â€
+Theorem exp_poly_cst_EVAL_THM = EVAL “exp_poly approx_stepsâ€
+
Definition exp_err_def:
- exp_err = ^(EVAL “inv (&FACT exp_steps * 2 pow (exp_steps - 1))†|> concl |> rhs)
+ exp_err = ^(exp_err_EVAL_THM |> concl |> rhs)
End
Definition exp_poly_cst:
- exp_poly_cst = ^(EVAL “exp_poly exp_steps†|> concl |> rhs)
+ exp_poly_cst = ^(exp_poly_cst_EVAL_THM |> concl |> rhs)
+End
+
+Triviality one_inv_one:
+ 1 * inv 1 = 1
+Proof
+ gs[REAL_MUL_RINV]
+QED
+
+Triviality mul_neg_one:
+ -1 * x = - x:real
+Proof
+ real_tac
+QED
+
+(** Still needs to be multiplied with x pow n **)
+Theorem cos_err_EVAL_THM = EVAL “inv (&FACT approx_steps)â€
+Theorem cos_poly_cst_EVAL_THM = EVAL “cos_poly approx_steps†|> SIMP_RULE std_ss [one_inv_one, REAL_MUL_LID, mul_neg_one]
+
+Definition cos_err_def:
+ cos_err iv = (max (abs (FST iv)) (abs (SND iv))) pow approx_steps * ^(cos_err_EVAL_THM |> concl |> rhs)
+End
+
+Definition cos_poly_cst_def:
+ cos_poly_cst = ^(cos_poly_cst_EVAL_THM |> concl |> rhs)
End
(**
@@ -31,8 +58,10 @@ Definition approx_def:
approx transc (iv:real#real) :(poly#real) option =
case transc of
| FUN tr VAR =>
- if (iv = (0, 0.5) ∧ tr = "exp") then
+ if (iv = (0, 1 * inv 2) ∧ tr = "exp") then
SOME (exp_poly_cst, exp_err)
+ else if (tr = "cos") then
+ SOME (cos_poly_cst, cos_err iv)
else NONE
| _ => NONE
End
@@ -260,52 +289,91 @@ Theorem checker_soundness:
Proof
rpt strip_tac
>> gs[checker_def, approx_def, CaseEq"option", CaseEq"prod", CaseEq"result",
- CaseEq"transc", interp_def, getFun_def]
+ CaseEq"transc"] >> rpt VAR_EQ_TAC
(* Step 1: Approximate the transcendental fun with its taylor series *)
>> irule REAL_LE_TRANS
- >> qexists_tac ‘abs (exp x - poly exp_poly_cst x) + abs (poly exp_poly_cst x - poly cert.poly x)’
+ >> qexists_tac ‘abs (interp (FUN tr VAR) x - poly transp x) + abs (poly transp x - poly cert.poly x)’
>> conj_tac
(* Approximation using triangle inequality *)
>- (
- ‘exp x - poly cert.poly x = (exp x - poly exp_poly_cst x) + (poly exp_poly_cst x - poly cert.poly x)’
- by real_tac
+ qmatch_goalsub_abbrev_tac ‘abs (transc_fun - poly _ _) ≤ _’
+ >> ‘transc_fun - poly cert.poly x = (transc_fun - poly transp x) + (poly transp x - poly cert.poly x)’
+ by real_tac
>> pop_assum $ rewrite_tac o single
>> irule REAL_ABS_TRIANGLE)
- >> ‘cert.eps = exp_err + (cert.eps - exp_err)’ by real_tac
+ (* Split the error into the error from Taylor series and the rest *)
+ >> ‘cert.eps = err + (cert.eps - err)’ by real_tac
>> pop_assum $ once_rewrite_tac o single
+ (* Split the proof into proving two separate approximations *)
>> irule REAL_LE_ADD2 >> reverse conj_tac
+ (* 1. error between Taylor series and certificate polynomial *)
>- (
gs[GSYM poly_compat, GSYM eval_simps]
>> rewrite_tac [poly_compat]
>> irule validateZeroesLeqErr_sound
- >> qexists_tac ‘diff (exp_poly_cst -p cert.poly)’ >> gs[]
- >> qexists_tac ‘(0, 1/2)’ >> gs[]
+ >> qexists_tac ‘diff (transp -p cert.poly)’ >> gs[]
+ >> qexists_tac ‘cert.iv’ >> gs[]
>> qexists_tac ‘cert.zeroes’ >> gs[]
+ >> ‘FST cert.iv ≤ SND cert.iv’ by real_tac
>> drule checkZeroes_sound
>> disch_then drule >> gs[])
- (* Now prove the Taylor series error
- TODO: Make separate soundness proof *)
- >> qspecl_then [‘x’, ‘exp_steps’] strip_assume_tac MCLAURIN_EXP_LE
+ (* 2. error between transcendental function and Taylor series *)
+ (* TODO: Make separate soundness proof *)
+ >> ‘(cert.iv = (0, 1 * inv 2) ∧ tr = "exp") ∨ tr = "cos"’ by (every_case_tac >> gs[])
+ (* exp function *)
+ >- (
+ gs[interp_def, getFun_def]
+ >> qspecl_then [‘x’, ‘approx_steps’] strip_assume_tac MCLAURIN_EXP_LE
+ >> pop_assum $ rewrite_tac o single
+ >> ‘poly exp_poly_cst x = evalPoly (exp_poly approx_steps) x’
+ by (gs[poly_compat] >> EVAL_TAC)
+ >> pop_assum $ rewrite_tac o single
+ >> gs[exp_sum_to_poly]
+ >> qmatch_goalsub_abbrev_tac ‘abs (exp_taylor + taylor_rem - exp_taylor) ≤ _’
+ >> ‘exp_taylor + taylor_rem - exp_taylor = taylor_rem’ by real_tac
+ >> pop_assum $ rewrite_tac o single
+ >> unabbrev_all_tac
+ >> ‘exp_err = inv (&FACT approx_steps * 2 pow (approx_steps - 1))’ by EVAL_TAC
+ >> pop_assum $ rewrite_tac o single
+ >> qspecl_then [‘approx_steps’, ‘x’,‘t’] mp_tac exp_remainder_bounded
+ >> impl_tac >> gs[]
+ >> real_tac)
+ (* cos function *)
+ >> gs[interp_def, getFun_def] >> rpt VAR_EQ_TAC
+ >> qspecl_then [‘x’, ‘approx_steps’] strip_assume_tac MCLAURIN_COS_LE
+ >> gs[]
>> pop_assum $ rewrite_tac o single
- >> ‘poly exp_poly_cst x = evalPoly (exp_poly exp_steps) x’
- by (gs[poly_compat] >> EVAL_TAC)
+ >> ‘poly cos_poly_cst x = evalPoly (cos_poly approx_steps) x’
+ by (rewrite_tac [GSYM approx_steps_def, cos_poly_cst_EVAL_THM]
+ >> gs[poly_compat, cos_poly_cst_def])
>> pop_assum $ rewrite_tac o single
- >> gs[exp_sum_to_poly]
- >> qmatch_goalsub_abbrev_tac ‘abs (exp_taylor + taylor_rem - exp_taylor) ≤ _’
- >> ‘exp_taylor + taylor_rem - exp_taylor = taylor_rem’ by real_tac
+ >> gs[cos_sum_to_poly]
+ >> qmatch_goalsub_abbrev_tac ‘abs (cos_taylor + taylor_rem - cos_taylor) ≤ _’
+ >> ‘cos_taylor + taylor_rem - cos_taylor = taylor_rem’ by real_tac
>> pop_assum $ rewrite_tac o single
>> unabbrev_all_tac
- >> ‘exp_err = inv (&FACT exp_steps * 2 pow (exp_steps - 1))’ by EVAL_TAC
+ >> ‘(x pow approx_steps) * cos t * inv (&FACT approx_steps) =
+ (cos t * ((x pow approx_steps) * inv (&FACT approx_steps)))’
+ by real_tac
+ >> ‘-(x pow approx_steps) * cos t * inv (&FACT approx_steps) =
+ -(cos t * ((x pow approx_steps) * inv (&FACT approx_steps)))’
+ by real_tac
+ >> rewrite_tac [GSYM approx_steps_def]
+ >> ntac 2 $ pop_assum $ rewrite_tac o single
+ >> rewrite_tac [ABS_NEG, Once ABS_MUL]
+ >> qmatch_goalsub_abbrev_tac ‘_ * err_cos_concr’
+ >> irule REAL_LE_TRANS
+ >> qexists_tac ‘ 1 * err_cos_concr’ >> conj_tac
+ >- (irule REAL_LE_RMUL_IMP >> unabbrev_all_tac >> gs[COS_BOUND, ABS_POS])
+ >> unabbrev_all_tac
+ >> rewrite_tac [REAL_MUL_LID, cos_err_def, ABS_MUL]
+ >> ‘abs (inv (&FACT approx_steps)) = inv (&FACT approx_steps)’
+ by (rewrite_tac[abs] >> EVAL_TAC >> gs[])
>> pop_assum $ rewrite_tac o single
- >> qspecl_then [‘exp_steps’, ‘x’,‘t’] mp_tac exp_remainder_bounded
- >> impl_tac >> gs[] >> reverse conj_tac
- >- gs[REAL_INV_1OVER]
- >> gs[abs] >> Cases_on ‘0 ≤ t’
- >- (
- irule REAL_LE_TRANS >> qexists_tac ‘x’ >> gs[REAL_INV_1OVER])
- >> irule REAL_LE_TRANS >> qexists_tac ‘-t’ >> conj_tac >- real_tac
- >> irule REAL_LE_TRANS >> qexists_tac ‘x’
- >> gs[REAL_INV_1OVER]
+ >> rewrite_tac[EVAL “(inv (&FACT approx_steps))â€, ABS_MUL]
+ >> irule REAL_LE_RMUL_IMP >> gs[GSYM POW_ABS]
+ >> irule POW_LE >> gs[ABS_POS]
+ >> irule RealSimpsTheory.maxAbs >> gs[]
QED
val _ = export_theory();
diff --git a/examples/mockUpScript.sml b/examples/mockUpScript.sml
index 7b3c341db231245cd8d650503a54903d894aba35..330b6935f55fea2fe0b6bd14522a99cafcd9cf38 100644
--- a/examples/mockUpScript.sml
+++ b/examples/mockUpScript.sml
@@ -1,5 +1,5 @@
-open realTheory realLib RealArith stringTheory;
-open renameTheory realPolyTheory checkerDefsTheory;
+open realTheory realLib RealArith stringTheory ;
+open renameTheory realPolyTheory checkerDefsTheory checkerTheory sturmComputeTheory;
open preamble;
val _ = new_theory "mockUp";
@@ -12,18 +12,111 @@ val _ = new_theory "mockUp";
Definition exp_example_def:
exp_example =
<|
- transc := (\ x. exp(x) );
- poly := [
- 36028796527503441 * inv (2 pow 55); (* c *)
- 2305847485481749147 * inv (2 pow 61); (* x *)
- 18445094142862447103 * inv (2 pow 65); (* x^2 *)
- 6162767460358145607 * inv (2 pow 65); (* x^3 *)
- 2970172589018140341 * inv (2 pow 66); (* x^4 *)
- 12670224417430857765 * inv (2 pow 70) (* x^5 *)
- ];
- omega := [];
- eps := 8443227763790560011 * inv (2 pow 89) ;
+ transc := FUN "cos" VAR;
+poly := [
+ 4289449735 * inv ( 2 pow 32 );
+ 139975391 * inv ( 2 pow 33 );
+ -2408138823 * inv ( 2 pow 32 );
+ 2948059219 * inv ( 2 pow 35 );
+ ];
+ eps := 582015 * inv (2 pow 31 );
+ zeroes := [
+ ( 1966918861 * inv (2 pow 33 ),
+ 3935725159 * inv (2 pow 34 ));
+ ( 2335336039 * inv (2 pow 32 ),
+ 1167903949 * inv (2 pow 31 ));
+ ( 1856818381 * inv (2 pow 31 ),
+ 3714108621 * inv (2 pow 32 ));
+ ];
+ iv :=
+ ( 858993459 * inv (2 pow 33 ),
+ 1 * inv (2 pow 0 ));
+ omega := [];
|>
End
+(* Definition narrowStep_def: *)
+(* narrowStep (l,h):(real#real) : (real#real) = *)
+(* let dist = (h - l); *)
+(* quart = dist / 8; *)
+(* lN = l + quart; *)
+(* hN = h - quart; *)
+(* in *)
+(* (lN,hN) *)
+(* End *)
+
+(* Definition doNarrow_def: *)
+(* doNarrow (l,h) p = *)
+(* if poly p l * poly p h ≤ 0 then *)
+(* let (lN, hN) = narrowStep (l,h) *)
+(* in *)
+(* if poly p lN * poly p hN ≤ 0 then *)
+(* (lN, hN) *)
+(* else (l,h) *)
+(* else (l,h) *)
+(* End *)
+
+(* Definition doNarrowExplicit_def: *)
+(* doNarrowExplicit (l,h) p = *)
+(* let (lN, hN) = doNarrow (l,h) p in *)
+(* if poly p l * poly p h ≤ 0 then *)
+(* let (lN, hN) = narrowStep (l,h) *)
+(* in *)
+(* if poly p lN * poly p hN ≤ 0 then *)
+(* SOME (SOME (lN, hN)) *)
+(* else SOME NONE *)
+(* else NONE *)
+(* End *)
+
+(*
+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 (transp, err) = EVAL “approx exp_example.transc exp_example.iv†|> concl |> rhs |> optionSyntax.dest_some |> pairSyntax.dest_pair
+val errorp = Parse.Term ‘^transp -p exp_example.poly’ |> EVAL |> rhs o concl
+val deriv1 = Parse.Term ‘diff ^errorp’ |> EVAL |> rhs o concl
+
+(* val test = Parse.Term ‘MAP (λ iv. doNarrowExplicit iv ^deriv1) exp_example.zeroes’ |> EVAL |> concl |> rhs *)
+(* val test2 = Parse.Term ‘MAP (λ ivopt. case ivopt of |SOME (SOME iv) => doNarrowExplicit iv ^deriv1 |_ => NONE) ^test’ |> EVAL |> rhs o concl *)
+(* val test3 = Parse.Term ‘MAP (λ ivopt. case ivopt of |SOME (SOME iv) => doNarrowExplicit iv ^deriv1 |_ => NONE) ^test2’ |> EVAL |> rhs o concl *)
+(* val test4 = Parse.Term ‘MAP (λ ivopt. case ivopt of |SOME (SOME iv) => doNarrowExplicit iv ^deriv1 |_ => NONE) ^test3’ |> EVAL |> rhs o concl *)
+(* val test4 = Parse.Term ‘MAP (λ ivopt. case ivopt of |SOME (SOME iv) => doNarrowExplicit iv ^deriv1 |_ => NONE) ^test4’ |> EVAL |> rhs o concl *)
+(* val test5 = Parse.Term ‘(MAP (λ ivopt. case ivopt of |SOME (SOME iv) => iv |_ => (0,0)) ^test4)’ |> EVAL |> rhs o concl *)
+
+(* val e1 = Parse.Term ‘getMaxWidth exp_example.zeroes’ |> EVAL |> rhs o concl *)
+(* val e2 = Parse.Term ‘getMaxWidth ^test5’ |> EVAL |> rhs o concl *)
+
+val sseq_aux = decls "sturm_seq_aux";
+val _ = computeLib.monitoring := SOME (same_const (hd sseq_aux));
+
+val sseq = Parse.Term ‘sturm_seq ^deriv1 (diff ^deriv1)’
+ |> EVAL
+ |> rhs o concl
+ |> optionSyntax.dest_some
+
+val res = Parse.Term ‘checkZeroes ^deriv1 (diff ^deriv1) exp_example.iv ^sseq exp_example.zeroes’
+ |> EVAL
+
+val res = Parse.Term ‘validateZeroesLeqErr ^errorp exp_example.iv exp_example.zeroes (exp_example.eps - ^err)’ |> EVAL
+
+val mAbs = Parse.Term ‘max (abs (FST exp_example.iv)) (abs (SND exp_example.iv))’ |> EVAL |> rhs o concl
+val B = Parse.Term ‘poly (MAP abs (diff ^errorp)) ^mAbs’ |> EVAL |> rhs o concl
+val e1 = Parse.Term ‘getMaxWidth ^test5’ |> EVAL |> rhs o concl
+val ub = Parse.Term ‘getMaxAbsLb ^errorp ^test5’ |> EVAL |> rhs o concl
+
+val err1 = Parse.Term ‘^ub + ^B * ^e1’ |> EVAL |> rhs o concl
+val err2 = Parse.Term ‘^err’ |> EVAL |> rhs o concl
+
+val res = Parse.Term ‘abs (poly ^errorp (FST exp_example.iv)) ≤ ^ub + ^B * ^e1’ |> EVAL
+val res2 = Parse.Term ‘abs (poly ^errorp (SND exp_example.iv)) ≤ ^ub + ^B * ^e1’ |> EVAL
+val res3 = Parse.Term ‘validBounds exp_example.iv exp_example.zeroes’ |> EVAL
+val res4 = Parse.Term ‘MAP (λ(u,v). poly (diff ^errorp) u * poly (diff ^errorp) v) exp_example.zeroes’ |> EVAL
+val res5 = Parse.Term ‘recordered (FST exp_example.iv) exp_example.zeroes (SND exp_example.iv)’ |> EVAL
+val res6 = Parse.Term ‘^ub + ^B * ^e1 ≤ (exp_example.eps - ^err)’ |> EVAL
+
+EVAL “checker exp_exampleâ€
+*)
+
val _ = export_theory();
diff --git a/examples/sollya_inputs/cos_approx.sollya b/examples/sollya_inputs/cos_approx.sollya
new file mode 100644
index 0000000000000000000000000000000000000000..77213db39635a7396e69712bb6860eb7ac149ec2
--- /dev/null
+++ b/examples/sollya_inputs/cos_approx.sollya
@@ -0,0 +1,61 @@
+oldDisplay=display;
+display = powers!; //putting ! after a command supresses its output
+
+approxPrec = 32;
+prec = approxPrec;
+deg = 3;
+f = cos(x);
+dom = [0.1; 1];
+p = fpminimax(f, deg, [|prec,prec...|], dom, absolute);
+
+derivativeZeros = findzeros(diff(p-f),dom); // here the derivativeZeros is a list of intervals that guarantee to contain the exact zeros
+//derivativeZeros = inf(dom).:derivativeZeros:.sup(dom);
+maximum=0;
+for t in derivativeZeros do {
+ r = evaluate(abs(p-f), t); // r is an interval, we should take its upper bound
+ if sup(r) > maximum then { maximum=sup(r); argmaximum=t; };
+ if (evaluate(diff(p-f),inf(t)) * evaluate (diff(p-f),sup(t)) <= 0 ) then {
+ print ("Ok zero:");
+ print (" (", mantissa (inf(t)), " * inv (2 pow", -exponent(inf(t)), "),");
+ print (" ", mantissa (sup(t)), " * inv (2 pow", -exponent(sup(t)), "));");
+ };
+};
+print("<|");
+print(" transc := FUN \"", f, "\" VAR;");
+print(" poly := [");
+for i from 0 to degree(p) do{
+ coeff_p = coeff(p, i);
+ print(" ", mantissa (coeff_p), " * inv ( 2 pow ", -exponent(coeff_p), ");");
+};
+print (" ];");
+print(" eps := ", mantissa(maximum), " * inv (2 pow", -exponent(maximum), ");");
+print(" zeroes := [");
+for t in derivativeZeros do{
+ print (" (", mantissa (inf(t)), " * inv (2 pow", -exponent(inf(t)), "),");
+ print (" ", mantissa (sup(t)), " * inv (2 pow", -exponent(sup(t)), "));");
+};
+print (" ];");
+print (" iv := ");
+print (" (", mantissa (inf(dom)), " * inv (2 pow", -exponent(inf(dom)), "),");
+print (" ", mantissa (sup(dom)), " * inv (2 pow", -exponent(sup(dom)), "));");
+print (" omega := [];");
+print("|>");
+
+/* Writing into a file */
+filename = "exp_coeffs_results.txt";
+write("Polynomial coefficients p_O,..., p_n represented as their mantissa and exponent, i.e. p_i=mantissa * 2^exponent\n") > filename;
+//attention, here ">" overwrites the file, use ">>" to append
+for i from 0 to degree(p) do{
+ coeff_p = coeff(p, i);
+ write(mantissa(coeff_p),"\t", exponent(coeff_p), "\n") >> filename;
+};
+write("The absolute error bound\n") >> filename;
+write(mantissa(maximum),"\t", exponent(maximum), "\n") >> filename;
+write(derivativeZeros, "\n") >> filename;
+
+
+/* Some sanity checks to verify that our eps is not smaller than the reliable bounds obtained with dedicated methods */
+display=oldDisplay!;
+print("We got epsilon:", maximum);
+print("Sanity check with supnorm:", supnorm(p, f, dom, absolute, 2^(-40)));
+print("Sanity check with infinity norm:", infnorm(p-f, dom, "proof.txt")); //here the last argument is rather for fun, it gives the proof of the bound, you can remove it.
diff --git a/examples/sollya_inputs/exp_approx2.sollya b/examples/sollya_inputs/exp_approx2.sollya
index 80d49cf295a1d0eb5fe2cb301e4ba5bedccf7e7a..433729cb18a4d76154b8c760d070282884279240 100644
--- a/examples/sollya_inputs/exp_approx2.sollya
+++ b/examples/sollya_inputs/exp_approx2.sollya
@@ -1,28 +1,49 @@
oldDisplay=display;
display = powers!; //putting ! after a command supresses its output
-approxPrec = 54;
-deg = 5;
-f = sin(x);
-dom = [0; 0.5];
-p = fpminimax(f, deg, [|D,D...|], dom, absolute);
+approxPrec = 32;
+prec = approxPrec;
+deg = 3;
+f = exp(x);
+dom = [0; 0.4];
+p = fpminimax(f, deg, [|prec,prec...|], dom, absolute);
-derivativeZeros = findzeros(diff(p-f),dom); // here the derivativeZeros is a list of intervals that guarantee to contain the exact zeros
-derivativeZeros = inf(dom).:derivativeZeros:.sup(dom);
+derivativeZeros = findzeros(diff(p-f),dom); // here the derivativeZeros is a list of intervals that guarantee to contain the exact zeros
+//derivativeZeros = inf(dom).:derivativeZeros:.sup(dom);
maximum=0;
for t in derivativeZeros do {
r = evaluate(abs(p-f), t); // r is an interval, we should take its upper bound
if sup(r) > maximum then { maximum=sup(r); argmaximum=t; };
+ if (evaluate(diff(p-f),inf(t)) * evaluate (diff(p-f),sup(t)) <= 0 ) then {
+ print ("Ok zero:");
+ print (" (", mantissa (inf(t)), " * inv (2 pow", -exponent(inf(t)), "),");
+ print (" ", mantissa (sup(t)), " * inv (2 pow", -exponent(sup(t)), "));");
+ };
};
print("<|");
-print(" transc := (\\ x. ", f, ");");
-print(" poly := ", p, ";");
-print(" eps := ", maximum, ";");
+print(" transc := FUN \"", f, "\" VAR;");
+print(" poly := [");
+for i from 0 to degree(p) do{
+ coeff_p = coeff(p, i);
+ print(" ", mantissa (coeff_p), " * inv ( 2 pow ", -exponent(coeff_p), ");");
+};
+print (" ];");
+print(" eps := ", mantissa(maximum), " * inv (2 pow", -exponent(maximum), ");");
+print(" zeroes := [");
+for t in derivativeZeros do{
+ print (" (", mantissa (inf(t)), " * inv (2 pow", -exponent(inf(t)), "),");
+ print (" ", mantissa (sup(t)), " * inv (2 pow", -exponent(sup(t)), "));");
+};
+print (" ];");
+print (" iv := ");
+print (" (", mantissa (inf(dom)), " * inv (2 pow", -exponent(inf(dom)), "),");
+print (" ", mantissa (sup(dom)), " * inv (2 pow", -exponent(sup(dom)), "));");
+print (" omega := [];");
print("|>");
/* Writing into a file */
filename = "exp_coeffs_results.txt";
-write("Polynomial coefficients p_O,..., p_n represented as their mantissa and exponent, i.e. p_i=mantissa * 2^exponent\n") > filename;
+write("Polynomial coefficients p_O,..., p_n represented as their mantissa and exponent, i.e. p_i=mantissa * 2^exponent\n") > filename;
//attention, here ">" overwrites the file, use ">>" to append
for i from 0 to degree(p) do{
coeff_p = coeff(p, i);
@@ -30,6 +51,7 @@ for i from 0 to degree(p) do{
};
write("The absolute error bound\n") >> filename;
write(mantissa(maximum),"\t", exponent(maximum), "\n") >> filename;
+write(derivativeZeros, "\n") >> filename;
/* Some sanity checks to verify that our eps is not smaller than the reliable bounds obtained with dedicated methods */
diff --git a/mcLaurinApproxScript.sml b/mcLaurinApproxScript.sml
index 2b4240d50c8fe61609bc3dc8737ce640bbef591d..ac6b00decefdb94b7f241e43b26dae6f2ea068bb 100644
--- a/mcLaurinApproxScript.sml
+++ b/mcLaurinApproxScript.sml
@@ -226,6 +226,24 @@ Proof
>> gs[]
QED
+(** Mclaurin series for sqrt fucntion **)
+Theorem MCLAURIN_SQRT_LE:
+ ∀ x n. 0 < x ⇒
+ ∃t. abs(t) ≤ abs (x) ∧
+ (\n. sqrt(1+x) = 1 +
+ sum(1,n)
+ (\m. -1 pow (m-1) *
+ inv(1+x) pow (((2*m - 1) DIV 2))
+ * inv(2 pow m) * &(FACT (m-1)) *
+ inv((2 pow (m-1)) * &(FACT (m-1))) +
+ -1 pow (n-1) *
+ inv((1+t) pow (((2*n - 1) DIV 2)))
+ * inv(2 pow n) * &(FACT (n-1)) *
+ inv ((2 pow (n-1)) * &(FACT (n-1))))) n
+Proof
+ cheat
+QED
+
(*** Prove lemma for bound on exp in the interval, x ∈ [0, 0.5]
based on John Harrison's paper **)
@@ -293,7 +311,7 @@ Proof
QED
Theorem REAL_EXP_MINUS1_BOUND_LEMMA:
- !x. &0 <= x /\ x <= inv(&2) ==> &1 - exp(-x) <= &2 * x
+ ∀x. &0 ≤ x ∧ x ≤ inv(&2) ⇒ &1 - exp(-x) ≤ &2 * x
Proof
REPEAT STRIP_TAC >> REWRITE_TAC[REAL_LE_SUB_RADD]
>> ONCE_REWRITE_TAC[REAL_ADD_SYM]
@@ -393,26 +411,72 @@ Proof
>> gs[REAL_INV_MUL']
QED
-(** Mclaurin series for sqrt fucntion **)
+Definition cos_poly_def:
+ cos_poly 0 = [] ∧
+ cos_poly (SUC n) =
+ if (EVEN n) then
+ cos_poly n ++ [-1 pow (n DIV 2) * inv (&FACT n)]
+ else cos_poly n ++ [0]
+End
-Theorem MCLAURIN_SQRT_LE:
- ∀ x n. 0 < x ⇒
- ∃t. abs(t) ≤ abs (x) ∧
- (\n. sqrt(1+x) = 1 +
- sum(1,n)
- (\m. -1 pow (m-1) *
- inv(1+x) pow (((2*m - 1) DIV 2))
- * inv(2 pow m) * &(FACT (m-1)) *
- inv((2 pow (m-1)) * &(FACT (m-1))) +
- -1 pow (n-1) *
- inv((1+t) pow (((2*n - 1) DIV 2)))
- * inv(2 pow n) * &(FACT (n-1)) *
- inv ((2 pow (n-1)) * &(FACT (n-1))))) n
+Theorem LENGTH_cos_poly:
+ LENGTH (cos_poly n) = n
Proof
- cheat
+ Induct_on ‘n’ >> gs[cos_poly_def]
+ >> cond_cases_tac >> gs[]
QED
+Theorem cos_sum_to_poly:
+ ∀ n x. evalPoly (cos_poly n) x =
+ sum(0,n)
+ (λm.
+ (&FACT m)â»Â¹ * x pow m *
+ if EVEN m then cos 0 * -1 pow (m DIV 2)
+ else sin 0 * -1 pow ((SUC m) DIV 2))
+Proof
+ Induct_on ‘n’ >> gs[sum, evalPoly_def, cos_poly_def]
+ >> cond_cases_tac
+ >> gs[evalPoly_app, COS_0, SIN_0, evalPoly_def, LENGTH_cos_poly]
+QED
-
+Theorem cos_even_remainder_bounded:
+ ∀ n.
+ EVEN n ⇒
+ inv (&FACT n) * cos t * x pow n * -1 pow (n DIV 2) ≤
+ abs(inv (&FACT n) * x pow n * -1 pow (n DIV 2))
+Proof
+ rpt strip_tac
+ >> ‘inv (&FACT n) * x pow n * -1 pow (n DIV 2) = inv (&FACT n) * 1 * x pow n * -1 pow (n DIV 2)’
+ by real_tac
+ >> pop_assum $ rewrite_tac o single
+ >> rewrite_tac[GSYM REAL_MUL_ASSOC]
+ >> once_rewrite_tac[REAL_ABS_MUL]
+ >> ‘0 ≤ inv (&FACT n)’
+ by gs[REAL_LE_INV]
+ >> ‘abs (inv (&FACT n)) = inv (&FACT n)’ by gs[abs]
+ >> pop_assum $ rewrite_tac o single
+ >> irule REAL_LE_LMUL_IMP
+ >> reverse conj_tac >- gs[]
+ >> once_rewrite_tac[REAL_ABS_MUL]
+ >> Cases_on ‘0 ≤ x pow n * -1 pow (n DIV 2)’
+ >- (
+ ‘abs (x pow n * -1 pow (n DIV 2)) = x pow n * -1 pow (n DIV 2)’ by gs[abs]
+ >> pop_assum $ rewrite_tac o single
+ >> irule REAL_LE_RMUL_IMP >> gs[COS_BOUNDS])
+ >> ‘x pow n * -1 pow (n DIV 2) < 0’ by real_tac
+ >> irule REAL_LE_TRANS
+ >> qexists_tac ‘-1 * (x pow n * -1 pow (n DIV 2))’
+ >> conj_tac
+ >- (
+ once_rewrite_tac [REAL_MUL_COMM]
+ >> drule REAL_LE_LMUL_NEG
+ >> disch_then $ qspecl_then [‘cos t’, ‘-1’] $ rewrite_tac o single
+ >> gs[COS_BOUNDS])
+ >> ‘∃ y. x pow n * -1 pow (n DIV 2) = -1 * y ∧ 0 ≤ y’
+ by (qexists_tac ‘-1 * x pow n * -1 pow (n DIV 2)’
+ >> real_tac)
+ >> qpat_x_assum `_ = -1 * y` $ rewrite_tac o single
+ >> gs[ABS_LE]
+QED
val _ = export_theory();