xref: /aosp_15_r20/external/fdlibm/readme (revision 1e651e1ef2b613db2c4b29ae59c1de74cf0222ae)
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