1# Usage: This files builds an error computation for one index value 2#only. To get a complete proof you need to run it through gappa for all 3#the possible values of the constants. 4 5# Running the ../maple/log-de.mpl Maple script 6# should create the required sed files and give the command to run. 7# Example sed -f ../maple/TEMPLOG/polynomials.sed -f ../maple/TEMPLOG/log-de_1.sed ../gappa/log-de-accurate.gappa | ~/gappa/src/gappa 8 9# Remove the following line when editing this file, to get back the warnings (all the current ones are harmless) 10#@ -Wno-hint-difference -Wno-null-denominator 11# Remark: asking for warnings resuires at least 512MB of memory 12 13#This is to prevent another warning which might let you think Gappa failed when it eventually succeeds. 14#@ -Wno-dichotomy-failure 15 16# Rounding operators and sequences definition 17@IEEEdouble = float<ieee_64,ne>; 18@IEEEext = float<x86_80,ne>; 19@Add22ext = add_rel<124>; 20@Mul22ext = mul_rel<124>; 21 22 23# This is to tell Gappa that E is an integer. 24E=int<ne>(dummyE); 25 26# polynomial coefficients, computed by Maple 27c14 = IEEEext(_c14h); 28c13 = IEEEext(_c13h); 29c12 = IEEEext(_c12h); 30c11 = IEEEext(_c11h); 31c10 = IEEEext(_c10h); 32c9 = IEEEext(_c9h); 33c8 = IEEEext(_c8h); 34 35c7h = IEEEext(_c7h); 36c7l = IEEEext(_c7l); 37c7=c7h+c7l; 38 39c6h = IEEEext(_c6h); 40c6l = IEEEext(_c6l); 41c6=c6h+c6l; 42 43c5h = IEEEext(_c5h); 44c5l = IEEEext(_c5l); 45c5=c5h+c5l; 46 47c4h = IEEEext(_c4h); 48c4l = IEEEext(_c4l); 49c4=c4h+c4l; 50 51c3h = IEEEext(_c3h); 52c3l = IEEEext(_c3l); 53c3=c3h+c3l; 54 55c2h = IEEEext(_c2h); 56c2l = IEEEext(_c2l); 57c2=c2h+c2l; 58 59c1h = IEEEext(_c1h); 60c1l = IEEEext(_c1l); 61c1=c1h+c1l; 62 63log2h = IEEEext(_log2h); 64log2l = IEEEext(_log2l); 65log2=log2h+log2l; 66logirh = IEEEext(_logirh); 67logirl = IEEEext(_logirl); 68logir=logirh+logirl; 69 70# Transcription of C code 71 72t13 IEEEext= c13 + z*c14; 73t12 IEEEext= c12 + z*t13; 74t11 IEEEext= c11 + z*t12; 75t10 IEEEext= c10 + z*t11; 76t9 IEEEext= c9 + z*t10; 77t8 IEEEext= c8 + z*t9; 78 79#Mul12ext(&p7h, &p7l, z, t8); 80p7=z*t8; 81#Add22ext(&t7h, &t7l, p7h,p7l, c7h,c7l); 82t7=Add22ext(c7,p7); 83 84#FMA22ext(&t6h, &t6l, z,0, t7h,t7l, c6h,c6l); 85p6=Mul22ext(z,t7); 86t6=Add22ext(c6,p6); 87 88#FMA22ext(&t5h, &t5l, z,0, t6h,t6l, c5h,c5l); 89p5=Mul22ext(z,t6); 90t5=Add22ext(c5,p5); 91 92#FMA22ext(&t4h, &t4l, z,0, t5h,t5l, c4h,c4l); 93p4=Mul22ext(z,t5); 94t4=Add22ext(c4,p4); 95 96#FMA22ext(&t3h, &t3l, z,0, t4h,t4l, c3h,c3l); 97p3=Mul22ext(z,t4); 98t3=Add22ext(c3,p3); 99 100#FMA22ext(&t2h, &t2l, z,0, t3h,t3l, c2h,c2l); 101p2=Mul22ext(z,t3); 102t2=Add22ext(c2,p2); 103 104#FMA22ext(&t1h, &t1l, z,0, t2h,t2l, c1h,c1l); 105p1=Mul22ext(z,t2); 106t1=Add22ext(c1,p1); 107 108#FMA22ext(&t0h, &t0l, z,0, t1h,t1l, argredtable[index].logirh, argredtable[index].logirl); 109p0=Mul22ext(z,t1); 110t0=Add22ext(logir,p0); 111 112#Mul22ext(&eh, &el, log2h,log2l, E, 0); 113e=Mul22ext(log2,E); 114 115#Add22ext(prh, prl, eh,el, t0h,t0l); 116logz=Add22ext(e,t0); 117 118#---------- What this code is supposed to approximate 119# Exact mathematical definition of the log 120Mlogz = E*Mlog2 + (Mlogir + Log1pz); 121 122PolyLog1pz = z*(c1+z*(c2+z*(c3+z*(c4+z*(c5+z*(c6+z*(c7+z*(c8+z*(c9+z*(c10+z*(c11+z*(c12+z*(c13+z*c14))))))))))))); 123 124epsilon = (logz - Mlogz)/Mlogz; 125 126# Auxiliary epsilons 127Alogz = e + (logir+p0); # exact, before approximating by the Add22 128eps1 = (Alogz - Mlogz)/Mlogz; 129eps2 = (logz - Alogz)/Alogz; 130 131epspoly=(p0-Log1pz)/Log1pz; 132epsolyRound=(p0-PolyLog1pz)/PolyLog1pz; 133epspolyApprox=(PolyLog1pz - Log1pz)/Log1pz; 134 135{ 136 |E| in [1,1024] 137/\ |z| in [1b-200, _zabsmax] 138/\ log2 - Mlog2 in [-1b-129, 1b-129] 139/\ logir - Mlogir in [-1b-129, 1b-129] 140/\ (PolyLog1pz - Log1pz)/Log1pz in [-_epsilonApproxAccurate, _epsilonApproxAccurate] 141/\ (PolyLog1pz - Log1pz) in [-_deltaApproxAccurate, _deltaApproxAccurate] 142-> 143epsilon in [-1b-119,1b-119] 144} 145 146epspoly -> epsolyRound + epspolyApprox + epsolyRound*epspolyApprox; 147 148Log1pz -> PolyLog1pz - (PolyLog1pz-Log1pz); 149 150#(Alogz2-Mlogz)/Mlogz -> ( (e + (logir +p0)) - (E*Mlog2 + (Mlogir + Log1pz)) ) / Mlogz; 151(Alogz-Mlogz)/Mlogz -> ( (e - Mlog2*E) + (logir-Mlogir) + Log1pz*((p0 - Log1pz)/Log1pz) ) / Mlogz; 152 153 154Mlog2 -> log2 - (log2-Mlog2); 155Mlogir -> logir - (logir - Mlogir); 156 157e-E*Mlog2 -> (e-E*log2) + (E*(log2-Mlog2)) ; 158 159epsilon-> eps1 + eps2 + eps1*eps2; 160 161epsilon $ E;