386/387 Maths Library

[Glenn Geers]

Hi,
	well the maths library I mentioned about a week ago is ready. Before
I post it (which group? - here, alt.sources, comp.sources.misc, comp.sources.
unix?) I'd like to get some feedback about posting libraries under the GNU
General Public Licence. I like GNU and support there ideals but I can also
appreciate software company X not wishing to provide source for a product
that has cost it Y dollars to produce. By putting my code under GNU I make
everyone who incorporates it in anything be obligated to provide source for
the product in which it is included. The alternative is to modify the GPL
in such a way that you can give away source code if you want but you are not
obligated to.
	Please *NO* flames. My aim is to deduce what licencing arrangements are best
for people in netland.

Anyway, I've included the README with this posting. It's changed quite a bit
since beta1.
					Enjoy,
						Glenn


Introduction
------------
The files in this directory consist of assembler and C source for an 
alternative maths library for Unices (including Xenix) running on an
80386/80387 combination and using gcc/gas as the compiler system. These
routines are from 5 to 10 times faster than those in the supplied maths
library; assuming that your library, like mine (ESIX rev. D and Xenix
2.3.2), does not use the '387 inbuilts to perform transcendental
calculations.

For those of you without a '387 you must use the full emulator and
consequently may not see any speed up. I haven't tried. Under Xenix you
must have a '387---some of the '387 instructions are not emulated. If you
have a '287 your in the same boat; some of the assembler routines won't
work.

I have coded the additional IEEE 754 required functions to provide a
conforming double precision implementation. 

You require gcc/gas in order to compile the assembler source code.

People who are using SUN 386i's (Roadrunner) may also find this useful.

File List
---------
CHANGELOG      atan.s         exp10.s        j1.c           pow.s
COPYING        atan2.s        exp2.s         jn.c           rint.s
COPYRIGHT      atanh.s        expm1.s        lgamma.c       scalb.s
Makefile       ceil.s         fabs.s         log.s          setcont.s
README         copysign.s     finite.s       log10.s        setinternal.s
TODO           cos.s          floor.s        log1p.s        sin.s
acos.s         cosh.s         fpumath.h      log2.s         sinh.s
acosh.s        drem.s         gamma.c        logb.s         sqrt.s
asin.s         erf.c          hypot.s        mathimpl.h     tan.s
asinh.s        exp.s          j0.c           paranoia.c     tanh.s

Comments
--------
The additional program paranoia.c attempts to tell how well your floating
point conforms to the IEEE standard. You may like to compare the output
when linked with your existing library and with this one. Paranoia cannot be
compiled using standard 'cc' on Esix 3.2 revision D.

The IEEE specified functions with which you may not be familiar are:

double copysign(x,y)
double x,y;
copysign returns x with the sign of y.

double drem(x)
double x;
drem returns the IEEE remainder of x/y - it may be slow.

int finite(x)
double x;
finite returns true if -inf < x < inf; false otherwise.

double logb(x)
double x;
logb returns the unbiased exponent of its argument.

double rint(x)
double x;
rint returns its argument rounded in the prevailing rounding mode.

double scalb(x, n)
double x;
int n;
scalb returns x*2^n.

Most of these are just hooks directly into functions provided as single
instructions in the 80387.

Additional Functions
--------------------
It is suggested in the book `Programming the 80386' by Crawford and Gelsinger
that all floating point exceptions except `invalid operation' be masked. 
That is the 80387's inbuilt exception handler should be used. If you want this 
behaviour call setinternal() at the start of your code. This function is defined
by:

int setinternal()

and is declared for your convenience in fpumath.h. The return value is the
current control word. 

You can set the control word using setcont(mode). Again, this is defined in
fpumath.h:

int setcont(mode)
int mode;

This is really only provided to reset the original mode obtained using 
setinternal(). The previous mode is returned.

Note that more general forms of these functions are provided in the standard 
library using the fpsetmask and fpgetmask routines.

If you use setinternal(), arithmetic operations like 1/0 will return
infinity - inf will be printed if you are printing the output. Note that 0/0 
will still produce a floating point exception; the value of this operation 
is undefined.

Acknowledgements
----------------
I'd like to thank the developers of the Berkley Software Distribution who
believe in software freedom and made my job easier by providing C source
for all the functions. I also thank the author of paranoia.c
********************************************************************************
``This product includes software developed by the University of California,
 Berkeley and its contributors''
********************************************************************************

I have decided to cover the parts I wrote by the GNU General Public Licence, a 
copy of which is included with this distribution. Although the BSD and GNU 
licences are different I don't really want people to split this distribution 
up. It's complete as it is.

In a nutshell, the Berkley code is covered by their licence. My code is covered
by GNU's. This infringes on nobodies rights.

Installation
------------
Edit the Makefile and change LIBDIR (default /lib) and INCDIR 
(default /usr/include) according to taste.
Type `make' to produce the library (libfpu.a) and paranoia.
Type `make install' to install libfpu.a in $(LIBDIR) and fpumath.h 
in $(INCDIR).



I'm interested in obtaining feedback on the performance of this library and
of being informed of any bugs. I will continue to support this as long as I
have a 32 bit Intel CPU running some form of UNIX and GNU C. My own
development system was ESIX System V.3.2 revision D running on a 20 MHz
386 (and 387 of course!).

Please mail comments and bug reports to glenn at qed.physics.su.oz.au.

			Share and Enjoy,
				Glenn

You can also ftp this stuff from suphys.physics.su.oz.au.

Glenn Geers                       | "So when it's over, we're back to people.
Department of Theoretical Physics |  Just to prove that human touch can have
The University of Sydney          |  no equal."
Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'

Ph: +61 2 692-3241



--
Glenn Geers                       | "So when it's over, we're back to people.
Department of Theoretical Physics |  Just to prove that human touch can have
The University of Sydney          |  no equal."
Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'