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