Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!swrinde!zaphod.mps.ohio-state.edu!caen!ox.com!emv From: glenn@qed.physics.su.OZ.AU (Glenn Geers) Newsgroups: comp.archives Subject: [sysv386...] 386/387 Math library available for ftp Keywords: Math library Message-ID: <1990Dec6.053508.3576@ox.com> Date: 6 Dec 90 05:35:08 GMT References: <1990Dec3.104538.16719@metro.ucc.su.OZ.AU> Sender: emv@ox.com (Edward Vielmetti) Reply-To: glenn@qed.physics.su.OZ.AU (Glenn Geers) Followup-To: comp.unix.sysv386,comp.unix.xenix.sco,comp.os.minix Organization: School of Physics, Uni of Sydney, Australia. Lines: 243 Approved: emv@ox.com (Edward Vielmetti) X-Original-Newsgroups: comp.unix.sysv386,comp.unix.xenix.sco,comp.os.minix Archive-name: math/ieee-float/387-library/1990-12-03 Archive-directory: suphys.physics.su.oz.au:/Esix/src/math/ [129.78.129.1] Original-posting-by: glenn@qed.physics.su.OZ.AU (Glenn Geers) Original-subject: 386/387 Math library available for ftp Reposted-by: emv@ox.com (Edward Vielmetti) Hi, the final version of my maths library is available for anonymous ftp from suphys.physics.su.oz.au (129.78.129.1) in the directory Esix/src/math. I have covered it by the GNU General Public Licence with one small modification.This mod is noted in the README and in each of my source files. I'd like to thank Dan Lau of Intel for some useful feedback. Hopefully this code is now bug free :-) The README is included with this posting. I will post this to alt.sources next week (I'm attending a conference this week)and to comp.sources.misc when that group comes to life again. Share and 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 you're 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 exp2.s j0.c rint.s COPYING atan2.s expm1.s j1.c scalb.s COPYRIGHT atanh.s fabs.s jn.c setcont.s Makefile ceil.s finite.s lgamma.c setinternal.s PROBLEMS copysign.s floor.s log.s sin.s README cos.s fmod.s log10.s sinh.s TODO cosh.s fpumath.h log1p.s sqrt.s _getsw.s d2dcomb.summ gamma.c log2.s sqrtp.s acos.s drem.s hypot.s logb.s tan.s acosh.s erf.c ieee_ext.s mathimpl.h tanh.s asin.s exp.s ieee_retro.c paranoia.c asinh.s exp10.s infinity.s pow.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 file d2dcomb.summ contains a summary of some test results. See that file for details. 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. IEEE denormal will be set if x is a denormal. 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. Does not raise any floating point exceptions. 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 the above are just hooks directly into functions provided as single instructions in the 80387. double infinity() infinity returns +inf. No exceptions are raised. int isnan(x) double x; isnan returns 1 if x is a nan; 0 otherwise. Ieee exceptions are not affected. int isnormal(x) double x; isnormal returns 1 if x is a normalized number; 0 otherwise. Ieee exceptions are not affected. int issubnormal(x) double x; issubnormal returns 1 if x is not a normalized number; 0 otherwise. Ieee exceptions are not affected. int iszero(x) double x; iszero returns 1 if x is +/- 0.0; 0 otherwise. Ieee exceptions are not affected. int isinf(x) double x; isinf returns 1 if x is +/- inf; 0 otherwise. Ieee exceptions are not affected. int signbit(x) double x; signbit return the sign of x. 1 if negative and 0 if positive. Ieee exceptions are not affected. The is... functions and signbit are in the file ieee_ext.s 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. 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. double sqrtp(x) double x; Returns the sqrt of x in 64 bit (double) mode irrespective of the current precision. 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. And also Dan Lau of Intel. ******************************************************************************** ``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 with my modification (see below), 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 with my modification (see below). This infringes on nobodies rights. 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: ******************************************************************************* 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. ******************************************************************************** This in no way changes the GNU licence as applicable to other programs. 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!), gcc 1.37.1 and gas 1.36 with COFF/stab patches. 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. Please mail comments and bug reports to glenn@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' Brought to you by Super Global Mega Corp .com