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;