Skip to content
Snippets Groups Projects
Commit e20ab296 authored by Heiko Becker's avatar Heiko Becker
Browse files

Make checker parametric in Taylor steps, add working exp example

parent b57629e2
No related branches found
No related tags found
No related merge requests found
Pipeline #53634 passed
......@@ -19,11 +19,13 @@ Theorem exp_err_big_EVAL_THM = EVAL “ inv (&FACT approx_steps * 2 pow approx_s
Theorem exp_poly_cst_EVAL_THM = EVAL exp_poly approx_steps
Definition exp_err_small_def:
exp_err_small = ^(exp_err_small_EVAL_THM |> concl |> rhs)
exp_err_small approxSteps = (* ^(exp_err_small_EVAL_THM |> concl |> rhs) *)
inv (&FACT approxSteps * 2 pow (approxSteps - 1))
End
Definition exp_err_big_def:
exp_err_big n = 2 pow n * &n pow approx_steps * ^(exp_err_big_EVAL_THM |> rhs o concl)
exp_err_big n approxSteps = 2 pow n * &n pow approxSteps * (* ^(exp_err_big_EVAL_THM |> rhs o concl)*)
inv (&FACT approxSteps * 2 pow approxSteps)
End
Definition exp_poly_cst:
......@@ -47,7 +49,8 @@ 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)
cos_err iv approxSteps = (max (abs (FST iv)) (abs (SND iv))) pow approxSteps * (* ^(cos_err_EVAL_THM |> concl |> rhs)*)
inv (&FACT approxSteps)
End
Definition cos_poly_cst_def:
......@@ -58,7 +61,8 @@ Theorem sin_err_EVAL_THM = EVAL “inv (&FACT approx_steps)”
Theorem sin_poly_cst_EVAL_THM = EVAL sin_poly approx_steps |> SIMP_RULE std_ss [one_inv_one, REAL_MUL_LID, mul_neg_one]
Definition sin_err_def:
sin_err iv = (max (abs (FST iv)) (abs (SND iv))) pow approx_steps * ^(sin_err_EVAL_THM |> concl |> rhs)
sin_err iv approxSteps = (max (abs (FST iv)) (abs (SND iv))) pow approxSteps * (* ^(sin_err_EVAL_THM |> concl |> rhs) *)
inv (&FACT approxSteps)
End
Definition sin_poly_cst_def:
......@@ -71,23 +75,23 @@ End
**)
(* TODO: Approximate transcendental functions on a case-by-case basis *)
Definition approxPoly_def:
approxPoly transc (iv:real#real) (hs:hint list) :(poly#real) option =
approxPoly transc (iv:real#real) (hs:hint list) approxSteps :(poly#real) option =
case transc of
| FUN tr VAR =>
if tr = "exp" then
case getExpHint hs of
| NONE =>
if iv = (0, 1 * inv 2) then
SOME (exp_poly_cst, exp_err_small)
SOME (exp_poly approxSteps, exp_err_small approxSteps)
else NONE
| SOME n =>
if iv = (0, &n * inv 2) then
SOME (exp_poly_cst, exp_err_big n)
SOME (exp_poly approxSteps, exp_err_big n approxSteps)
else NONE
else if tr = "cos" then
SOME (cos_poly_cst, cos_err iv)
SOME (cos_poly approxSteps, cos_err iv approxSteps)
else if tr = "sin" then
SOME (sin_poly_cst, sin_err iv)
SOME (sin_poly approxSteps, sin_err iv approxSteps)
else NONE
| _ => NONE
End
......@@ -163,8 +167,10 @@ End
Overall certificate checker combines all of the above functions into one that
runs over the full certificate **)
Definition checker_def:
checker (cert:certificate) :result =
case approxPoly cert.transc cert.iv cert.hints of
checker (cert:certificate) approxSteps :result =
if ~ EVEN approxSteps ~ EVEN (approxSteps DIV 2) approxSteps = 0
then Invalid "Need even number of approximation steps"
else case approxPoly cert.transc cert.iv cert.hints approxSteps of
| NONE => Invalid "Could not find appropriate approximation"
| SOME (transp, err) =>
let errorp = transp -p cert.poly;
......@@ -320,18 +326,20 @@ Proof
QED
Theorem checker_soundness:
cert.
checker cert = Valid
cert approxSteps.
checker cert approxSteps = Valid
x.
FST(cert.iv) x x SND (cert.iv)
abs (interp cert.transc x - poly cert.poly x) cert.eps
Proof
rpt strip_tac
>> gs[checker_def, approxPoly_def, CaseEq"option", CaseEq"prod", CaseEq"result",
CaseEq"transc"] >> rpt VAR_EQ_TAC
rpt gen_tac >> gs[checker_def]
>> cond_cases_tac
>> gs[checker_def, approxPoly_def,
CaseEq"option", CaseEq"prod", CaseEq"result", CaseEq"transc"]
>> rpt strip_tac >> rpt VAR_EQ_TAC
(* Step 1: Approximate the transcendental fun with its taylor series *)
>> irule REAL_LE_TRANS
>> qexists_tac abs (interp (FUN tr VAR) x - poly transp x) + abs (poly transp x - poly cert.poly x)
>> qexists_tac abs (interp cert.transc x - poly transp x) + abs (poly transp x - poly cert.poly x)
>> conj_tac
(* Approximation using triangle inequality *)
>- (
......@@ -367,19 +375,18 @@ Proof
(* exp function, 0 to 1/2 *)
>- (
gs[interp_def, getFun_def]
>> qspecl_then [x, approx_steps] strip_assume_tac MCLAURIN_EXP_LE
>> qspecl_then [x, approxSteps] 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)
>> poly transp x = evalPoly (exp_poly approxSteps) x
by (gs[poly_compat] (* >> EVAL_TAC *))
>> pop_assum $ rewrite_tac o single
>> rewrite_tac[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_small = 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_small
>> exp_err_small approxSteps = inv (&FACT approxSteps * 2 pow (approxSteps - 1)) by EVAL_TAC
>> qspecl_then [approxSteps, x,t] mp_tac exp_remainder_bounded_small
>> impl_tac >> gs[]
>> real_tac)
(* exp function, 0 to 1 *)
......@@ -390,33 +397,33 @@ Proof
>> CCONTR_TAC
>> pop_assum $ mp_tac o SIMP_RULE std_ss []
>> rewrite_tac[REAL_INV_INJ] >> real_tac)
>> err = exp_err_big n transp = exp_poly_cst by gs[]
>> err = exp_err_big n approxSteps transp = exp_poly approxSteps by gs[]
>> rpt VAR_EQ_TAC
>> rewrite_tac[GSYM poly_compat, eval_simps]
>> exp_poly_cst = exp_poly approx_steps by EVAL_TAC
(* >> ‘exp_poly_cst = exp_poly approxSteps’ by EVAL_TAC *)
>> pop_assum $ rewrite_tac o single
>> rewrite_tac[exp_sum_to_poly]
>> qspecl_then [x, approx_steps] strip_assume_tac MCLAURIN_EXP_LE
>> qspecl_then [x, approxSteps] strip_assume_tac MCLAURIN_EXP_LE
>> pop_assum $ rewrite_tac o single
>> 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_big n = 2 pow n * &n pow approx_steps * inv (&FACT approx_steps * 2 pow approx_steps)
by (rewrite_tac[approx_steps_def] >> EVAL_TAC)
>> exp_err_big n approxSteps = 2 pow n * &n pow approxSteps * inv (&FACT approxSteps * 2 pow approxSteps)
by (rewrite_tac[] >> EVAL_TAC)
>> pop_assum $ rewrite_tac o single
>> qspecl_then [approx_steps, n, x,t] mp_tac exp_remainder_bounded_big
>> qspecl_then [approxSteps, n, x,t] mp_tac exp_remainder_bounded_big
>> impl_tac
>- (rpt conj_tac >> gs[] >> real_tac)
>> rewrite_tac[approx_steps_def])
>> rewrite_tac[])
(* cos function *)
>- (
gs[interp_def, getFun_def] >> rpt VAR_EQ_TAC
>> qspecl_then [x, approx_steps] strip_assume_tac MCLAURIN_COS_LE
>> qspecl_then [x, approxSteps] strip_assume_tac MCLAURIN_COS_LE
>> gs[]
>> pop_assum $ rewrite_tac o single
>> poly cos_poly_cst x = evalPoly (cos_poly approx_steps) x
by (rewrite_tac [GSYM approx_steps_def, cos_poly_cst_EVAL_THM]
>> poly (cos_poly approxSteps) x = evalPoly (cos_poly approxSteps) x
by (rewrite_tac [cos_poly_cst_EVAL_THM]
>> gs[poly_compat, cos_poly_cst_def])
>> pop_assum $ rewrite_tac o single
>> gs[cos_sum_to_poly]
......@@ -424,35 +431,47 @@ Proof
>> cos_taylor + taylor_rem - cos_taylor = taylor_rem by real_tac
>> pop_assum $ rewrite_tac o single
>> unabbrev_all_tac
>> (x pow approx_steps) * cos t * inv (&FACT approx_steps) =
(cos t * ((x pow approx_steps) * inv (&FACT approx_steps)))
>> (x pow approxSteps) * cos t * inv (&FACT approxSteps) =
(cos t * ((x pow approxSteps) * inv (&FACT approxSteps)))
by real_tac
>> -(x pow approx_steps) * cos t * inv (&FACT approx_steps) =
-(cos t * ((x pow approx_steps) * inv (&FACT approx_steps)))
>> -(x pow approxSteps) * cos t * inv (&FACT approxSteps) =
-(cos t * ((x pow approxSteps) * inv (&FACT approxSteps)))
by real_tac
>> rewrite_tac [GSYM approx_steps_def]
>> rewrite_tac []
>> ntac 2 $ pop_assum $ rewrite_tac o single
>> rewrite_tac [ABS_NEG, Once ABS_MUL]
>> qmatch_goalsub_abbrev_tac _ * err_cos_concr
>> rewrite_tac [GSYM REAL_MUL_ASSOC]
>> qmatch_goalsub_abbrev_tac abs (cos _ * 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])
>> qexists_tac 1 * abs err_cos_concr >> conj_tac
>- (rewrite_tac[ABS_MUL] >> irule REAL_LE_RMUL_IMP >> unabbrev_all_tac >> gs[COS_BOUND, ABS_POS])
>> rewrite_tac[REAL_MUL_LID]
>> abs err_cos_concr = err_cos_concr
by (unabbrev_all_tac
>> rewrite_tac[ABS_REFL]
>> irule REAL_LE_MUL >> conj_tac
>- (irule REAL_LE_INV >> gs[REAL_POS])
>> irule REAL_LE_MUL >> conj_tac
>> gs[REAL_POW_GE0])
>> pop_assum $ rewrite_tac o single
>> unabbrev_all_tac
>> rewrite_tac [REAL_MUL_LID, cos_err_def, ABS_MUL]
>> abs (inv (&FACT approx_steps)) = inv (&FACT approx_steps)
>> rewrite_tac [cos_err_def]
(* >> ‘abs (inv (&FACT approxSteps)) = inv (&FACT approxSteps)’
by (rewrite_tac[abs] >> EVAL_TAC >> gs[])
>> pop_assum $ rewrite_tac o single
>> rewrite_tac[EVAL (inv (&FACT approx_steps)), ABS_MUL, real_div, REAL_MUL_LID]
>> irule REAL_LE_RMUL_IMP >> gs[GSYM POW_ABS]
>> pop_assum $ rewrite_tac o single *)
>> imp_res_tac EVEN_ODD_EXISTS >> gs[POW_MINUS1]
(* >> rewrite_tac[ABS_MUL, real_div, REAL_MUL_LID] *)
>> irule REAL_LE_LMUL_IMP >> gs[GSYM POW_ABS]
>> irule REAL_LE_TRANS
>> qexists_tac abs (x pow (2 * m)) >> gs[ABS_LE, GSYM POW_ABS]
>> irule POW_LE >> gs[ABS_POS]
>> irule RealSimpsTheory.maxAbs >> gs[])
(* sin *)
>> gs[interp_def, getFun_def] >> rpt VAR_EQ_TAC
>> qspecl_then [x, approx_steps] strip_assume_tac MCLAURIN_SIN_LE
>> qspecl_then [x, approxSteps] strip_assume_tac MCLAURIN_SIN_LE
>> gs[]
>> pop_assum $ rewrite_tac o single
>> poly sin_poly_cst x = evalPoly (sin_poly approx_steps) x
by (rewrite_tac [GSYM approx_steps_def, sin_poly_cst_EVAL_THM]
>> poly (sin_poly approxSteps) x = evalPoly (sin_poly approxSteps) x
by (rewrite_tac [sin_poly_cst_EVAL_THM]
>> gs[poly_compat, sin_poly_cst_def])
>> pop_assum $ rewrite_tac o single
>> gs[sin_sum_to_poly]
......@@ -460,26 +479,38 @@ Proof
>> sin_taylor + taylor_rem - sin_taylor = taylor_rem by real_tac
>> pop_assum $ rewrite_tac o single
>> unabbrev_all_tac
>> (x pow approx_steps) * inv (&FACT approx_steps) * sin t =
(sin t * ((x pow approx_steps) * inv (&FACT approx_steps)))
>> inv (&FACT approxSteps) * sin t * x pow approxSteps * -1 pow (approxSteps DIV 2) =
(sin t * ((x pow approxSteps) * inv (&FACT approxSteps) * -1 pow (approxSteps DIV 2)))
by real_tac
>> -(x pow approx_steps) * inv (&FACT approx_steps) * sin t =
-(sin t * ((x pow approx_steps) * inv (&FACT approx_steps)))
>> -(x pow approxSteps) * inv (&FACT approxSteps) * sin t =
-(sin t * ((x pow approxSteps) * inv (&FACT approxSteps)))
by real_tac
>> rewrite_tac [GSYM approx_steps_def]
>> rewrite_tac []
>> ntac 2 $ pop_assum $ rewrite_tac o single
>> rewrite_tac [ABS_NEG, Once ABS_MUL]
>> rewrite_tac[GSYM REAL_MUL_ASSOC]
>> qmatch_goalsub_abbrev_tac _ * err_sin_concr
>> rewrite_tac [ABS_NEG, Once ABS_MUL]
>> irule REAL_LE_TRANS
>> qexists_tac 1 * err_sin_concr >> conj_tac
>> qexists_tac 1 * abs err_sin_concr >> conj_tac
>- (irule REAL_LE_RMUL_IMP >> unabbrev_all_tac >> gs[SIN_BOUND, ABS_POS])
>> unabbrev_all_tac
>> rewrite_tac [REAL_MUL_LID, sin_err_def, ABS_MUL]
>> abs (inv (&FACT approx_steps)) = inv (&FACT approx_steps)
by (rewrite_tac[abs] >> EVAL_TAC >> gs[])
>> abs err_sin_concr = err_sin_concr
by (unabbrev_all_tac
>> rewrite_tac[ABS_REFL]
>> irule REAL_LE_MUL >> conj_tac
>> gs[REAL_POW_GE0]
>> irule REAL_LE_MUL >> gs[REAL_POS, REAL_POW_GE0])
>> pop_assum $ rewrite_tac o single
>> rewrite_tac[EVAL (inv (&FACT approx_steps)), ABS_MUL, real_div, REAL_MUL_LID]
>> irule REAL_LE_RMUL_IMP >> gs[GSYM POW_ABS]
>> unabbrev_all_tac
>> rewrite_tac [sin_err_def]
(* >> ‘abs (inv (&FACT approxSteps)) = inv (&FACT approxSteps)’
by (rewrite_tac[abs] >> EVAL_TAC >> gs[])
>> pop_assum $ rewrite_tac o single *)
>> imp_res_tac EVEN_ODD_EXISTS >> gs[POW_MINUS1]
(* >> rewrite_tac[ABS_MUL, real_div, REAL_MUL_LID] *)
>> irule REAL_LE_LMUL_IMP >> gs[GSYM POW_ABS]
>> irule REAL_LE_TRANS
>> qexists_tac abs (x pow (2 * m)) >> gs[ABS_LE, GSYM POW_ABS]
>> irule POW_LE >> gs[ABS_POS]
>> irule RealSimpsTheory.maxAbs >> gs[]
QED
......
......@@ -38,7 +38,7 @@ 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]
Theorem checkerSucceeds:
checker cos_example = Valid
checker cos_example 16 = Valid
Proof
EVAL_TAC
QED
......
open realTheory realLib RealArith stringTheory ;
open renameTheory realPolyTheory checkerDefsTheory checkerTheory sturmComputeTheory;
open preamble;
val _ = new_theory "mockUpExp";
Definition exp_example_def:
exp_example =
<|
transc := FUN "exp" VAR;
poly := [
9002292208500749 * inv ( 2 pow 53 );
4578369858298151 * inv ( 2 pow 52 );
7596726129252049 * inv ( 2 pow 54 );
1260902011754487 * inv ( 2 pow 52 );
];
eps := 1226771469587 * inv (2 pow 51 );
zeroes := [
( 81979129 * inv (2 pow 29 ),
40989565 * inv (2 pow 28 ));
( 275130831 * inv (2 pow 29 ),
17195677 * inv (2 pow 25 ));
( 461584775 * inv (2 pow 29 ),
57698097 * inv (2 pow 26 ));
];
iv :=
( 0 * inv (2 pow 0 ),
1 * inv (2 pow 0 ));
hints := [ EXP_UB_SPLIT 2];
|>
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]
Theorem checkerSucceeds:
checker exp_example 16 = Valid
Proof
EVAL_TAC
QED
(*
fun debugChecker cert_def =
let
val eval = fn t => Parse.Term t |> EVAL
val certN = cert_def |> concl |> lhs
val approxSteps = Parse.Term ‘14:num’
val (transp, err) = eval ‘approxPoly ^(certN).transc ^(certN).iv ^(certN).hints ^approxSteps’
|> concl |> rhs |> optionSyntax.dest_some |> pairSyntax.dest_pair
val errorp = eval ‘^transp -p ^(certN).poly’ |> rhs o concl
val deriv1 = eval ‘diff ^errorp’ |> rhs o concl
(* uncomment for debugging *)
(* val sseq_aux = decls "sturm_seq_aux";
val _ = computeLib.monitoring := SOME (same_const (hd sseq_aux)); *)
val sseq = eval ‘sturm_seq ^deriv1 (diff ^deriv1)’
|> rhs o concl
|> optionSyntax.dest_some
val res = eval ‘checkZeroes ^deriv1 (diff ^deriv1) ^(certN).iv ^sseq ^(certN).zeroes’
val res = Parse.Term ‘validateZeroesLeqErr ^errorp ^(certN).iv ^(certN).zeroes (^(certN).eps - ^err)’ |> EVAL
val mAbs = Parse.Term ‘max (abs (FST ^(certN).iv)) (abs (SND ^(certN).iv))’ |> EVAL |> rhs o concl
val B = Parse.Term ‘poly (MAP abs (diff ^errorp)) ^mAbs’ |> EVAL |> rhs o concl
val e1 = Parse.Term ‘getMaxWidth ^(certN).zeroes’ |> EVAL |> rhs o concl
val ub = Parse.Term ‘getMaxAbsLb ^errorp ^(certN).zeroes’ |> EVAL |> rhs o concl
val res = Parse.Term ‘abs (poly ^errorp (FST ^(certN).iv)) ≤ ^ub + ^B * ^e1’ |> EVAL
val res2 = Parse.Term ‘abs (poly ^errorp (SND ^(certN).iv)) ≤ ^ub + ^B * ^e1’ |> EVAL
val res3 = Parse.Term ‘validBounds ^(certN).iv ^(certN).zeroes’ |> EVAL
val res4 = Parse.Term ‘MAP (λ(u,v). poly (diff ^errorp) u * poly (diff ^errorp) v) ^(certN).zeroes’ |> EVAL
val res5 = Parse.Term ‘recordered (FST ^(certN).iv) ^(certN).zeroes (SND ^(certN).iv)’ |> EVAL
val res6 = Parse.Term ‘^ub + ^B * ^e1 ≤ (^(certN).eps - ^err)’ |> EVAL
val err1 = Parse.Term ‘^ub + ^B * ^e1’ |> EVAL |> rhs o concl
val err2 = Parse.Term ‘^(certN).eps - ^err’ |> EVAL |> rhs o concl *)
val _ = export_theory();
......@@ -37,7 +37,7 @@ 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]
Theorem checkerSucceeds:
checker sin_example = Valid
checker sin_example 16 = Valid
Proof
EVAL_TAC
QED
......
display = powers;
prec = 64;
oldDisplay=display;
display = powers!; //putting ! after a command supresses its output
diam = 1b-30;
approxPrec = 53;
prec = approxPrec;
deg = 3;
f = exp(x);
d = [0; 0.5];
p = remez(f,5,d);
derivativeZeros = dirtyfindzeros(diff(p-f),d);
derivativeZeros = inf(d).:derivativeZeros:.sup(d);
dom = [0; 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);
if r > maximum then { maximum=r; argmaximum=t; };
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("|>");
\ No newline at end of file
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("|>");
/* 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.
oldDisplay=display;
display = powers!; //putting ! after a command supresses its output
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);
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.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment