diff --git a/checkerScript.sml b/checkerScript.sml index f6ac4a09da447107340a05186339c52351e87e1f..139db6733c2c22deb51499afc6fc04e03487c3d8 100644 --- a/checkerScript.sml +++ b/checkerScript.sml @@ -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 diff --git a/examples/mockUpCosScript.sml b/examples/mockUpCosScript.sml index 2cb2c3748e7eb3c85494ecebff261055723551e5..f2b3d2f051a8d287d3a65ef5a260caac26ae1e92 100644 --- a/examples/mockUpCosScript.sml +++ b/examples/mockUpCosScript.sml @@ -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 diff --git a/examples/mockUpExpScript.sml b/examples/mockUpExpScript.sml new file mode 100644 index 0000000000000000000000000000000000000000..ce9cc669f8573d848bc131f15f58bae452a1b54a --- /dev/null +++ b/examples/mockUpExpScript.sml @@ -0,0 +1,76 @@ +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(); diff --git a/examples/mockUpSinScript.sml b/examples/mockUpSinScript.sml index c52732b07fd39f96356c237ba64ec3041aac87fe..477c1c2894a9d3f8957240313408d476bcb6ba88 100644 --- a/examples/mockUpSinScript.sml +++ b/examples/mockUpSinScript.sml @@ -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 diff --git a/examples/sollya_inputs/exp_approx.sollya b/examples/sollya_inputs/exp_approx.sollya index f8b4e371e0ceb92a4a9138cda7f8fbd68e1dee98..bed4a9e1f8d00b1e0d0fc942d9f3df9d2290f5be 100644 --- a/examples/sollya_inputs/exp_approx.sollya +++ b/examples/sollya_inputs/exp_approx.sollya @@ -1,17 +1,61 @@ -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. diff --git a/examples/sollya_inputs/exp_approx2.sollya b/examples/sollya_inputs/exp_approx2.sollya deleted file mode 100644 index 433729cb18a4d76154b8c760d070282884279240..0000000000000000000000000000000000000000 --- a/examples/sollya_inputs/exp_approx2.sollya +++ /dev/null @@ -1,61 +0,0 @@ -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.