1# Usage: This files builds an error computation for the case when both the 2#�index and the exponent are zero. 3 4# sed -f ../maple/TEMPLOG/polynomials.sed -f ../maple/TEMPLOG/log-de_0.sed ../gappa/log-de-index0-E0.gappa | $GAPPA 5 6 7 8# NOTATION CONVENTION 9# Variables that will be replaced with Maple-computed constants begin with an underscore 10# Variables that correspond to double-precision variables in the code begin with a small letter 11# variable eps_xxx denote the relative rounding error of the last machine operation that computes xxx 12 13# Remove the following line when editing this file, to get back the warnings (all the current ones are harmless) 14#@ -Wno-hint-difference -Wno-null-denominator 15# Remark: asking for warnings resuires at least 512MB of memory 16 17@IEEEdouble = float<ieee_64,ne>; 18@IEEEext = float<x86_80,ne>; 19 20 21# polynomial coefficients, computed by Maple 22#c1 = IEEEext(_c1); 23c2 = IEEEext(_c2); 24c3 = IEEEext(_c3); 25c4 = IEEEext(_c4); 26c5 = IEEEext(_c5); 27c6 = IEEEext(_c6); 28c7 = IEEEext(_c7); 29log2h = IEEEext(_log2h); 30log2l = IEEEext(_log2l); 31#r = IEEEext(_rval); 32#logirh = IEEEext(_logirh); 33#logirl = IEEEext(_logirl); 34 35 36# Transcription of the code, NOT using FMA 37 38# We have removed logirh and Elog2 39z2 IEEEext= z*z; 40p67 IEEEext= c6 + z*c7; 41p45 IEEEext= c4 + z*c5; 42p23 IEEEext= c2 + z*c3; 43#p01 IEEEext= z; 44z4 IEEEext= z2*z2; 45p47 IEEEext= p45 + z2*p67; 46p03 IEEEext= z + z2*p23; # suppressed by hand the exact additions of 0 47logz IEEEext= p03 + z4*p47; 48 49#---------- What this code is supposed to approximate 50# Exact mathematical definition of the log 51Mlogz = Log1pz; 52 53############# Now come auxiliary definitions 54 55# Give the polynomial in Estrin form (no need to transcribe the intermediate steps) 56# Since the argument reduction is exact, Mz = z. 57Mz2 = z*z; 58Mz4 = Mz2*Mz2; 59Mz3 = z*z*z; 60Mz5 = z*z*z*z*z; 61Mz6 = z*z*z*z*z*z; 62eps_z2 = (z2 - z*z)/(z*z); 63eps_z4 = (z4 - z2*z2)/(z2*z2); 64t3 = IEEEext(z*c3); 65t5 = IEEEext(z*c5); 66t7 = IEEEext(z*c7); 67t23 = IEEEext(z2*p23); 68t47 = IEEEext(z4*p47); 69t67 = IEEEext(z2*p67); 70eps_t3 = (t3 - z*c3)/(z*c3); 71eps_t5 = (t5 - z*c5)/(z*c5); 72eps_t7 = (t7 - z*c7)/(z*c7); 73eps_t23 = (t23 - z2*p23)/(z2*p23); 74eps_t67 = (t67 - z2*p67)/(z2*p67); 75eps_t47 = (t47 - z4*p47)/(z4*p47); 76eps_p03 = (p03 - (z+t23)) / (z+t23); 77eps_p23 = (p23 - (c2 + t3)) / (c2 + t3); 78eps_p45 = (p45 - (c4 + t5)) / (c4 + t5); 79eps_p67 = (p67 - (c6 + t7)) / (c6 + t7); 80eps_p47 = (p47 - (p45 + t67)) / (p45 + t67); 81eps_logz = (logz - (p03 + t47))/(p03 + t47); 82 83# The following are just for lightening notations 84eps2_z4 = (1+eps_z2)*(1+eps_z2)*(1+eps_z4)-1; 85eps2_p67 = (1+eps_z2)*(1+eps_p67) -1; 86eps2_p47 = (1+eps2_z4)*(1+eps_p47) -1; 87 88PolyLog1pz = z + Mz2*(c2+z*c3) + Mz4*( (c4+z*c5) + Mz2*(c6+z*c7) ); 89 90epsilon = (logz - Mlogz)/Mlogz; 91 92epsilon_approx = (PolyLog1pz - Log1pz)/Log1pz; 93epsilon_round = (logz-PolyLog1pz)/PolyLog1pz; 94 95# The theorem. Remark that the case z=0 has to be proven by hand (fortunately it is easy) 96 97{ 98 (z in [1b-200, _zabsmax] \/ z in [-_zabsmax, -1b-200]) 99/\ epsilon_approx in [-_epsilonApproxQuick, _epsilonApproxQuick] 100/\ (PolyLog1pz - Log1pz) in [-_deltaApproxQuick, _deltaApproxQuick] 101-> 102epsilon in [-1b-63, 1b-63] 103} 104 105epsilon -> epsilon_approx + epsilon_round + epsilon_approx*epsilon_round; 106 107#We need to factor z by using : 108(logz-PolyLog1pz)/PolyLog1pz -> ((logz/z) - PolyLog1pz/z) / (PolyLog1pz/z); 109# Here the denominator will be bound properly: 110PolyLog1pz/z -> 1 + z*(c2+z*c3) + z*z*z*( (c4+z*c5) + z*z*(c6+z*c7) ) ; 111# so all we have is to cancel as much as possible from the numerator. 112# As all the polynomials come close to zero, we do this by forcing appearance of rel error 113 114# (1) logz = IEEEext(p03 + t47). Since p03+t47 comes close to zero, force appearance of relative error 115logz/z -> (1+ (IEEEext(p03 + t47) - (p03 + t47))/(p03 + t47) ) * ((p03 + t47)/z) ; 116# ------------------rel error----------------- --closer to PolyLog1pz/z--- 117 118# Now if we try to print out the rel error above, it is undefined. 119# This is because Gappa is not able to prove that p03+t47 does not contain 0. 120# So the first job is to get a better enclosure of p03+t47, without 0. 121 122# First get a good enclosure on p03 (so that it does not contain zero): same process as (1) 123p03 -> (1+ eps_p03) * (z+t23); 124# Again, the relative error above exists only if z+t23 does not contain 0 125z+t23-> z*(1+t23/z); 126# now we force the relative error of the rounding of t23 127t23/z -> eps_t23 * ((z2/z)*p23) + (z2/z)*p23; 128z2/z -> ( (z2 - z*z)/(z*z) ) * z + z; # that's eps_z2 129# Now p03 is OK, but p03 + t47 still contains 0: rewrite it, first by rewriting t47: 130t47 -> (1 + eps_t47) * (z4*p47) ; 131# now we may rewrite p03+t47, with the purpose to factor z 132p03+t47 -> z*( (1+eps_p03)*(1+t23/z) + (1+eps_t47)*(z4*p47/z)); 133# remains to bound z4*p47/z, a piece of cake (we even recycle z2/z) 134z4*p47/z -> (1 + (IEEEext(z2*z2) - z2*z2)/(z2*z2)) * p47*z2*(z2/z) ; 135 136 137# At this point the relative error in (1) exists at last. Beside we have almost built a good approx of: 138(p03+t47)/z -> (1+eps_p03)*(1+t23/z) + (1+eps_t47)*(z4*p47/z) ; 139# The previous gives a bound of 2^-8 for the approximation error. Let us refine it. 140# Replacing the terms with the best expression we have for them gives: 141(logz-PolyLog1pz)/PolyLog1pz -> 142 ( (1+ eps_logz)*( (1+eps_p03)*(1+t23/z) 143 + (1+eps_t47)*(z4*p47/z) ) 144 - (1 + z*(c2+z*c3) + z*z*z*( (c4+z*c5) + z*z*(c6+z*c7) )) 145 ) / (PolyLog1pz/z); 146# We develop and simplify the logz/z part of the expression in Maple using expand(). Need to set: interface(prettyprint=0); 147#expand( (1+ eps_logz)*( (1+eps_p03)*(1+t23/z) + (1+eps_t47)*(z4*p47/z) )); 148# Then we may cancel the 1s. 149# Then we regroup the next significant terms: t23/z - z*(c2+z*c3), and 150# (also add parentheses around (z4*p47/z) and (t23/z) so that the remainder is bound properly using previous hints) 151(logz-PolyLog1pz)/PolyLog1pz -> 152 ( 153 (t23/z - z*(c2+z*c3)) 154 +((z4*p47/z) - z*z*z*( (c4+z*c5) + z*z*(c6+z*c7))) 155 + (eps_p03 + eps_p03*(t23/z) + + (z4*p47/z)*eps_t47+eps_logz+eps_logz*(t23/z) 156 + eps_logz*eps_p03 + eps_logz*eps_p03*(t23/z) + eps_logz*(z4*p47/z) + eps_logz*(z4*p47/z)*eps_t47 ) 157 ) / (PolyLog1pz/z); 158# A quick check gives a bound of 2^-62 for the third term (sum of epsilons), which is OK for now. 159# Now let us work on this next significant term, starting with the expression we already have for t23/z: 160t23/z - z*(c2+z*c3) -> eps_t23 * ((z2/z)*p23) + ( (z2/z)*p23 - z*(c2+z*c3)) ; 161# ----------small------- 162# insert rel error of p23, and replace z2/z with the hint already given, then simplify and look for the most significant term : 163# expand((eps_z2 * z + z) * (1+eps_p23)*(c2+t3) - z*(c2+z*c3)); 164(z2/z)*p23 - z*(c2+z*c3) -> (z*t3 - z*z*c3) + eps_z2*z*c2 + eps_z2*z*t3 + eps_z2*z*eps_p23*c2 + eps_z2*z*eps_p23*t3 + z*eps_p23*c2 + z*eps_p23*t3 ; 165(z*t3 - z*z*c3) -> z*(z*c3*eps_t3); 166 167# At this point we have a bound of 2^-71 for t23/z - z*(c2+z*c3), perfect. 168 169# Let us attack (z4*p47/z) - z*z*z*( (c4+z*c5) + z*z*(c6+z*c7)) (currently 2^-25 -- interval arithmetic on z^3) 170 171# Rewrite z4 as z2*z2*(1+eps_z4), then inside this rewrite z2 as z*z*(1+eps_z2). 172z4/z -> z*z*z + z*z*z*((1+eps_z2)*(1+eps_z2)*(1+eps_z4) -1); #OK 173# Also rewrite p47 all the way down : we get 174# p47 -> (1+eps_p47)*(p45 + (1+eps_t67)*(z2*p67)); 175# p45 -> (1+eps_p45)*(c4+(1+eps_t5)*(z*c5)); 176# z2*p67 -> (1+eps_z2)*z*z*(1+eps_p67)*(c6+(1+eps_t7)*(c7*z)); 177# Putting together: 178(z4*p47/z) - z*z*z*( (c4+z*c5) + z*z*(c6+z*c7)) -> 179 z*z*z*( (1+eps2_p47)*(((1+eps_p45)*(c4+(1+eps_t5)*(z*c5))) + (1+eps_t67)*(z*z*(1+eps2_p67)*(c6+(1+eps_t7)*c7*z))) - c4 - z*c5 - z*z*(c6+z*c7)) ; 180 181 182# This last hint is computed by one expand() and one collect(,z) in maple 183 (1+eps2_p47)*(((1+eps_p45)*(c4+(1+eps_t5)*(z*c5))) + (1+eps_t67)*(z*z*(1+eps2_p67)*(c6+(1+eps_t7)*c7*z))) - c4 - z*c5 - z*z*(c6+z*c7) -> 184(eps2_z4*c7*eps_t7+eps2_z4*eps2_p67*c7+c7*eps_t7+eps2_z4*c7+eps2_p67*c7+eps2_p67*c7*eps_t7+eps2_z4*eps2_p67*c7*eps_t7+eps_t67* 185c7+eps_p47*c7+eps_t67*c7*eps_t7+eps_p47*c7*eps_t7+eps_p47*eps_t67*c7+eps2_z4*eps_p47*c7+eps2_z4*eps_t67*c7+eps_t67*eps2_p67*c7 186+eps_p47*eps2_p67*c7+eps_p47*eps_t67*c7*eps_t7+eps2_z4*eps_p47*c7*eps_t7+eps2_z4*eps_p47*eps_t67*c7+eps2_z4*eps_t67*c7*eps_t7+ 187eps_t67*eps2_p67*c7*eps_t7+eps2_z4*eps_p47*eps2_p67*c7+eps_p47*eps2_p67*c7*eps_t7+eps_p47*eps_t67*eps2_p67*c7+eps2_z4*eps_t67* 188eps2_p67*c7+eps2_z4*eps_p47*eps_t67*c7*eps_t7+eps2_z4*eps_p47*eps2_p67*c7*eps_t7+eps2_z4*eps_p47*eps_t67*eps2_p67*c7+eps2_z4* 189eps_p47*eps_t67*eps2_p67*c7*eps_t7+eps_p47*eps_t67*eps2_p67*c7*eps_t7+eps2_z4*eps_t67*eps2_p67*c7*eps_t7)*Mz3+(eps_t67* 190eps2_p67*c6+eps_p47*eps2_p67*c6+eps_p47*c6+eps2_z4*eps_p47*c6+eps_p47*eps_t67*eps2_p67*c6+eps_t67*c6+eps2_z4*eps_p47*eps2_p67* 191c6+eps2_p67*c6+eps2_z4*eps_t67*eps2_p67*c6+eps2_z4*eps_p47*eps_t67*eps2_p67*c6+eps2_z4*eps_t67*c6+eps2_z4*eps_p47*eps_t67*c6+ 192eps2_z4*c6+eps2_z4*eps2_p67*c6+eps_p47*eps_t67*c6)*Mz2+(eps2_z4*eps_p47*eps_p45*c5+c5*eps_t5+eps2_z4*eps_p45*c5*eps_t5+eps_p47 193*c5+eps_p45*c5+eps_p47*eps_p45*c5+eps_p45*c5*eps_t5+eps_p47*eps_p45*c5*eps_t5+eps_p47*c5*eps_t5+eps2_z4*eps_p47*c5+eps2_z4* 194eps_p45*c5+eps2_z4*eps_p47*eps_p45*c5*eps_t5+eps2_z4*eps_p47*c5*eps_t5+eps2_z4*c5*eps_t5+eps2_z4*c5)*z+eps_p47*c4+eps2_z4* 195eps_p45*c4+eps2_z4*eps_p47*c4+eps2_z4*c4+eps_p47*eps_p45*c4+eps2_z4*eps_p47*eps_p45*c4+eps_p45*c4 196;