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