1#######################################################################
2# This file is part of the crlibm library, and is distributed under
3# the  LGPL.
4# To use:
5# restart; read "log-de.mpl";
6#TODO output scripts for the Gappa proof out for Paterson/Stockmayer
7Digits := 100:
8
9interface(quiet=true):
10
11read "common-procedures.mpl":
12read "double-extended.mpl":
13mkdir("TEMPLOG"):
14
15
16log2h,log2l := hiloExt(log(2)):
17
18
19L := 7: # number of bits used to address the table
20
21MAXINDEX    := round(2^L * (sqrt(2)-1)):
22
23for i from 0 to MAXINDEX-1 do
24    center[i] := 1 + i*2^(-L): # center[i] in [1, 2[
25    # We want it to fit on 11 bits of mantissa
26    r[i] :=  round(evalf(  (1/center[i]) * 2^(11)) ) / 2^(11) ;
27
28od:
29for i from MAXINDEX to 2^L do
30    # y has been divided by two, center[i] in [0.5, 1[
31    center[i]:=(1 + i*2^(-L)) / 2:
32    # We want it to fit on 11 bits of mantissa,
33    r[i] :=  round(evalf(  (1/center[i]) * 2^(10)) ) / 2^(10) ;
34od:
35
36# Note that we go up to 2^L although the case 2^L is wrapped to zero
37# in the C code. It could be important for zmax (but it turns out not).
38
39for i from 0 to 2^L do
40    (logirh[i], logirl[i]) := hiloExt(-log(r[i])):
41od:
42
43#Computation of ZMax.
44for i from 0 to MAXINDEX-1 do
45    t_x := center[i] + 2^(-L-1) :
46    zmax[i] := (t_x*r[i]-1) :
47    t_x := center[i] - 2^(-L-1) :
48    zmin[i] := (t_x*r[i]-1) :
49    zabsmax[i] := max(abs(zmin[i]), abs(zmax[i])):
50od:
51for i from MAXINDEX to 2^L do
52    t_x := center[i] + 2^(-L-2) :
53    zmax[i] := (t_x*r[i]-1) :
54    t_x := center[i] - 2^(-L-2) :
55    zmin[i] := (t_x*r[i]-1) :
56    zabsmax[i] := max(abs(zmin[i]), abs(zmax[i])):
57od:
58
59zmaxmax:=0:
60zminmin:=0:
61for i from 0 to 2^L do
62    if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi:
63    if zmin[i] < zminmin then zminmin := zmin[i]: fi:
64od:
65printf("zminmin = -2^(%2f)   zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ) :
66
67
68PolyDegreeQuick:=7:
69
70#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof
71
72if(1+1=3) then # The polynomial used to be computed here,
73polyQuick:= polyExact2Ext(x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'), 0):
74
75epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax):
76deltaApproxQuick := numapprox[infnorm]( polyQuick-log(1+x), x=zminmin..zmaxmax):
77else
78# magical polynomial given by Sylvain's tool
79polyQuick := x * (1 + (x * ((-0.50000000000000000000000000000000000000000000000000000000000000e0) + (x * (0.33333333333333342585191871876304503530263900756835937500000000e0 + (x * ((-0.24999999998423708125194764306797878816723823547363281250000000e0) + (x * (0.19999999996748479835773082413652446120977401733398437500000000e0 + (x * ((-0.16666957260954737285452154083031928166747093200683593750000000e0) + (x * 0.14286056338555042088955815415829420089721679687500000000000000e0)))))))))))):
80
81# validated infinite norms given by Christoph's tool
82epsilonApproxQuick := 2^(-64.119667587221):
83deltaApproxQuick := 2^(-72.230387592107):
84fi:
85printf(" rel approximation error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
86printf(" abs approximation error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :
87
88# For the Paterson/Stockmayer method
89polyQuick0:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'):
90
91## Création des nouveaux coefficients
92polyQuick0:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-1,0], 1 ,  'deltaApprox'):
93
94## Création des nouveaux coefficients
95a7 := convert(coeff(polyQuick0,x,7),rational):
96a6 := convert(coeff(polyQuick0,x,6)/a7,rational):
97a5 := convert(coeff(polyQuick0,x,5)/a7,rational):
98a4 := convert(coeff(polyQuick0,x,4)/a7,rational):
99
100a3 := convert(coeff(polyQuick0,x,3),rational):
101a2 := convert(coeff(polyQuick0,x,2),rational):
102a1 := convert(coeff(polyQuick0,x,1),rational):
103a0 := convert(coeff(polyQuick0,x,0),rational):
104
105alpha0 := a5 - 1:
106alpha1 := a6:
107alpha2 := a4 - alpha0 * a6:
108
109##
110a7_hi,a7_lo         := hiloExt(a7):
111a3_hi,a3_lo         := hiloExt(a3):
112a2_hi,a2_lo         := hiloExt(a2):
113a0_hi,a0_lo         := hiloExt(a0):
114
115c_hi,c_lo         := hiloExt(alpha0):
116alpha_hi,alpha_lo := hiloExt(alpha1):
117beta_hi,beta_lo   := hiloExt(alpha2):
118
119#coef_c := c_hi + c_lo:
120#coef_alpha := alpha_hi + alpha_lo:
121#coef_beta := beta_hi + beta_lo:
122#coef_c7 := a7_hi + a7_lo:
123#coef_c3 := a3_hi + a3_lo:
124#coef_c2 := a2_hi + a2_lo:
125#coef_c0 := a0_hi + a0_lo:
126
127coef_c := c_hi:
128coef_alpha := alpha_hi:
129coef_beta := beta_hi:
130coef_c7 := a7_hi:
131coef_c3 := a3_hi:
132coef_c2 := a2_hi:
133coef_c0 := a0_hi:
134
135polyQuickPat := ((x^2 + coef_c)*(x + coef_alpha) + (x + coef_beta)) * (coef_c7*x^4) + ((coef_c3 * x + coef_c2)*x^2 + (x)):
136
137epsilonApproxQuickPat := numapprox[infnorm]( 1-polyQuickPat/log(1+x), x=zminmin..zmaxmax):
138printf("   approximation rel error for Paterson/Stockmayer in the quick phase is 2^(%2f)\n", log2(epsilonApproxQuickPat) ) :
139deltaApproxQuickPat := numapprox[infnorm]( polyQuickPat-log(1+x), x=zminmin..zmaxmax):
140printf("   approximation abs error for Paterson/Stockmayer in the quick phase is 2^(%2f)\n", log2(deltaApproxQuickPat) ) :
141
142
143
144PolyDegreeAccurate:=14:
145
146printf("Computing the polynomial for accurate phase (this may take some time...)\n"):
147pe:= x  * numapprox[minimax](  log(1+x)/x,  x=-zmaxmax..zmaxmax,  [PolyDegreeAccurate-1,0], 1 ,  'deltaApprox'):
148
149
150MaxDegreeDDE:=8:  #
151
152polyAccurate := polyExact2Ext(pe, MaxDegreeDDE):
153deltaApproxAccurate := numapprox[infnorm](polyAccurate-log(1+x), x=-zmaxmax..zmaxmax):
154epsilonApproxAccurate := numapprox[infnorm]( 1-polyAccurate/log(1+x), x=-zmaxmax..zmaxmax):
155printf("   approximation error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :
156
157
158
159
160filename:="TEMPLOG/log-de.h":
161fd:=fopen(filename, WRITE, TEXT):
162
163fprintf(fd, "/*File generated by maple/log-de.mpl*/\n\n"):
164# some constants
165fprintf(fd, "#define L        %d\n", L):
166fprintf(fd, "#define MAXINDEX %d\n", MAXINDEX):
167fprintf(fd, "#define INDEXMASK %d\n", 2^L-1):
168fprintf(fd, "static const double two64 = %1.25e ;\n", evalf(2^64)):
169
170fprintf(fd, "#if 1\n#define ESTRIN\n#else\n#define PATERSON\n#endif\n"):
171
172fprintf(fd, "#if defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)\n\n"):
173
174fprintf(fd, "#ifdef PATERSON\n"):
175  #  polynomial for quick phase
176  fprintf(fd,"static const long double c7 = %1.50e;\n",coef_c7):
177  fprintf(fd,"static const long double c3 = %1.50e;\n",coef_c3):
178  fprintf(fd,"static const long double c2 = %1.50e;\n",coef_c2):
179  fprintf(fd,"static const long double ps_alpha = %1.50e;\n",coef_alpha):
180  fprintf(fd,"static const long double ps_beta = %1.50e;\n",coef_beta):
181  fprintf(fd,"static const long double ps_c = %1.50e;\n",coef_c):
182fprintf(fd, "#endif/* PATERSON*/\n\n"):
183
184fprintf(fd, "#ifdef ESTRIN\n"):
185 fprintf(fd, "  /* Coefficients are read directly from the array thanks to the following macros */ \n"):
186  for i from PolyDegreeQuick to 2 by -1 do
187    fprintf(fd, "#define c%d  c[%d]\n", i, PolyDegreeQuick-i):
188  od:
189   fprintf(fd, "static const double c[%d] =  {\n",PolyDegreeQuick):
190   for i from PolyDegreeQuick to 2 by -1 do
191     fprintf(fd, "   /* c%d  = */ %1.50eL, \n", i, coeff(polyQuick,x,i)):
192   od:
193  fprintf(fd, "}; \n \n"):
194fprintf(fd, "#endif/* ESTRIN */\n\n"):
195
196
197
198  for i from PolyDegreeAccurate to 1 by -1 do
199    fprintf(fd, "#define c%dh  ch[%d]\n", i, PolyDegreeAccurate-i):
200  od:
201
202  for i from MaxDegreeDDE-1 to 1 by -1 do
203    fprintf(fd, "#define c%dl  cl[%d]\n", i, MaxDegreeDDE-1-i):
204  od:
205  fprintf(fd, "#define PREFETCH_POLY_ACCURATE \n"):
206  fprintf(fd, "\n#else /* not(defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64)),\n   assuming Itanium, otherwise we shouldn't be there */ \n\n"):
207  fprintf(fd, "#define PREFETCH_POLY_QUICK "):
208  for i from PolyDegreeQuick to 2 by -1 do
209    fprintf(fd, "c%d=c[%d]; ", i, PolyDegreeQuick-i):
210  od:
211  fprintf(fd, "\n#define PREFETCH_POLY_ACCURATE "):
212  for i from PolyDegreeAccurate to 1 by -1 do
213    fprintf(fd, "c%dh=ch[%d]; ", i, PolyDegreeAccurate-i):
214    if i mod 4 =0 then  fprintf(fd, "\\\n        "): fi:
215  od:
216  fprintf(fd, "\\\n        "):
217  for i from MaxDegreeDDE-1 to 1 by -1 do
218    fprintf(fd, "c%dl=cl[%d]; ", i, MaxDegreeDDE-1-i):
219  od:
220
221  fprintf(fd, "\n#endif /* defined(CRLIBM_TYPECPU_X86) || defined(CRLIBM_TYPECPU_AMD64) */ \n\n"):
222
223  # Various constants
224  fprintf(fd, "static const long double log2H  = %1.50eL ;\n", log2h):
225  fprintf(fd, "static const long double log2L  = %1.50eL ;\n", log2l):
226
227  # The polynomials
228  #  polynomial for quick phase
229#  for i from PolyDegreeQuick to 1 by -1 do
230#    fprintf(fd, "static const long double c%d =    %1.50eL ;\n", i, coeff(polyQuick,x,i)):
231#  od:
232
233  #  polynomial for accurate phase
234  fprintf(fd, "static const long double ch[%d] =  {\n",PolyDegreeAccurate):
235   for i from PolyDegreeAccurate to 1 by -1 do
236     (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
237     fprintf(fd, "   /* ch%d  = */ %1.50eL, \n", i, ch):
238   od:
239  fprintf(fd, "}; \n \n"):
240
241  fprintf(fd, "static const long double cl[%d] =  {\n", MaxDegreeDDE):
242  for i from MaxDegreeDDE-1 to 1 by -1 do
243     (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
244     fprintf(fd, "   /* cl%d  = */ %1.50eL, \n", i, cl):
245   od:
246  fprintf(fd, "}; \n \n"):
247
248
249  # The tables
250
251
252  fprintf(fd, "typedef struct rri_tag {float r; long double logirh;  long double logirl; } rri ;  \n"):
253  fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L):
254  for i from 0 to 2^L-1 do
255      fprintf(fd, "  { \n"):
256      fprintf(fd, "    %1.50eL,   /* r[%d] */ \n", r[i], i):
257      fprintf(fd, "    %1.50eL, /* logirh[%d] */ \n", logirh[i], i):
258      fprintf(fd, "    %1.50eL, /* logirl[%d] */ \n", logirl[i], i):
259      fprintf(fd, "  } "):
260      if(i<2^L-1) then  fprintf(fd, ", \n"): fi
261  od:
262fprintf(fd, "}; \n \n"):
263
264fclose(fd):
265
266printf("\n************ DONE TEMPLOG/log-de.h ************\n");
267
268
269
270
271filename:="TEMPLOG/polynomials.sed":
272fd:=fopen(filename, WRITE, TEXT):
273for i from PolyDegreeAccurate to 1 by -1 do
274    (ch, cl) := hiloExt(coeff(polyAccurate,x,i)):
275    fprintf(fd, " s/_c%dh/%1.40e/g\n", i, ch):
276    if(i<MaxDegreeDDE) then
277        fprintf(fd, " s/_c%dl/%1.40e/g\n", i, cl):
278    fi:
279od:
280for i from PolyDegreeQuick to 1 by -1 do
281    fprintf(fd, " s/_c%d/%1.40e/g\n", i, coeff(polyQuick,x,i)):
282od:
283fprintf(fd, " s/_deltaApproxQuick/%1.40e/g\n", deltaApproxQuick):
284fprintf(fd, " s/_epsilonApproxQuick/%1.40e/g\n", epsilonApproxQuick):
285fprintf(fd, " s/_deltaApproxAccurate/%1.40e/g\n", deltaApproxAccurate):
286fprintf(fd, " s/_epsilonApproxAccurate/%1.40e/g\n", epsilonApproxAccurate):
287fprintf(fd, " s/_log2h/%1.40e/g\n", log2h):
288fprintf(fd, " s/_log2l/%1.40e/g\n", log2l):
289fclose(fd):
290
291printf("\n************ DONE TEMPLOG/polynomials.sed ************\n");
292
293for j from 0 to 2^L-1 do
294    filename:=cat("TEMPLOG/log-de_",j,".sed"):
295    fd:=fopen(filename, WRITE, TEXT):
296    fprintf(fd, " s/_rval/%1.40e/g\n",     r[j]):
297    fprintf(fd, " s/_logirh/%1.40e/g\n", logirh[j]):
298    fprintf(fd, " s/_logirl/%1.40e/g\n", logirl[j]):
299    fprintf(fd, " s/_zabsmax/%1.40e/g\n", zabsmax[j]):
300    fclose(fd):
301  od:
302
303
304
305printf("\n************ DONE TEMPLOG/log-de*.sed ************\n");
306
307# A shell script to use them
308filename:="../gappa/log-de/run-log-de-proof.sh":
309fd:=fopen(filename, WRITE, TEXT):
310fprintf(fd, "#!/bin/sh\n"):
311fprintf(fd, "# You probably need to edit the path to the gappa executable\n"):
312fprintf(fd, "GAPPA=~/local/src/gappa/src/gappa\n"):
313fprintf(fd, "CRLIBMDIR=../..\n"):
314
315fprintf(fd, "echo Accurate phase, case E=0, index=0: 1>&2\n"):
316fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index0-E0.gappa | $GAPPA \n"):
317
318fprintf(fd, "echo Accurate phase, case E!=0, index=0 1>&2\n"):
319fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index0-E1N.gappa | $GAPPA \n"):
320
321fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
322fprintf(fd, "  echo Accurate phase, case E=0, index=$num 1>&2\n"):
323fprintf(fd, "  sed -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index1N-E0.gappa | $GAPPA \n"):
324fprintf(fd, "  echo 1>&2\n"):
325fprintf(fd, "done\n"):
326
327fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
328fprintf(fd, "  echo Accurate phase, case E!=0, index = $num 1>&2 \n"):
329fprintf(fd, "  sed -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-acc-index1N-E1N.gappa | $GAPPA \n"):
330fprintf(fd, "  echo 1>&2\n"):
331fprintf(fd, "done\n\n"):
332
333
334fprintf(fd, "echo Quick phase, case E=0, index=0  1>&2\n"):
335fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-index0-E0.gappa | $GAPPA \n"):
336fprintf(fd, "  echo 1>&2\n"):
337
338fprintf(fd, "echo Quick phase, case E!=0, index=0 1>&2 \n"):
339fprintf(fd, "sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_0.sed $CRLIBMDIR/gappa/log-de/log-de-index0-E1N.gappa | $GAPPA \n"):
340fprintf(fd, "  echo 1>&2\n"):
341
342
343fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
344fprintf(fd, "  echo Quick phase, for all E,  index=$num 1>&2\n"):
345fprintf(fd, "  sed  -f $CRLIBMDIR/maple/TEMPLOG/polynomials.sed  -f $CRLIBMDIR/maple/TEMPLOG/log-de_$num.sed $CRLIBMDIR/gappa/log-de/log-de-index1N-E0N.gappa | $GAPPA \n"):
346fprintf(fd, "  echo 1>&2\n"):
347fprintf(fd, "done\n"):
348
349fclose(fd):
350
351printf("\n************ DONE TEMPLOG/run-log-de-proof.sh ************\n"):
352printf("Now you should run \n"):
353printf(" sh ../../gappa/log-de/run-log-de-proof.sh  2> ../../maple/TEMPLOG/Gappa.out\n"):
354printf("  (You probably need to edit the path to the gappa executable within run-log-de-proof.sh)\n"):
355printf("Then look at TEMPLOG/Gappa.out. It shouldn't contain 'some enclosures were not satisfied'.\n If it does, report it !\n"):
356
357
358printf("----DONE---\n") :
359
360