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 "exp-td.mpl";
26Digits := 120:
27
28interface(quiet=true):
29
30read "common-procedures.mpl":
31read "triple-double.mpl":
32mkdir("TEMPEXP"):
33
34L := 12:
35
36printf("   memory requirement for L = %d and two triple-double tables: %d bytes\n",L,48*2^(ceil(L/2)));
37
38rmax := log(2) / (2^(L+1)):
39
40printf("   maximal absolute value for rmax = 2^(%f)\n",log[2](rmax)):
41
42
43MsLog2Div2L := evalf(-log(2)/(2^L)):
44
45msLog2Div2Lh, msLog2Div2Lm, msLog2Div2Ll := hi_mi_lo(MsLog2Div2L):
46
47epsMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm + msLog2Div2Ll) - MsLog2Div2L)/MsLog2Div2L)):
48epsDDMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm) - MsLog2Div2L)/MsLog2Div2L)):
49
50printf("   error made by storing MsLog2Div2L as a double-double: 2^(%f)\n",log[2](epsDDMsLog2Div2L)):
51printf("   error made by storing MsLog2Div2L as a triple-double: 2^(%f)\n",log[2](epsMsLog2Div2L)):
52
53gap := -floor(-log[2](abs(msLog2Div2Lm/msLog2Div2Lh))):
54
55printf("   |msLog2Div2Lm| <= 2^(%f) * |msLog2Div2Lh|\n",gap):
56
57
58log2InvMult2L := nearest(2^L / (log(2))):
59
60shiftConst := 2^(52) + 2^(51):
61
62indexmask1 := 2^(L/2) - 1:
63indexmask2 := indexmask1 * 2^(L/2):
64
65largest := 2^(1023) * ((2^(53) - 1) / 2^(52)):
66smallest := 2^(-1023) * 1 * 2^(-51):
67
68overflowbound := nearest(log(largest)):
69
70overflowboundHex := ieeehexa(overflowbound):
71overflowSimplebound := convert(overflowboundHex[1],decimal,hex):
72
73underflowbound := nearest(log(2^(-1075))):
74
75denormbound := nearest(log(2^(-1022) * 1)):
76
77
78overUnderflowboundHex := ieeehexa(min(abs(underflowbound),min(abs(overflowbound),abs(denormbound)))):
79overUnderflowSimplebound := convert(overUnderflowboundHex[1],decimal,hex):
80
81twoPowerM1000 := 2^(-1000):
82twoPower1000 := 2^(1000):
83
84twoM52 := 2^(-52):
85mTwoM53 := - 2^(-53):
86
87for i from 0 to 2^(L/2) - 1 do
88	twoPowerIndex1hi[i], twoPowerIndex1mi[i], twoPowerIndex1lo[i] := hi_mi_lo(evalf(2^(i/(2^L)))):
89	twoPowerIndex2hi[i], twoPowerIndex2mi[i], twoPowerIndex2lo[i] := hi_mi_lo(evalf(2^(i/(2^(L/2))))):
90od:
91
92
93
94PolyDegreeQuick:=4:
95printf("   degree of the polynomial used in the quick phase is %d\n",PolyDegreeQuick);
96
97
98
99polyQuick:= poly_exact(1 + x + 0.5*x^2 + x^3 * (numapprox[minimax](((exp(x) - (1 + x + 0.5*x^2))/x^3),
100				x=-rmax..rmax, [PolyDegreeQuick-3,0], 1 ,  'deltaApprox'))):
101
102epsilonApproxQuick := numapprox[infnorm]( ((polyQuick-1)/(exp(x)-1))-1, x=-rmax..rmax):
103printf("   approximation rel error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
104deltaApproxQuick := numapprox[infnorm]( polyQuick-exp(x), x=-rmax..rmax):
105printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :
106
107
108
109PolyDegreeAccurate:=7:
110
111printf("   degree of the polynomial used in the accurate phase is %d\n",PolyDegreeAccurate):
112
113DDNumberAccu:=5:
114
115
116printf("   number of double doubles used for the coefficients is %d\n",DDNumberAccu):
117
118
119
120polyAccurate:= poly_exact2(1 + x + 0.5*x^2 + x^3 * (numapprox[minimax](((exp(x) - (1 + x + 0.5*x^2))/x^3),
121				x=-rmax..rmax, [PolyDegreeAccurate-3,0], 1 ,  'deltaApprox')),
122				DDNumberAccu):
123
124
125epsilonApproxAccurate := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-rmax..rmax):
126printf("   approximation rel error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :
127deltaApproxAccurate := numapprox[infnorm]( polyAccurate-exp(x), x=-rmax..rmax):
128printf("   approximation abs error for the accurate phase is 2^(%2f)\n", log2(deltaApproxAccurate) ) :
129
130
131epsilonApproxRmAccurate := numapprox[infnorm]( (x/(exp(x)-1))-1, x=-rmax*2^(-52)..rmax*2^(-52)):
132epsilonApproxRlAccurate := numapprox[infnorm]( (x/(exp(x)-1))-1, x=-rmax*2^(-105)..rmax*2^(-105)):
133
134printf("   approximation rel error for approximating exp(rm) - 1 by rm is 2^(%2f)\n", log2(abs(epsilonApproxRmAccurate))):
135printf("   approximation rel error for approximating exp(rl) - 1 by rl is 2^(%2f)\n", log2(abs(epsilonApproxRlAccurate))):
136
137epsilonApproxAccurateSpecial := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-2^(-30)..2^(-30)):
138printf("   approximation rel error for the accurate phase in the special interval (|r| \\leq 2^(-30)) is 2^(%2f)\n",
139log2(epsilonApproxAccurateSpecial) ) :
140
141epsilonApproxAccurateSpecial2 := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-2^(-18)..2^(-18)):
142printf("   approximation rel error for the accurate phase in the special interval (|r| \\leq 2^(-18)) is 2^(%2f)\n",
143log2(epsilonApproxAccurateSpecial2) ) :
144
145
146
147epsilon_quick := 2^(-64): # The Gappa proof will show this bound
148
149
150
151
152
153#-------------------------------------------------------------------
154# Output
155
156filename:="TEMPEXP/exp-td.h":
157fd:=fopen(filename, WRITE, TEXT):
158
159fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):
160
161fprintf(fd, "\n/*File generated by maple/exp-td.mpl*/\n"):
162
163fprintf(fd, "\#define L %d\n",L):
164fprintf(fd, "\#define LHALF %d\n",L/2):
165fprintf(fd, "\#define log2InvMult2L %1.50e\n",log2InvMult2L):
166fprintf(fd, "\#define msLog2Div2Lh %1.50e\n",msLog2Div2Lh):
167fprintf(fd, "\#define msLog2Div2Lm %1.50e\n",msLog2Div2Lm):
168fprintf(fd, "\#define msLog2Div2Ll %1.50e\n",msLog2Div2Ll):
169fprintf(fd, "\#define shiftConst %1.50e\n",shiftConst):
170fprintf(fd, "\#define INDEXMASK1 0x%08x\n",indexmask1):
171fprintf(fd, "\#define INDEXMASK2 0x%08x\n",indexmask2):
172fprintf(fd, "\#define OVRUDRFLWSMPLBOUND 0x%08x\n",overUnderflowSimplebound):
173fprintf(fd, "\#define OVRFLWBOUND %1.50e\n",overflowbound):
174fprintf(fd, "\#define LARGEST %1.50e\n",largest):
175fprintf(fd, "\#define SMALLEST %1.50e\n",smallest):
176fprintf(fd, "\#define DENORMBOUND %1.50e\n",denormbound):
177fprintf(fd, "\#define UNDERFLWBOUND %1.50e\n",underflowbound):
178fprintf(fd, "\#define twoPowerM1000 %1.50e\n",twoPowerM1000):
179fprintf(fd, "\#define twoPower1000 %1.50e\n",twoPower1000):
180fprintf(fd, "\#define ROUNDCST %1.50e\n", compute_rn_constant(epsilon_quick)):
181fprintf(fd, "\#define RDROUNDCST %1.50e\n", epsilon_quick):
182fprintf(fd, "\#define twoM52 %1.50e\n", twoM52):
183fprintf(fd, "\#define mTwoM53 %1.50e\n", mTwoM53):
184
185fprintf(fd,"\n\n"):
186
187for i from 3 to PolyDegreeQuick do
188	fprintf(fd, "\#define c%d %1.50e\n",i,coeff(polyQuick,x,i)):
189od:
190
191
192for i from 3 to DDNumberAccu-1 do
193	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
194	fprintf(fd, "\#define accPolyC%dh %1.50e\n",i,hi):
195	fprintf(fd, "\#define accPolyC%dl %1.50e\n",i,lo):
196od:
197
198for i from DDNumberAccu to PolyDegreeAccurate do
199	fprintf(fd, "\#define accPolyC%d %1.50e\n",i,coeff(polyAccurate,x,i)):
200od:
201
202fprintf(fd,"\n\n"):
203
204# Print the tables
205fprintf(fd, "typedef struct tPi_t_tag {double hi; double mi; double lo;} tPi_t;  \n"):
206fprintf(fd, "static const tPi_t twoPowerIndex1[%d] = {\n", 2^(L/2)):
207for i from 0 to 2^(L/2)-1 do
208      fprintf(fd, "  { \n"):
209      fprintf(fd, "    %1.50e, /* twoPowerIndex1hi[%d] */ \n", twoPowerIndex1hi[i], i):
210      fprintf(fd, "    %1.50e, /* twoPowerIndex1mi[%d] */ \n", twoPowerIndex1mi[i], i):
211      fprintf(fd, "    %1.50e, /* twoPowerIndex1lo[%d] */ \n", twoPowerIndex1lo[i], i):
212      fprintf(fd, "  } "):
213      if(i<2^(L/2)-1) then  fprintf(fd, ", \n"): fi
214od:
215fprintf(fd, "}; \n \n"):
216fprintf(fd, "static const tPi_t twoPowerIndex2[%d] = {\n", 2^(L/2)):
217for i from 0 to 2^(L/2)-1 do
218      fprintf(fd, "  { \n"):
219      fprintf(fd, "    %1.50e, /* twoPowerIndex2hi[%d] */ \n", twoPowerIndex2hi[i], i):
220      fprintf(fd, "    %1.50e, /* twoPowerIndex2mi[%d] */ \n", twoPowerIndex2mi[i], i):
221      fprintf(fd, "    %1.50e, /* twoPowerIndex2lo[%d] */ \n", twoPowerIndex2lo[i], i):
222      fprintf(fd, "  } "):
223      if(i<2^(L/2)-1) then  fprintf(fd, ", \n"): fi
224od:
225fprintf(fd, "}; \n \n"):
226
227fprintf(fd, "\n\n"):
228
229fclose(fd):
230
231filename:="TEMPEXP/exp-td-accurate.sed":
232fd:=fopen(filename, WRITE, TEXT):
233fprintf(fd, " s/_rhmax/%1.50e/g\n", rmax):
234fprintf(fd, " s/_rmmax/%1.50e/g\n", rmax*2^(-52)):
235fprintf(fd, " s/_rlmax/%1.50e/g\n", rmax*2^(-105)):
236fprintf(fd, " s/_epsilonApproxAccurate/%1.50e/g\n", epsilonApproxAccurate):
237fprintf(fd, " s/_epsilonApproxRmAccurate/%1.50e/g\n", epsilonApproxRmAccurate):
238fprintf(fd, " s/_epsilonApproxRlAccurate/%1.50e/g\n", epsilonApproxRlAccurate):
239fprintf(fd, " s/_epsilonApproxSpecial/%1.50e/g\n", epsilonApproxAccurateSpecial):
240
241for i from 3 to DDNumberAccu-1 do
242	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
243	fprintf(fd, "s/_accPolyC%dh/%1.50e/g\n",i,hi):
244	fprintf(fd, "s/_accPolyC%dl/%1.50e/g\n",i,lo):
245od:
246
247for i from DDNumberAccu to PolyDegreeAccurate do
248	fprintf(fd, "s/_accPolyC%d/%1.50e/g\n",i,coeff(polyAccurate,x,i)):
249od:
250
251fclose(fd):
252
253
254
255printf("----DONE---\n") :
256
257
258