Alpha release of Freely-Distributable libm available for testing

[David Hough]

K-C Ng has built "fdlibm" to the point where he would like to obtain outside
testing assistance.   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 it's had some testing.   
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 262K bytes, including an unconscionable
amount devoted to meeting the requirements of various contradictory 
standards test suites.

Here's K-C's Readme:

	***************************************
	* Announcing FDLIBM Version 1.0 alpha *
	***************************************
============================================================
			FDLIBM
============================================================
  (developed at SunPro, a Sun Microsystems, Inc. business.)
	      Version 1.0 alpha, 11/19/92


FDLIBM (Freely Distributable LIBM) is a C math library 
for machines that support IEEE floating-point arithmetic. 
In this release, only double precision is supported.

FDLIBM is intended to provide a reasonably portable, reference 
quality math library (libm.a). For an early release of FDLIBM, 
please send e-mail to 
		david.houghaEng.sun.com

--------------
1. ASSUMPTIONS
--------------
FDLIBM assumes:
 a.  IEEE style (if not precise compliance) arithmetic
 b.  32 bit integer arithmetic,
 c.  each floating-point number must be in IEEE 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 an -infinity is returned by 
default.

Similarly, exp(-huge) will be mapped to tiny*tiny to generate
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 underscore
prepend 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 fdlim: default
					 to IEEE
    make "CFLAGS = -D_XOPEN_SOURCE"  --- multi-standard fdlim: default
					 to X/OPEN
    make "CFLAGS = -D_POSIX_SOURCE"  --- multi-standard fdlim: default
					 to POSIX/ANSI
    make "CFLAGS = -D_SVID3_SOURCE"  --- multi-standard fdlim: 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. PROBLEMS ?
--------------
Please send comments and bug report to: 
		fdlibm-commentsasunpro.eng.sun.com