1#######################################################################
2# This file is part of crlibm, the correctly rounded mathematical library,
3# which has been developed by the Arénaire project at École normale supérieure
4# de Lyon.
5#
6# Copyright (C) 2004-2011 David Defour, Catherine Daramy-Loirat,
7# Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Quirin Lauter,
8# and Jean-Michel Muller
9#
10# This library is free software; you can redistribute it and/or
11# modify it under the terms of the GNU Lesser General Public
12# License as published by the Free Software Foundation; either
13# version 2.1 of the License, or (at your option) any later version.
14#
15# This library is distributed in the hope that it will be useful,
16# but WITHOUT ANY WARRANTY; without even the implied warranty of
17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
18# Lesser General Public License for more details.
19#
20# You should have received a copy of the GNU Lesser General Public
21# License along with this library; if not, write to the Free Software
22# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23
24# To use:
25# restart; read "log2-td.mpl";
26Digits := 120:
27
28interface(quiet=true):
29
30read "common-procedures.mpl":
31read "triple-double.mpl":
32mkdir("TEMPLOG"):
33
34
35# We want log2h + log2m + log2l + delta = log(2) such that
36# log2h and log2m have at least 11 trailing zeros
37# in order to have an exact multiplication with E, which is lower than 1024 in
38# magnitude
39# The resting accuracy is enough for both quick and accurate phases.
40
41log2acc := log(2):
42log2h := round(log2acc * 2**(floor(-log[2](abs(log2acc))) + (53 - 11))) / 2**(floor(-log[2](abs(log2acc))) + (53 - 11)):
43log2m := round((log2acc - log2h) * 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11))) / 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11)):
44log2l := log2acc - (log2h + log2m):
45
46
47L := 7: # number of bits used to address the table
48
49MAXINDEX    := round(2^L * (sqrt(2)-1)):
50
51for i from 0 to MAXINDEX-1 do
52    center[i] := 1 + i*2^(-L): # center[i] in [1, 2[
53    t :=  evalf(1/center[i]):
54    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
55    (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))):
56od:
57for i from MAXINDEX to 2^L do
58    # y has been divided by two, center[i] in [0.5, 1[
59    center[i]:=(1 + i*2^(-L)) / 2:
60    t :=  evalf(1/center[i]):
61    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
62    (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))):
63od:
64
65
66
67
68#Computation of ZMax.
69for i from 0 to MAXINDEX-1 do
70    __x := center[i] + 2^(-L-1) :
71    zmax[i] := (__x*r[i]-1) :
72    __x := center[i] - 2^(-L-1) :
73    zmin[i] := (__x*r[i]-1) :
74od:
75for i from MAXINDEX to 2^L do
76    __x := center[i] + 2^(-L-2) :
77    zmax[i] := (__x*r[i]-1) :
78    __x := center[i] - 2^(-L-2) :
79    zmin[i] := (__x*r[i]-1) :
80od:
81
82zmaxmax:=0:
83zminmin:=0:
84for i from 0 to 2^L do
85    tabulated_value := logih[i] + logim[i] + logil[i]:
86    poly_approx_min := evalf(log(1+zmin[i])):
87    poly_approx_max := evalf(log(1+zmax[i])):
88
89    # Test if we have a case where we cancellate a lot
90    # i.e. the polynomial approximation value is greater than half the tabulated value
91    # the tabulated value is not exactly zero and we are of opposite sign
92
93    if ((abs(poly_approx_min) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_min * tabulated_value < 0)) then
94	printf("Polynomial approximation is greater in magnitude in zmin[%d] than half the tabluated value\n",i):
95
96	printf("The tabulated value is %1.50e\n",tabulated_value):
97	if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi:
98	printf("The value of polynomial in zmin[%d] is %1.50e\n",i,poly_approx_min):
99        if (poly_approx_min <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_min)))) fi:
100
101	summe := poly_approx_min + tabulated_value:
102	printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))):
103
104
105    fi:
106
107    if ((abs(poly_approx_max) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_max * tabulated_value <0)) then
108	printf("Polynomial approximation is greater in magnitude in zmax[%d] than half the tabluated value\n",i):
109
110	printf("The tabulated value is %1.50e\n",tabulated_value):
111	if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi:
112	printf("The value of polynomial in zmax[%d] is %1.50e\n",i,poly_approx_max):
113        if (poly_approx_max <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_max)))) fi:
114
115	summe := poly_approx_max + tabulated_value:
116	printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))):
117
118
119    fi:
120
121
122    if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi:
123    if zmin[i] < zminmin then zminmin := zmin[i]: fi:
124od:
125printf("zminmin = -2^(%2f)   zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ) :
126
127PolyDegreeQuick:=7:
128
129printf("   degree of the polynomial used in the quick phase is %d\n",PolyDegreeQuick);
130
131DDNumberQuick:=3:
132
133printf("   number of double doubles used for the coefficients is %d\n",DDNumberQuick);
134
135#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof
136#and constrain the first two coefficients to 1 and -1/2 in order to save up a full multiplication and a rounding error
137polyQuick:= poly_exact2(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x,
138			x=-zmaxmax..zmaxmax,  [PolyDegreeQuick-3,0], 1 ,  'deltaApprox')))), DDNumberQuick):
139
140#Try to verify the bound for using double double arithmetic.
141#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits)
142#with the maximum value of the first term which is calculated in double precision only
143
144p := unapply(polyQuick,x):
145printf("   using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n",
146	DDNumberQuick,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(DDNumberQuick)),5)):
147
148
149epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax):
150printf("   approximation rel error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
151deltaApproxQuick := numapprox[infnorm]( polyQuick-log(1+x), x=zminmin..zmaxmax):
152printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :
153
154
155PolyDegreeAccurate:=14:
156
157printf("   degree of the polynomial used in the accurate phase is %d\n",PolyDegreeAccurate);
158
159DDNumberAccu:=7:
160TDNumberAccu:=3:
161
162printf("   number of triple doubles used for the coefficients is %d\n",TDNumberAccu);
163printf("   number of double doubles used for the coefficients is %d\n",DDNumberAccu);
164
165
166#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof
167polyAccurate:= poly_exact32(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x,
168				x=-zmaxmax..zmaxmax,  [PolyDegreeAccurate-3,0], 1 ,  'deltaApprox')))),
169				TDNumberAccu, DDNumberAccu):
170
171#Try to verify the bound for using double double arithmetic.
172#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits)
173#with the maximum value of the first term which is calculated in double precision only
174
175pp := unapply(polyAccurate,x):
176printf("   using only %d triple doubles should be fine since the mi z ulp should affect the result starting from bit %f\n",
177	TDNumberAccu,evalf(106 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu)),5)):
178printf("   using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n",
179	DDNumberAccu,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu + DDNumberAccu)),5)):
180
181
182epsilonApproxAccurate := numapprox[infnorm]( 1-polyAccurate/log(1+x), x=zminmin..zmaxmax):
183printf("   approximation rel error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :
184deltaApproxAccurate := numapprox[infnorm]( polyAccurate-log(1+x), x=zminmin..zmaxmax):
185printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxAccurate) ) :
186
187
188
189#Compute now the inverse of ln(2) for the final addition of ln(x) with this constant for obtaining log2(x)
190#Compute also the relative error of the constant stored as a double double.
191Log2inv := evalf(1 / log(2)):
192
193(log2invh, log2invl) := hi_lo(Log2inv):
194
195Log2invhl := log2invh + log2invl:
196
197epsilonLog2invhl := evalf(abs((Log2invhl - Log2inv) / Log2inv)):
198
199printf("   Log2inv = 1 / ln(2) stored as a double-double is exact with a relative error of 2^(%2f)\n",
200	evalf(log[2](epsilonLog2invhl))):
201
202
203#-------------------------------------------------------------------
204# Output
205
206
207filename:="TEMPLOG/log2-td.h":
208fd:=fopen(filename, WRITE, TEXT):
209
210fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):
211
212fprintf(fd, "\n/*File generated by maple/log2-td.mpl*/\n"):
213
214fprintf(fd, "\n\#define L %d\n\n",L):
215fprintf(fd, "\#define MAXINDEX %d\n\n", MAXINDEX):
216fprintf(fd, "\#define INDEXMASK %d\n", 2^L-1):
217fprintf(fd, "\#define two52 %1.50e\n", 2^(52)):
218fprintf(fd, "\#define log2h %1.50e\n", log2h):
219fprintf(fd, "\#define log2m %1.50e\n", log2m):
220fprintf(fd, "\#define log2l %1.50e\n", log2l):
221fprintf(fd, "\#define log2invh %1.50e\n",log2invh):
222fprintf(fd, "\#define log2invl %1.50e\n",log2invl):
223
224epsilon_quick_1 := 2^(-61): # The Gappa proof will show this bound
225epsilon_quick_2 := 2^(-61): # The Gappa proof will show this bound
226fprintf(fd, "\#define ROUNDCST1 %1.50e\n", compute_rn_constant(epsilon_quick_1)):
227fprintf(fd, "\#define ROUNDCST2 %1.50e\n", compute_rn_constant(epsilon_quick_2)):
228fprintf(fd, "\#define RDROUNDCST1 %1.50e\n", epsilon_quick_1):
229fprintf(fd, "\#define RDROUNDCST2 %1.50e\n", epsilon_quick_2):
230
231
232fprintf(fd, "\n\n"):
233
234
235# Print the defines for the define statements
236
237for i from 3 to PolyDegreeQuick do
238	fprintf(fd, "\#define c%d %1.50e\n",i,coeff(polyQuick,x,i)):
239od:
240
241fprintf(fd, "\n\n"):
242
243for i from 3 to (DDNumberAccu + TDNumberAccu -1) do
244	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
245	fprintf(fd, "\#define accPolyC%dh %1.50e\n",i,hi):
246	fprintf(fd, "\#define accPolyC%dl %1.50e\n",i,lo):
247od:
248
249for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do
250	fprintf(fd, "\#define accPolyC%d %1.50e\n",i,coeff(polyAccurate,x,i)):
251od:
252
253fprintf(fd, "\n\n"):
254
255# Print the table
256fprintf(fd, "typedef struct rri_tag {float ri; double logih; double logim; double logil;} rri;  \n"):
257fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L):
258for i from 0 to 2^L-1 do
259      fprintf(fd, "  { \n"):
260      fprintf(fd, "    %1.50e,   /* r[%d] */ \n", r[i], i):
261      fprintf(fd, "    %1.50e, /* logih[%d] */ \n", logih[i], i):
262      fprintf(fd, "    %1.50e, /* logim[%d] */ \n", logim[i], i):
263      fprintf(fd, "    %1.50e, /* logil[%d] */ \n", logil[i], i):
264      fprintf(fd, "  } "):
265      if(i<2^L-1) then  fprintf(fd, ", \n"): fi
266od:
267fprintf(fd, "}; \n \n"):
268
269fclose(fd):
270
271for j from 0 to 2^L-1 do
272    filename:=cat("TEMPLOG/log2-td_",j,".sed"):
273    fd:=fopen(filename, WRITE, TEXT):
274    fprintf(fd, " s/_log2h/%1.50e/g\n", log2h):
275    fprintf(fd, " s/_log2m/%1.50e/g\n", log2m):
276    fprintf(fd, " s/_log2l/%1.50e/g\n", log2l):
277    fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]):
278    fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]):
279    fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]):
280    fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]):
281    fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]):
282    for i from 3 to PolyDegreeQuick do
283        fprintf(fd, " s/_c%d/%1.50e/g\n", i, coeff(polyQuick,x,i)):
284    od:
285    fprintf(fd, " s/_epsilonApproxQuick/%1.50e/g\n", epsilonApproxQuick):
286    fprintf(fd, " s/_epsilonLog2invhl/%1.50e/g\n", epsilonLog2invhl):
287    fclose(fd):
288  od:
289
290for j from 0 to 2^L-1 do
291    filename:=cat("TEMPLOG/log2-td-accurate_",j,".sed"):
292    fd:=fopen(filename, WRITE, TEXT):
293    fprintf(fd, " s/_log2h/%1.50e/g\n", log2h):
294    fprintf(fd, " s/_log2m/%1.50e/g\n", log2m):
295    fprintf(fd, " s/_log2l/%1.50e/g\n", log2l):
296    fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]):
297    fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]):
298    fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]):
299    fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]):
300    fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]):
301    for i from 3 to (DDNumberAccu + TDNumberAccu -1) do
302	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
303        fprintf(fd, " s/_accPolyC%dh/%1.50e/g\n", i, hi):
304        fprintf(fd, " s/_accPolyC%dl/%1.50e/g\n", i, lo):
305    od:
306    for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do
307        fprintf(fd, " s/_accPolyC%d/%1.50e/g\n", i, coeff(polyAccurate,x,i)):
308    od:
309    fprintf(fd, " s/_epsilonApproxAccurate/%1.50e/g\n", epsilonApproxAccurate):
310    fprintf(fd, " s/_epsilonLog2invhl/%1.50e/g\n", epsilonLog2invhl):
311    fclose(fd):
312  od:
313
314# A shell script to use them
315filename:="run-log2-td-proof.sh":
316fd:=fopen(filename, WRITE, TEXT):
317fprintf(fd, "#!/bin/sh\n"):
318fprintf(fd, "# You probably need to edit the path to the gappa executable\n"):
319fprintf(fd, "GAPPA=~/ble/gappa-0.4.5/src/gappa\n"):
320fprintf(fd, "# Test all the possible table value for E=1\n"):
321fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1):
322fprintf(fd, "  echo $num, E=1:\n"):
323fprintf(fd, "  sed -f ./TEMPLOG/log2-td_$num.sed log2-td.gappa | $GAPPA > /dev/null\n"):
324fprintf(fd, "  echo\n"):
325fprintf(fd, "done\n"):
326fprintf(fd, "# For the case E=0 we first handle the cases 0 and %d using log2-td-E0-logir0.gappa\n", 2^L):
327fprintf(fd, "echo 0 and %d, E=0:\n", 2^L):
328fprintf(fd, "sed -f log2-td_0.sed log2-td-E0-logir0.gappa | $GAPPA > /dev/null\n"):
329fprintf(fd, "# then the other cases where logirh <> 0\n"):
330fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
331fprintf(fd, "  echo $num, E=0:\n"):
332fprintf(fd, "  sed -f ./TEMPLOG/log2-td_$num.sed log2-td-E0.gappa | $GAPPA > /dev/null\n"):
333fprintf(fd, "  echo\n"):
334fprintf(fd, "done\n"):
335fprintf(fd, "# Accurate phase: Test all the possible table value for E=1\n"):
336fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1):
337fprintf(fd, "  echo Accurate phase: $num, E=1:\n"):
338fprintf(fd, "  sed -f ./TEMPLOG/log2-td-accurate_$num.sed log2-td-accurate.gappa | $GAPPA > /dev/null\n"):
339fprintf(fd, "  echo\n"):
340fprintf(fd, "done\n"):
341fprintf(fd, "# Accurate phase: For the case E=0 we first handle the cases 0 and %d using log2-td-accurate-E0-logir0.gappa\n", 2^L):
342fprintf(fd, "echo 0 and %d, E=0:\n", 2^L):
343fprintf(fd, "sed -f ./TEMPLOG/log2-td-accurate_0.sed log2-td-accurate-E0-logir0.gappa | $GAPPA > /dev/null\n"):
344fprintf(fd, "# Accurate phase: then the other cases where logirh <> 0\n"):
345fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1):
346fprintf(fd, "  echo $num, E=0:\n"):
347fprintf(fd, "  sed -f ./TEMPLOG/log2-td-accurate_$num.sed log2-td-accurate-E0.gappa | $GAPPA > /dev/null\n"):
348fprintf(fd, "  echo\n"):
349fprintf(fd, "done\n"):
350fclose(fd):
351
352printf("----DONE---\n") :
353
354
355