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