Beta release of Freely-Distributable libm available for testing
KC Ng
sun!Eng!Kwok.Ng
Fri Feb 12 13:23:16 PST 1993
Subject: Beta release of Freely-Distributable libm available for testing
The SunPro "fdlibm" is now ready for further outside testing. Persons
who are interested in helping to test may obtain the source code.
Although the source code is not legally enchained, consider yourself
morally enjoined not to incorporate it into any products or further
distribute it until external testing is complete. That's why this
announcement is only to numeric-interest.
This library will also be compiled and eventually distributed as the
bundled /usr/ccs/lib/libm.a (what you used to call /usr/lib/libm.a)
with Solaris 2.x, for some x to be determined, for SPARC and x86.
The extra "Sun value added" that was in SunOS 4.x libm.a, beyond that
required for external standards, will continue to be distributed by
SunPro as part of its unbundled compiler products.
The fdlibm source code amounts to some 303K bytes, including an
unconscionable amount devoted to meeting the requirements of various
contradictory standards test suites.
Here's the Readme file in FDLIBM:
**************************************
* Announcing FDLIBM Version 2 (beta) *
**************************************
============================================================
FDLIBM
============================================================
(developed at SunPro, a Sun Microsystems, Inc. business.)
Version 2 (beta), 2/9/93
FDLIBM (Freely Distributable LIBM) is a C math library
for machines that support IEEE 754 floating-point arithmetic.
In this release, only double precision is supported.
FDLIBM is intended to provide a reasonably portable (see
assumption below), reference quality (below one ulp for
major functions like sin,cos,exp,log) math library
(libm.a). For an early release of FDLIBM, please send
e-mail to
kwok.ngaEng.sun.com
--------------
1. ASSUMPTIONS
--------------
FDLIBM assumes:
a. IEEE 754 style (if not precise compliance) arithmetic
b. 32 bit integer arithmetic,
c. each floating-point number must be in IEEE 754 double format,
and that each number can be retrieved as two 32-bit integers;
Example: let y = 2.0
double fp number y: 2.0
IEEE double format: 0x4000000000000000
Referencing y as two integers:
*(int*)&y,*(1+(int*)&y) = {0x40000000,0x0} (on sparc)
{0x0,0x40000000} (on 386)
Note: FDLIBM will detect at run time the correct ordering of
the high and low part of a floating-point number.
d. IEEE exceptions may trigger "signals" as is common in Unix
implementations.
-------------------
2. EXCEPTION CASES
-------------------
All exception cases in the FDLIBM functions will be mapped
to one of the following four exceptions:
+-huge*huge, +-tiny*tiny, +-1.0/0.0, +-0.0/0.0
(overflow) (underflow) (divided-by-zero) (invalid)
For example, log(0) is a singularity and is thus mapped to
-1.0/0.0 = -infinity.
That is, FDLIBM's log will compute -one/zero and return the
computed value. On an IEEE machine, this will trigger the
divided-by-zero exception and a negative infinity is returned by
default.
Similarly, exp(-huge) will be mapped to tiny*tiny to generate
an underflow signal.
--------------------------------
3. STANDARD CONFORMANCE WRAPPER
--------------------------------
The default FDLIBM functions (compiled with -D_IEEE_LIBM flag)
are of IEEE spirit (i.e., return the most reasonable result in
floating-point arithmetic). If one wants FDLIBM to comply
standards like SVID, X/OPEN, or POSIX/ANSI, then one can
create a multi-standard compliant FDLIBM. In this case, each
function in FDLIBM is actually a standard compliant wrapper
function.
Here is how it works. For each function (e.g., exp) in libm.a,
there is a core (IEEE) function which has three underscores
prepended to it (i.e., ___exp). The wrapper function will twist
the result of the core function to comply to the standard
specified by the value of a global variable __lib_version:
if __lib_version = _IEEE_, return the result of the core function;
if __lib_version = _XOPEN_, return XOPEN result;
if __lib_version = _SVID3_, return SVID result;
if __lib_version = _POSIX_, return POSIX/ANSI result.
(See lib_version.h for definition of the macro.)
--------------------------------
4. HOW TO CREATE FDLIBM's libm.a
--------------------------------
There are two types of libm.a. One is IEEE only, and the other is
multi-standard compliant (supports IEEE,XOPEN,POSIX/ANSI,SVID).
To create the IEEE only libm.a, use
make "CFLAGS = -D_IEEE_LIBM"
This will create an IEEE libm.a, which is smaller in size, and
somewhat faster.
To create a multi-standard compliant libm, use
make "CFLAGS = -D_IEEE_SOURCE" --- multi-standard fdlibm: default
to IEEE
make "CFLAGS = -D_XOPEN_SOURCE" --- multi-standard fdlibm: default
to X/OPEN
make "CFLAGS = -D_POSIX_SOURCE" --- multi-standard fdlibm: default
to POSIX/ANSI
make "CFLAGS = -D_SVID3_SOURCE" --- multi-standard fdlibm: default
to SVID
There is a global variable (__lib_version) to control which
standard to follow. See lib_version.h.
Example:
Here is how one makes a SVID compliant libm.
First make the library by
make "CFLAGS = -D_SVID3_SOURCE".
The libm.a of FDLIBM will be multi-standard compliant and the
variable __lib_version is initialized to the value _SVID_ (see
lib_version.h).
=========== begin of sample.c ==============
main()
{
double y0();
printf("y0(1e300) = %1.20e\n",y0(1e300));
exit(0);
}
=========== end of sample.c ==============
Compile and run the sample.c to get the (default) SVID required result:
% cc sample.c libm.a
% a.out
y0: TLOSS error
y0(1e300) = 0.00000000000000000000e+00
Given a multi-standard libm.a, one can change the default standard as
follow:
xopen_version.c
===============================
#include "lib_version.h" /* FDLIBM's include files */
int __lib_version = _XOPEN_;
===============================
% cc xopen_version.c sample.c libm.a
% a.out
y0(1e300) = 0.00000000000000000000e+00 /* XOPEN answer */
ieee_version.c
===============================
#include "lib_version.h" /* FDLIBM's include files */
int __lib_version = _IEEE_;
===============================
% cc ieee_version.c sample.c libm.a
% a.out
y0(1e300) = -1.36813604503424810557e-151 /* IEEE answer */
Warning: changing the value of __lib_version affects only the
multi-standard compliant libm.a. It will have no effect on an
IEEE libm.a (created by setting the compiler flags -D_IEEE_LIBM).
------------------------------------------
5. IMPROVEMENTS IN BETA RELEASE (VS ALPHA)
------------------------------------------
+ fixed bugs in pow, tanh
+ fixed lgamma's monotonicity problem,
+ better algorithm for erf, erfc
+ better algorithm for acos and asin (error < 1ulp)
+ tuned up asinh, acosh, atanh, sinh, tanh, and jn
--------------
6. PROBLEMS ?
--------------
Please send comments and bug report to:
fdlibm-commentsasunpro.eng.sun.com
More information about the Numeric-interest
mailing list