986
|
1 |
|
|
2 |
*********************************
|
|
3 |
* Announcing FDLIBM Version 5.3 *
|
|
4 |
*********************************
|
|
5 |
============================================================
|
|
6 |
FDLIBM
|
|
7 |
============================================================
|
|
8 |
developed at Sun Microsystems, Inc.
|
|
9 |
|
|
10 |
What's new in FDLIBM 5.3?
|
|
11 |
|
|
12 |
CONFIGURE
|
|
13 |
To build FDLIBM, edit the supplied Makefile or create
|
|
14 |
a local Makefile by running "sh configure"
|
|
15 |
using the supplied configure script contributed by Nelson Beebe
|
|
16 |
|
|
17 |
BUGS FIXED
|
|
18 |
|
|
19 |
1. e_pow.c incorrect results when
|
|
20 |
x is very close to -1.0 and y is very large, e.g.
|
|
21 |
pow(-1.0000000000000002e+00,4.5035996273704970e+15) = 0
|
|
22 |
pow(-9.9999999999999978e-01,4.5035996273704970e+15) = 0
|
|
23 |
Correct results are close to -e and -1/e.
|
|
24 |
|
|
25 |
2. k_tan.c error was > 1 ulp target for FDLIBM
|
|
26 |
5.2: Worst error at least 1.45 ulp at
|
|
27 |
tan(1.7765241907548024E+269) = 1.7733884462610958E+16
|
|
28 |
5.3: Worst error 0.96 ulp
|
|
29 |
|
|
30 |
NOT FIXED YET
|
|
31 |
|
|
32 |
3. Compiler failure on non-standard code
|
|
33 |
Statements like
|
|
34 |
*(1+(int*)&t1) = 0;
|
|
35 |
are not standard C and cause some optimizing compilers (e.g. GCC)
|
|
36 |
to generate bad code under optimization. These cases
|
|
37 |
are to be addressed in the next release.
|
|
38 |
|
|
39 |
FDLIBM (Freely Distributable LIBM) is a C math library
|
|
40 |
for machines that support IEEE 754 floating-point arithmetic.
|
|
41 |
In this release, only double precision is supported.
|
|
42 |
|
|
43 |
FDLIBM is intended to provide a reasonably portable (see
|
|
44 |
assumptions below), reference quality (below one ulp for
|
|
45 |
major functions like sin,cos,exp,log) math library
|
|
46 |
(libm.a). For a copy of FDLIBM, please see
|
|
47 |
http://www.netlib.org/fdlibm/
|
|
48 |
or
|
|
49 |
http://www.validlab.com/software/
|
|
50 |
|
|
51 |
--------------
|
|
52 |
1. ASSUMPTIONS
|
|
53 |
--------------
|
|
54 |
FDLIBM (double precision version) assumes:
|
|
55 |
a. IEEE 754 style (if not precise compliance) arithmetic;
|
|
56 |
b. 32 bit 2's complement integer arithmetic;
|
|
57 |
c. Each double precision floating-point number must be in IEEE 754
|
|
58 |
double format, and that each number can be retrieved as two 32-bit
|
|
59 |
integers through the using of pointer bashing as in the example
|
|
60 |
below:
|
|
61 |
|
|
62 |
Example: let y = 2.0
|
|
63 |
double fp number y: 2.0
|
|
64 |
IEEE double format: 0x4000000000000000
|
|
65 |
|
|
66 |
Referencing y as two integers:
|
|
67 |
*(int*)&y,*(1+(int*)&y) = {0x40000000,0x0} (on sparc)
|
|
68 |
{0x0,0x40000000} (on 386)
|
|
69 |
|
|
70 |
Note: Four macros are defined in fdlibm.h to handle this kind of
|
|
71 |
retrieving:
|
|
72 |
|
|
73 |
__HI(x) the high part of a double x
|
|
74 |
(sign,exponent,the first 21 significant bits)
|
|
75 |
__LO(x) the least 32 significant bits of x
|
|
76 |
__HIp(x) same as __HI except that the argument is a pointer
|
|
77 |
to a double
|
|
78 |
__LOp(x) same as __LO except that the argument is a pointer
|
|
79 |
to a double
|
|
80 |
|
|
81 |
To ensure obtaining correct ordering, one must define __LITTLE_ENDIAN
|
|
82 |
during compilation for little endian machine (like 386,486). The
|
|
83 |
default is big endian.
|
|
84 |
|
|
85 |
If the behavior of pointer bashing is undefined, one may hack on the
|
|
86 |
macro in fdlibm.h.
|
|
87 |
|
|
88 |
d. IEEE exceptions may trigger "signals" as is common in Unix
|
|
89 |
implementations.
|
|
90 |
|
|
91 |
-------------------
|
|
92 |
2. EXCEPTION CASES
|
|
93 |
-------------------
|
|
94 |
All exception cases in the FDLIBM functions will be mapped
|
|
95 |
to one of the following four exceptions:
|
|
96 |
|
|
97 |
+-huge*huge, +-tiny*tiny, +-1.0/0.0, +-0.0/0.0
|
|
98 |
(overflow) (underflow) (divided-by-zero) (invalid)
|
|
99 |
|
|
100 |
For example, log(0) is a singularity and is thus mapped to
|
|
101 |
-1.0/0.0 = -infinity.
|
|
102 |
That is, FDLIBM's log will compute -one/zero and return the
|
|
103 |
computed value. On an IEEE machine, this will trigger the
|
|
104 |
divided-by-zero exception and a negative infinity is returned by
|
|
105 |
default.
|
|
106 |
|
|
107 |
Similarly, exp(-huge) will be mapped to tiny*tiny to generate
|
|
108 |
an underflow signal.
|
|
109 |
|
|
110 |
|
|
111 |
--------------------------------
|
|
112 |
3. STANDARD CONFORMANCE WRAPPER
|
|
113 |
--------------------------------
|
|
114 |
The default FDLIBM functions (compiled with -D_IEEE_LIBM flag)
|
|
115 |
are in "IEEE spirit" (i.e., return the most reasonable result in
|
|
116 |
floating-point arithmetic). If one wants FDLIBM to comply with
|
|
117 |
standards like SVID, X/OPEN, or POSIX/ANSI, then one can
|
|
118 |
create a multi-standard compliant FDLIBM. In this case, each
|
|
119 |
function in FDLIBM is actually a standard compliant wrapper
|
|
120 |
function.
|
|
121 |
|
|
122 |
File organization:
|
|
123 |
1. For FDLIBM's kernel (internal) function,
|
|
124 |
File name Entry point
|
|
125 |
---------------------------
|
|
126 |
k_sin.c __kernel_sin
|
|
127 |
k_tan.c __kernel_tan
|
|
128 |
---------------------------
|
|
129 |
2. For functions that have no standards conflict
|
|
130 |
File name Entry point
|
|
131 |
---------------------------
|
|
132 |
s_sin.c sin
|
|
133 |
s_erf.c erf
|
|
134 |
---------------------------
|
|
135 |
3. Ieee754 core functions
|
|
136 |
File name Entry point
|
|
137 |
---------------------------
|
|
138 |
e_exp.c __ieee754_exp
|
|
139 |
e_sinh.c __ieee754_sinh
|
|
140 |
---------------------------
|
|
141 |
4. Wrapper functions
|
|
142 |
File name Entry point
|
|
143 |
---------------------------
|
|
144 |
w_exp.c exp
|
|
145 |
w_sinh.c sinh
|
|
146 |
---------------------------
|
|
147 |
|
|
148 |
Wrapper functions will twist the result of the ieee754
|
|
149 |
function to comply to the standard specified by the value
|
|
150 |
of _LIB_VERSION
|
|
151 |
if _LIB_VERSION = _IEEE_, return the ieee754 result;
|
|
152 |
if _LIB_VERSION = _SVID_, return SVID result;
|
|
153 |
if _LIB_VERSION = _XOPEN_, return XOPEN result;
|
|
154 |
if _LIB_VERSION = _POSIX_, return POSIX/ANSI result.
|
|
155 |
(These are macros, see fdlibm.h for their definition.)
|
|
156 |
|
|
157 |
|
|
158 |
--------------------------------
|
|
159 |
4. HOW TO CREATE FDLIBM's libm.a
|
|
160 |
--------------------------------
|
|
161 |
There are two types of libm.a. One is IEEE only, and the other is
|
|
162 |
multi-standard compliant (supports IEEE,XOPEN,POSIX/ANSI,SVID).
|
|
163 |
|
|
164 |
To create the IEEE only libm.a, use
|
|
165 |
make "CFLAGS = -D_IEEE_LIBM"
|
|
166 |
This will create an IEEE libm.a, which is smaller in size, and
|
|
167 |
somewhat faster.
|
|
168 |
|
|
169 |
To create a multi-standard compliant libm, use
|
|
170 |
make "CFLAGS = -D_IEEE_MODE" --- multi-standard fdlibm: default
|
|
171 |
to IEEE
|
|
172 |
make "CFLAGS = -D_XOPEN_MODE" --- multi-standard fdlibm: default
|
|
173 |
to X/OPEN
|
|
174 |
make "CFLAGS = -D_POSIX_MODE" --- multi-standard fdlibm: default
|
|
175 |
to POSIX/ANSI
|
|
176 |
make "CFLAGS = -D_SVID3_MODE" --- multi-standard fdlibm: default
|
|
177 |
to SVID
|
|
178 |
|
|
179 |
|
|
180 |
Here is how one makes a SVID compliant libm.
|
|
181 |
Make the library by
|
|
182 |
make "CFLAGS = -D_SVID3_MODE".
|
|
183 |
The libm.a of FDLIBM will be multi-standard compliant and
|
|
184 |
_LIB_VERSION is initialized to the value _SVID_ .
|
|
185 |
|
|
186 |
example1:
|
|
187 |
---------
|
|
188 |
main()
|
|
189 |
{
|
|
190 |
double y0();
|
|
191 |
printf("y0(1e300) = %1.20e\n",y0(1e300));
|
|
192 |
exit(0);
|
|
193 |
}
|
|
194 |
|
|
195 |
% cc example1.c libm.a
|
|
196 |
% a.out
|
|
197 |
y0: TLOSS error
|
|
198 |
y0(1e300) = 0.00000000000000000000e+00
|
|
199 |
|
|
200 |
|
|
201 |
It is possible to change the default standard in multi-standard
|
|
202 |
fdlibm. Here is an example of how to do it:
|
|
203 |
example2:
|
|
204 |
---------
|
|
205 |
#include "fdlibm.h" /* must include FDLIBM's fdlibm.h */
|
|
206 |
main()
|
|
207 |
{
|
|
208 |
double y0();
|
|
209 |
_LIB_VERSION = _IEEE_;
|
|
210 |
printf("IEEE: y0(1e300) = %1.20e\n",y0(1e300));
|
|
211 |
_LIB_VERSION = _XOPEN_;
|
|
212 |
printf("XOPEN y0(1e300) = %1.20e\n",y0(1e300));
|
|
213 |
_LIB_VERSION = _POSIX_;
|
|
214 |
printf("POSIX y0(1e300) = %1.20e\n",y0(1e300));
|
|
215 |
_LIB_VERSION = _SVID_;
|
|
216 |
printf("SVID y0(1e300) = %1.20e\n",y0(1e300));
|
|
217 |
exit(0);
|
|
218 |
}
|
|
219 |
|
|
220 |
% cc example2.c libm.a
|
|
221 |
% a.out
|
|
222 |
IEEE: y0(1e300) = -1.36813604503424810557e-151
|
|
223 |
XOPEN y0(1e300) = 0.00000000000000000000e+00
|
|
224 |
POSIX y0(1e300) = 0.00000000000000000000e+00
|
|
225 |
y0: TLOSS error
|
|
226 |
SVID y0(1e300) = 0.00000000000000000000e+00
|
|
227 |
|
|
228 |
Note: Here _LIB_VERSION is a global variable. If global variables
|
|
229 |
are forbidden, then one should modify fdlibm.h to change
|
|
230 |
_LIB_VERSION to be a global constant. In this case, one
|
|
231 |
may not change the value of _LIB_VERSION as in example2.
|
|
232 |
|
|
233 |
---------------------------
|
|
234 |
5. NOTES ON PORTING FDLIBM
|
|
235 |
---------------------------
|
|
236 |
Care must be taken when installing FDLIBM over existing
|
|
237 |
libm.a.
|
|
238 |
All co-existing function prototypes must agree, otherwise
|
|
239 |
users will encounter mysterious failures.
|
|
240 |
|
|
241 |
So far, the only known likely conflict is the declaration
|
|
242 |
of the IEEE recommended function scalb:
|
|
243 |
|
|
244 |
double scalb(double,double) (1) SVID3 defined
|
|
245 |
double scalb(double,int) (2) IBM,DEC,...
|
|
246 |
|
|
247 |
FDLIBM follows Sun definition and use (1) as default.
|
|
248 |
If one's existing libm.a uses (2), then one may raise
|
|
249 |
the flags _SCALB_INT during the compilation of FDLIBM
|
|
250 |
to get the correct function prototype.
|
|
251 |
(E.g., make "CFLAGS = -D_IEEE_LIBM -D_SCALB_INT".)
|
|
252 |
NOTE that if -D_SCALB_INT is raised, it won't be SVID3
|
|
253 |
conformant.
|
|
254 |
|
|
255 |
--------------
|
|
256 |
6. PROBLEMS ?
|
|
257 |
--------------
|
|
258 |
Please send comments and bug reports to the electronic mail address
|
|
259 |
suggested by:
|
|
260 |
fdlibm-comments AT sun.com
|
|
261 |
|