Alternative 80386 math lib v2.0 part06/06

Mon Dec 17 08:09:16 AEST 1990

Submitted-by: root at trantor
Archive-name: mathlib2.0/part06

---- Cut Here and feed the following to sh ----
#!/bin/sh
# this is mathlib.06 (part 6 of mathlib2.0)
# do not concatenate these parts, unpack them in order with /bin/sh
# file COPYING continued
#
if test ! -r _shar_seq_.tmp; then
	echo 'Please unpack part 1 first!'
	exit 1
fi
(read Scheck
 if test "$Scheck" != 6; then
	echo Please unpack part "$Scheck" next!
	exit 1
 else
	exit 0
 fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
	echo 'x - still skipping COPYING'
else
echo 'x - continuing file COPYING'
sed 's/^X//' << 'SHAR_EOF' >> 'COPYING' &&
X
X    This program is free software; you can redistribute it and/or modify
X    it under the terms of the GNU General Public License as published by
X    the Free Software Foundation; either version 1, or (at your option)
X    any later version.
X
X    This program is distributed in the hope that it will be useful,
X    but WITHOUT ANY WARRANTY; without even the implied warranty of
X    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
X    GNU General Public License for more details.
X
X    You should have received a copy of the GNU General Public License
X    along with this program; if not, write to the Free Software
X    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
X
Also add information on how to contact you by electronic and paper mail.
X
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
X
X    Gnomovision version 69, Copyright (C) 19xx name of author
X    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
X    This is free software, and you are welcome to redistribute it
X    under certain conditions; type `show c' for details.
X
The hypothetical commands `show w' and `show c' should show the
appropriate parts of the General Public License.  Of course, the
commands you use may be called something other than `show w' and `show
c'; they could even be mouse-clicks or menu items--whatever suits your
program.
X
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary.  Here a sample; alter the names:
X
X  Yoyodyne, Inc., hereby disclaims all copyright interest in the
X  program `Gnomovision' (a program to direct compilers to make passes
X  at assemblers) written by James Hacker.
X
X  <signature of Ty Coon>, 1 April 1989
X  Ty Coon, President of Vice
X
That's all there is to it!
SHAR_EOF
echo 'File COPYING is complete' &&
chmod 0644 COPYING ||
echo 'restore of COPYING failed'
Wc_c="`wc -c < 'COPYING'`"
test 12488 -eq "$Wc_c" ||
	echo 'COPYING: original size 12488, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= COPYRIGHT ==============
if test -f 'COPYRIGHT' -a X"$1" != X"-c"; then
	echo 'x - skipping COPYRIGHT (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting COPYRIGHT (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'COPYRIGHT' &&
All the assembler source files and the C source files ieee_retro.c, nextafter.c
and nextafterf.c are copyright 1990 to G. Geers.
X
The following C source files were written by W.J. Cody of Argonne National
Laboratory: erf.c, erff.c, gamma.c, gammaf.c, j0.c, j0f.c, j1.c, j1f.c,
lgamma.c and lgammaf.c. The float versions (names with f appended) were 
derived from the originals by G. Geers. These files are used with permission.
X
The original BASIC version of paranoia is copyright Professor W. M. Kahan.
SHAR_EOF
chmod 0644 COPYRIGHT ||
echo 'restore of COPYRIGHT failed'
Wc_c="`wc -c < 'COPYRIGHT'`"
test 504 -eq "$Wc_c" ||
	echo 'COPYRIGHT: original size 504, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= Makefile ==============
if test -f 'Makefile' -a X"$1" != X"-c"; then
	echo 'x - skipping Makefile (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting Makefile (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'Makefile' &&
# This file is part of the 80386 alternative math library and is covered by
# the GNU General Public Licence
X
GCC=gcc
CC=cc
CFLAGS=-O -fstrength-reduce -DIEEE
LIBDIR=/lib
INCDIR=/usr/include
# if you want to use setinternal in paranoia; define PTEST
#PTEST=
PTEST=-DTEST
VERS=2.0
X
s.o:
X	$(GCC) -c $<
X
c.o:
X	$(GCC) $(CFLAGS) -c $<
X
CSRC=j0.c lgamma.c gamma.c j1.c erf.c nextafter.c
X
CFSRC=j0f.c lgammaf.c gammaf.c j1f.c erff.c nextafterf.c
X
SSRC=acos.s copysign.s drem.s fabs.s hypot.s logb.s scalb.s tan.s \
X    asin.s ceil.s exp.s expm1.s finite.s log.s log1p.s \
X    pow.s sin.s atan.s cos.s exp10.s floor.s log10.s rint.s sqrt.s \
X    exp2.s log2.s sinh.s cosh.s tanh.s \
X    asinh.s acosh.s atanh.s atan2.s fmod.s ieee_ext.s \
X    infinity.s sqrtp.s ieee_values.s
X
SFSRC=acosf.s copysignf.s fabsf.s log10f.s sinf.s acoshf.s cosf.s \
X	finitef.s log1pf.s sinhf.s asinf.s coshf.s floorf.s log2f.s \
X	sqrtf.s asinhf.s dremf.s fmodf.s logbf.s sqrtpf.s atan2f.s \
X	exp10f.s hypotf.s logf.s tanf.s atanf.s exp2f.s ieee_extf.s \
X	powf.s tanhf.s atanhf.s expf.s ieee_valuesf.s rintf.s ceilf.s \
X	expm1f.s infinityf.s scalbf.s
X
COMMONSRC=ieee_retro.c _getsw.s setcont.s setinternal.s
X
ADD=CHANGELOG COPYING COPYRIGHT Makefile PROBLEMS README TODO d2dcomb.summ
X
ALLSRC=$(CSRC) $(CFSRC) $(SSRC) $(SFSRC) $(COMMONSRC) $(TESTSRC) $(ADD)
X
LIBOBJ=acos.o atanh.o erf.o hypot.o log10.o sin.o \
X       acosh.o ceil.o exp.o j0.o logb.o sinh.o \
X       asin.o copysign.o exp10.o expm1.o j1.o sqrt.o \
X       asinh.o cos.o fabs.o pow.o tan.o \
X       atan.o cosh.o finite.o lgamma.o rint.o tanh.o \
X       atan2.o drem.o floor.o log.o scalb.o log1p.o gamma.o \
X       exp2.o log2.o fmod.o ieee_ext.o \
X       infinity.o sqrtp.o ieee_values.o nextafter.o 
X
LIBFOBJ=acosf.o atanhf.o erff.o hypotf.o log10f.o sinf.o \
X       acoshf.o ceilf.o expf.o j0f.o logbf.o sinhf.o \
X       asinf.o copysignf.o exp10f.o expm1f.o j1f.o sqrtf.o \
X       asinhf.o cosf.o fabsf.o powf.o tanf.o \
X       atanf.o coshf.o finitef.o lgammaf.o rintf.o tanhf.o \
X       atan2f.o dremf.o floorf.o logf.o scalbf.o log1pf.o gammaf.o \
X       exp2f.o log2f.o fmodf.o ieee_extf.o \
X       infinityf.o sqrtpf.o ieee_valuesf.o nextafterf.o 
X
COMMONOBJ=ieee_retro.o _getsw.o setcont.o setinternal.o
X
TESTSRC=paranoia.c
X
all: libfpu libffpu libcfpu test ftest
X
libfpu: $(LIBOBJ) $(COMMONOBJ)
X	ar vr libfpu.a $(LIBOBJ) $(COMMONOBJ)
X
libffpu: $(LIBFOBJ) $(COMMONOBJ)
X	ar vr libffpu.a $(LIBFOBJ) $(COMMONOBJ)
X
libcfpu: $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ)
X	ar vr libcfpu.a $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ)
X
test: libfpu $(TESTOBJ)
X	$(GCC) -o paranoia $(TESTSRC) -DTEST libfpu.a
X
ftest: libffpu $(TESTOBJ)
X	$(GCC) -o fparanoia -DTEST -DSingle $(TESTSRC) libffpu.a
X
paranoia: test ftest
X
install: libfpu libffpu libcfpu
X	-cp libfpu.a $(LIBDIR)/libfpu.a
X	-chown bin $(LIBDIR)/libfpu.a
X	-chgrp bin $(LIBDIR)/libfpu.a
X	-chmod 444 $(LIBDIR)/libfpu.a
X	-cp libffpu.a $(LIBDIR)/libffpu.a
X	-chown bin $(LIBDIR)/libffpu.a
X	-chgrp bin $(LIBDIR)/libffpu.a
X	-chmod 444 $(LIBDIR)/libffpu.a
X	-cp libcfpu.a $(LIBDIR)/libcfpu.a
X	-chown bin $(LIBDIR)/libcfpu.a
X	-chgrp bin $(LIBDIR)/libcfpu.a
X	-chmod 444 $(LIBDIR)/libcfpu.a
X	-cp fpumath.h $(INCDIR)/fpumath.h
X
clean:
X	rm -f $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ) $(TESTOBJ) core a.out
X
clobber:
X	rm -f $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ) $(TESTOBJ) libfpu.a \
X	libffpu.a libcfpu.a paranoia fparanoia core a.out mathlib.0[1-6]
X
shar:
X	shar -c -L 55 -a -n mathlib$(VERS) -o mathlib $(ALLSRC) fpumath.h
X
# Dependencies generated using gcc -MM *.c
erf.o : erf.c 
gamma.o : gamma.c 
j0.o : j0.c
j1.o : j1.c
lgamma.o : lgamma.c
X	$(CC) -O -c lgamma.c
ieee_retro.o : ieee_retro.c
nextafter.o : nextafter.c fpumath.h
erff.o : erff.c fpumath.h
gammaf.o : gammaf.c  fpumath.h
j0f.o : j0f.c fpumath.h
j1f.o : j1f.c fpumath.h
nextafterf.o : nextafterf.c fpumath.h
lgammaf.o : lgammaf.c fpumath.h
paranoia.o : paranoia.c 
SHAR_EOF
chmod 0644 Makefile ||
echo 'restore of Makefile failed'
Wc_c="`wc -c < 'Makefile'`"
test 3879 -eq "$Wc_c" ||
	echo 'Makefile: original size 3879, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= PROBLEMS ==============
if test -f 'PROBLEMS' -a X"$1" != X"-c"; then
	echo 'x - skipping PROBLEMS (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting PROBLEMS (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'PROBLEMS' &&
exp() (and related exp functions) and pow() give rise to LOS exceptions. I
can't think of any other way of calculating these functions. The results I get
agree with the results from a SUN 4 and with Abramowitz and Stegun 
(to 15 places). If anyone knows of another algorithn please let me know.
X
There may still be problems with some functions for pathological arguments.
In some cases I should break the range of the argument up and use
approximations for large and/or small arguments. In fact, the ESIX supplied
exp() is just as bad (good?!?) as mine :-)
SHAR_EOF
chmod 0644 PROBLEMS ||
echo 'restore of PROBLEMS failed'
Wc_c="`wc -c < 'PROBLEMS'`"
test 557 -eq "$Wc_c" ||
	echo 'PROBLEMS: original size 557, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= README ==============
if test -f 'README' -a X"$1" != X"-c"; then
	echo 'x - skipping README (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README' &&
VERSION: This is version 2.0
--------
This is version 2.0 of the alternative 386 math library. It supplants
all previous versions.
X
Version Notes
-------------
The main change between 1.1 and 2.0 is that all functions are now available
in float versions. To use them you need gcc/gas in order to create the
library (libffpu.a) and an ANSI C compiler to compile your source code because
K&R-1 compilers pass all floating point argumnets to functions as doubles and
this is not the desired behaviour.
X
The float function names are formed from the double versions by appending
an 'f', for example sin -> sinf. These functions are fully prototyped in
fpumath.h.
X
Note that jn and yn are not available at this time.
X
Three libraries are now built:
libfpu.a - double version
libffpu.a - float version
libcfpu.a - combined version 
X
The speed difference between floats and doubles is dependent on the
instruction mix - from memory operands or on stack operands. Floats are
faster to load from memory but once on the stack, all operations occur on
the 80 bit temporary real representation.
X
To use the float version effectively you really need an ANSI stdio
implementation.
X
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 typically from 2 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.
X
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.
X
I have coded the additional IEEE 754 required functions to provide a
conforming double precision implementation. 
X
You require gcc/gas in order to compile the assembler source code.
X
People who are using SUN 386i's (Roadrunner) may also find this useful.
X
I have found that writing good low level numerical code is far harder
than I initially thought it would be. But I've learnt a lot, and I still 
enjoy doing it.
X
File List
---------
CHANGELOG       atanhf.s        expm1f.s        j0.c            rint.s
COPYING         ceil.s          fabs.s          j0f.c           rintf.s
COPYRIGHT       ceilf.s         fabsf.s         j1.c            scalb.s
Makefile        copysign.s      finite.s        j1f.c           scalbf.s
PROBLEMS        copysignf.s     finitef.s       lgamma.c        setcont.s
README          cos.s           floor.s         lgammaf.c       setinternal.s
TODO            cosf.s          floorf.s        log.s           sin.s
_getsw.s        cosh.s          fmod.s          log10.s         sinf.s
acos.s          coshf.s         fmodf.s         log10f.s        sinh.s
acosf.s         d2dcomb.summ    fpumath.h       log1p.s         sinhf.s
acosh.s         drem.s          gamma.c         log1pf.s        sqrt.s
acoshf.s        dremf.s         gammaf.c        log2.s          sqrtf.s
asin.s          erf.c           hypot.s         log2f.s         sqrtp.s
asinf.s         erff.c          hypotf.s        logb.s          sqrtpf.s
asinh.s         exp.s           ieee_ext.s      logbf.s         tan.s
asinhf.s        exp10.s         ieee_extf.s     logf.s          tanf.s
atan.s          exp10f.s        ieee_retro.c    nextafter.c     tanh.s
atan2.s         exp2.s          ieee_values.s   nextafterf.c    tanhf.s
atan2f.s        exp2f.s         ieee_valuesf.s  paranoia.c
atanf.s         expf.s          infinity.s      pow.s
atanh.s         expm1.s         infinityf.s     powf.s
X
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.
X
The file d2dcomb.summ contains a summary of some test results. See that
file for details.
X
The IEEE specified functions with which you may not be familiar are:
X
double copysign(x,y)
double x,y;
copysign returns x with the sign of y. IEEE denormal will be set if x is a
denormal.
X
double drem(x)
double x;
drem returns the IEEE remainder of x/y - it may be slow.
X
int finite(x)
double x;
finite returns true if -inf < x < inf; false otherwise. Does not raise any
floating point exceptions.
X
double logb(x)
double x;
logb returns the unbiased exponent of its argument.
X
double rint(x)
double x;
rint returns its argument rounded in the prevailing rounding mode.
X
double scalb(x, n)
double x;
int n;
scalb returns x*2^n.
X
Most of the above are just hooks directly into functions provided as single
instructions in the 80387.
X
double infinity()
infinity returns +inf. No exceptions are raised.
X
int isnan(x)
double x;
isnan returns 1 if x is a nan; 0 otherwise. Ieee exceptions are not
affected.
X
int isnormal(x)
double x;
isnormal returns 1 if x is a normalized number; 0 otherwise. Ieee exceptions 
are not affected.
X
int issubnormal(x)
double x;
issubnormal returns 1 if x is not a normalized number; 0 otherwise. Ieee 
exceptions are not affected.
X
int iszero(x)
double x;
iszero returns 1 if x is +/- 0.0; 0 otherwise. Ieee exceptions are not affected.
X
int isinf(x)
double x;
isinf returns 1 if x is +/- inf; 0 otherwise. Ieee exceptions are not affected.
X
int signbit(x)
double x;
signbit return the sign of x. 1 if negative and 0 if positive. Ieee
exceptions are not affected.
X
The is... functions and signbit are in the file ieee_ext.s
X
void ieee_retrospective(f)
FILE *f;
ieee_retrospective prints a list of IEEE exceptions that are currently
active on the 80387 in the file pointed to by f.
X
double
max_normal()
max_normal returns the maximum +ve normalized number. No exceptions are
raised.
X
double
min_normal()
mmin_normal returns the minimum +ve normalized number. No exceptions are
raised.
X
double
max_subnormal()
max_subnormal returns the largest +ve denormalized number. This raises the
denormal exception because of the way floating point numbers are returned.
X
double
min_subnormal()
min_subnormal returns the minimum +ve denormalized number. This raises the
denormal exception because of the way floating point numbers are returned.
X
double
quiet_nan(n)
long n;
quiet_nan returns a quiet nan. The argument is ignored. No exceptions are
raised.
X
double
signaling_nan(n)
long n;
signaling_nan returns a signaling nan. The argument is ignored. This raises
the invalid operation exception because of the way floating point numbers
are returned.
X
double
nextafter(x,y)
double x,y;
nextafter returns the next representable floating point number from x in
the direction of y. Exceptions may be raised depending on the arguments.
X
double erfcx(x)
double x;
erfcx returns x^2 * erfc(x).
X
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:
X
int setinternal()
X
and is declared for your convenience in fpumath.h. The return value is the
current control word. 
X
You can set the control word using setcont(mode). Again, this is defined in
fpumath.h:
X
int setcont(mode)
int mode;
X
This is really only provided to reset the original mode obtained using 
setinternal(). The previous mode is returned.
X
Note that more general forms of these functions are provided in the standard 
library using the fpsetmask and fpgetmask routines.
X
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.
X
double sqrtp(x)
double x;
Returns the sqrt of x in 64 bit (double) mode irrespective of the current
precision.
X
Acknowledgements
----------------
I'd like to thank the developers of the Berkeley 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. And also Dan Lau
of Intel. Also thanks go to W.J. Cody of Argonne National Laboratory
(USA).
X
And to everybody who has sent me feedback.
X
X
I have decided to cover the parts I wrote by the GNU General Public Licence 
with my modification (see below), a copy of which is included with this 
distribution.
X
In a nutshell, the Berkeley code is covered by their licence. My code is covered
by GNU's with my modification (see below). This infringes on nobodies rights.
X
Amendment to the GNU General Public Licence (*ONLY* applicable to my code)
-------------------------------------------
I may choose to cover this code by the Free Software Foundation's General
Public Library Licence when it becomes available. At present I make the 
following amendment to the licencing agreement:
X
*******************************************************************************
You may use this library in products which are distributed in binary format
only as long as you provide source code for the library with the distribution
or will supply the source code for the library for a period of 3 years from the
date of providing the binary files.
********************************************************************************
X
This in no way changes the GNU licence as applicable to other programs.
X
Also I have not modified the GPL it is ditributed in its original form as
it must be.
X
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).
X
X
X
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!), gcc 1.37.1 and gas 1.36 with COFF/stab patches.
X
I have also modfied the f2c libF77 to use setcont() and ieee_retrospective()
which makes the whole f2c environment on a 386 look like SUN FORTRAN. If
anyone is interested in these (trivial) changes drop me a line.
X
Please mail comments and bug reports to glenn at qed.physics.su.oz.au.
X
X			Share and Enjoy,
X				Glenn
X
You can also ftp this stuff from suphys.physics.su.oz.au.
X
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'
X
Ph: +61 2 692-3241
SHAR_EOF
chmod 0644 README ||
echo 'restore of README failed'
Wc_c="`wc -c < 'README'`"
test 11151 -eq "$Wc_c" ||
	echo 'README: original size 11151, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= TODO ==============
if test -f 'TODO' -a X"$1" != X"-c"; then
	echo 'x - skipping TODO (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting TODO (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'TODO' &&
0. Keep modifying the sources to improve performance.
1. Modify the assembler source so that it works with the standard 'as'.
2. Add some more esoteric functions (suggestions?).
3. Provide jn() and yn().
SHAR_EOF
chmod 0644 TODO ||
echo 'restore of TODO failed'
Wc_c="`wc -c < 'TODO'`"
test 204 -eq "$Wc_c" ||
	echo 'TODO: original size 204, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= d2dcomb.summ ==============
if test -f 'd2dcomb.summ' -a X"$1" != X"-c"; then
	echo 'x - skipping d2dcomb.summ (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting d2dcomb.summ (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'd2dcomb.summ' &&
Comments
--------
These tests were done using routines distributed with the Portable Math 
Library by Fred Fish. I have a sneaking feeling that his exp, atanh and 
atan return bad results (double precision FORTRAN library under TOPS-20) 
since the '387 and SUN 4 results agree. I have partially verified this 
using the extended tables in Abramowitz & Stegun.
X
Results from libfpu.a
dabs:	maximum relative error 0.000000e+00
acos:	maximum relative error 1.395807e-16
acosh:	maximum relative error 1.686042e-16
asin:	maximum relative error 1.378420e-16
asinh:	maximum relative error 1.181272e-13
atan:	maximum relative error 6.406394e-12
atan2:	maximum relative error 2.261728e-15
atanh:	maximum relative error 3.731406e-13
cos:	maximum relative error 0.000000e+00
cosh:	maximum relative error 3.828768e-15
exp:	maximum relative error 4.165239e-13
log:	maximum relative error 0.000000e+00
log10:	maximum relative error 0.000000e+00
sin:	maximum relative error 1.855096e-16
sinh:	maximum relative error 3.828768e-15
sqrt:	maximum relative error 1.252752e-16
tan:	maximum relative error 1.110223e-16
tanh:	maximum relative error 1.201235e-16
X
Results from libfpu.a with extended precision
dabs:	maximum relative error 0.000000e+00
acos:	maximum relative error 1.892807e-16
acosh:	maximum relative error 0.000000e+00
asin:	maximum relative error 1.431811e-16
asinh:	maximum relative error 7.000127e-16
atan:	maximum relative error 6.406394e-12
atan2:	maximum relative error 2.261728e-15
atanh:	maximum relative error 3.731406e-13
cos:	maximum relative error 0.000000e+00
cosh:	maximum relative error 0.000000e+00
exp:	maximum relative error 4.176205e-13
log:	maximum relative error 0.000000e+00
log10:	maximum relative error 0.000000e+00
sin:	maximum relative error 1.855096e-16
sinh:	maximum relative error 0.000000e+00
sqrt:	maximum relative error 1.252752e-16
tan:	maximum relative error 1.110223e-16
tanh:	maximum relative error 1.226565e-16
X
Results form libm.a (Esix standard) - deleted functions do not exist
dabs:	maximum relative error 0.000000e+00
acos:	maximum relative error 1.892807e-16
asin:	maximum relative error 1.431811e-16
atan:	maximum relative error 6.406394e-12
atan2:	maximum relative error 2.261728e-15
cos:	maximum relative error 0.000000e+00
cosh:	maximum relative error 3.828768e-15
exp:	maximum relative error 4.165239e-13
log:	maximum relative error 0.000000e+00
log10:	maximum relative error 0.000000e+00
sin:	maximum relative error 1.855096e-16
sinh:	maximum relative error 3.828768e-15
sqrt:	maximum relative error 1.252752e-16
tan:	maximum relative error 1.110223e-16
tanh:	maximum relative error 1.201235e-16
X
Results form libm.a (SUN 4)
dabs:	maximum relative error 0.000000e+00
acos:	maximum relative error 0.000000e+00
acosh:	maximum relative error 1.686042e-16
asin:	maximum relative error 1.378420e-16
asinh:	maximum relative error 7.000127e-16
atan:	maximum relative error 6.406394e-12
atan2:	maximum relative error 2.261728e-15
atanh:	maximum relative error 3.731406e-13
cos:	maximum relative error 1.156277e-16
cosh:	maximum relative error 1.651640e-16
exp:	maximum relative error 4.176205e-13
log:	maximum relative error 0.000000e+00
log10:	maximum relative error 0.000000e+00
sin:	maximum relative error 1.855096e-16
sinh:	maximum relative error 0.000000e+00
sqrt:	maximum relative error 1.252753e-16
tan:	maximum relative error 1.425732e-16
tanh:	maximum relative error 0.000000e+00
SHAR_EOF
chmod 0644 d2dcomb.summ ||
echo 'restore of d2dcomb.summ failed'
Wc_c="`wc -c < 'd2dcomb.summ'`"
test 3424 -eq "$Wc_c" ||
	echo 'd2dcomb.summ: original size 3424, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= fpumath.h ==============
if test -f 'fpumath.h' -a X"$1" != X"-c"; then
	echo 'x - skipping fpumath.h (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting fpumath.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fpumath.h' &&
/*
** This file is covered by the GNU General Public Licence
** 
** This file should be installed somewhere in your standard include
** search path and given the name fpumath.h.
*/
X
/*
** This is a good way of flagging whether 
** you are using the alternate stuff
*/
#ifndef	__FPUMATH_H
#define	__FPUMATH_H
X
#define __FPU__
X
X
#if defined(__GNUC__) || defined(__STDC__)
X
extern double j0(double), j1(double), jn(int,double), y0(double), y1(double);
extern double yn(int,double);
extern double erf(double), erfc(double), erfcx(double);
extern double exp(double), log(double), log10(double), pow(double,double);
extern double sqrt(double);
extern double expm1(double), log1p(double), exp10(double);
extern double exp2(double), log2(double);
extern double floor(double), ceil(double), fabs(double);
extern double lgamma(double), gamma(double);
extern double sinh(double), cosh(double), tanh(double);
extern double asinh(double), acosh(double), atanh(double);
extern double sin(double), cos(double), tan(double), asin(double);
extern double acos(double), atan(double), atan2(double,double);
extern double hypot(double,double);
extern double sqrtp(double);
X
#else
X
extern double j0(), j1(), jn(), y0(), y1(), yn();
extern double erf(), erfc(), erfcx();
extern double exp(), log(), log10(), pow(), sqrt();
extern double expm1(), log1p(), exp10();
extern double exp2(), log2();
extern double floor(), ceil(), fabs();
extern double lgamma(), gamma();
extern double sinh(), cosh(), tanh();
extern double asinh(), acosh(), atanh();
extern double sin(), cos(), tan(), asin();
extern double acos(), atan(), atan2(), hypot();
extern double sqrtp();
X
#endif
X
/*
** IEEE functions
*/
X
#if defined(__GNUC__) || defined(__STDC__)
X
extern double rint(double), drem(double,double), scalb(double,int);
extern double logb(double), copysign(double,double);
extern double infinity(void), nextafter(double,double);
extern int isnan(double), iszero(double), isinf(double), isnormal(double);
extern int issubnormal(double);
extern int signbit(double), finite(double);
extern double max_normal(void), min_normal(void), min_subnormal(void);
extern double max_subnormal(void);
extern double quiet_nan(long int),signaling_nan(long int);
X
#else
X
extern double rint(), drem(), scalb(), logb(), copysign();
extern double infinity(), nextafter();
extern int isnan(), iszero(), isinf(), isnormal(), issubnormal();
extern int signbit(), finite();
extern double max_normal(), min_normal(), min_subnormal(), max_subnormal();
extern double quiet_nan(),signaling_nan();
X
#endif
X
/*
** Float versions
*/
#if defined(__GNUC__) || defined(__STDC__)
X
extern float j0f(float), j1f(float), jnf(int,float), y0f(float);
extern float y1f(float), ynf(int,float);
extern float erff(float), erfcf(float), erfcxf(float), sqrtf(float);
extern float expf(float), logf(float), log10f(float), powf(float, float);
extern float expm1f(float), log1pf(float), exp10f(float);
extern float exp2f(float), log2f(float);
extern float floorf(float), ceilf(float), fabsf(float);
extern float lgammaf(float), gammaf(float);
extern float sinhf(float), coshf(float), tanhf(float);
extern float asinhf(float), acoshf(float), atanhf(float);
extern float sinf(float), cosf(float), tanf(float), asinf(float);
extern float acosf(float), atanf(float), atan2f(float, float);
extern float hypotf(float, float);
extern float sqrtpf(float);
X
extern float rintf(float), dremf(float), scalbf(float, int);
extern float logbf(float), copysignf(float, float);
extern float infinityf(void), nextafterf(float, float);
extern int isnanf(float), iszerof(float), isinff(float);
extern int isnormalf(float), issubnormalf(float);
extern int signbitf(float), finitef(float);
extern float max_normalf(void), min_normalf(void);
extern float min_subnormalf(void), max_subnormalf(void);
extern float quiet_nanf(long int), signaling_nanf(long int);
X
#endif
X
/*
** Extensions
*/
X
#if defined(__GNUC__) || defined(__STDC__)
X
extern int setinternal(void), setcont(int);
X
#else
X
extern int setinternal(), setcont();
X
#endif
X
/*
** The argument of ieee_retrospective is a FILE * 
** but we don't want to include stdio.h.
*/
X
extern void ieee_retrospective();
X
/*
** Useful defines for setcont
*/
X
#define MASK_ALL0 0x107f	/* single precision */
#define MASK_ALL 0x127f		/* double precision */
#define MASK_ALL1 0x137f 	/* extra precision */
X
#define __infinity infinity
X
#define M_E	2.7182818284590452354
#define M_LOG2E	1.4426950408889634074
#define M_LOG10E	0.43429448190325182765
#define M_LN2	0.69314718055994530942
#define M_LN10	2.30258509299404568402
#define M_PI	3.14159265358979323846
#define M_PI_2	1.57079632679489661923
#define M_PI_4	0.78539816339744830962
#define M_1_PI	0.31830988618379067154
#define M_2_PI	0.63661977236758134308
#define M_2_SQRTPI	1.12837916709551257390
#define M_SQRT2	1.41421356237309504880
#define M_SQRT1_2	0.70710678118654752440
X
#ifndef IEEE
#define MAXFLOAT	((float)3.40282346638528860e+38)
#define	MINFLOAT	((float)1.40129846432481707e-45)
#else
#define	MAXFLOAT	(max_normalf())
#define	MINFLOAT	(min_subnormalf())
#endif
X
#ifndef IEEE
#define	MAXDOUBLE	1.79769313486231570e+308  /* wrong in math.h */
#define	MINDOUBLE	4.94065645841246544e-324
#else
#define	MAXDOUBLE	(max_normal())
#define	MINDOUBLE	(min_subnormal())
#endif
X
#define ULP		1.1102230246251565e-16
#define ULPF		5.9604644775390625e-08
X
#define	HUGE_VAL	MAXFLOAT
X
#ifdef IEEE
#define HUGE	(__infinity())
#else
#define HUGE 	MAXDOUBLE
#endif
X
#endif
SHAR_EOF
chmod 0644 fpumath.h ||
echo 'restore of fpumath.h failed'
Wc_c="`wc -c < 'fpumath.h'`"
test 5452 -eq "$Wc_c" ||
	echo 'fpumath.h: original size 5452, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
rm -f _shar_seq_.tmp
echo You have unpacked the last part
exit 0
--
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'