glenn@extro.ucc.su.OZ.AU (Glenn Geers) (12/08/90)
Here is my alternative math library for the 80386/80387. You need gcc/gas
in order to compile it but you can link with code compiled with cc.
I will post this to comp.sources.misc when that group starts running at speed
again.
Glenn
Submitted-by: glenn@trantor
Archive-name: libfpu/part01
#!/bin/sh
# This is libfpu, a shell archive (shar 3.47)
# made 12/07/1990 20:56 UTC by glenn@trantor
# Source directory /usr1/src/math/dist
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 3702 -rw-r--r-- CHANGELOG
# 12488 -rw-r--r-- COPYING
# 316 -rw-r--r-- COPYRIGHT
# 1984 -rw-r--r-- Makefile
# 295 -rw-r--r-- PROBLEMS
# 9579 -rw-r--r-- README
# 244 -rw-r--r-- TODO
# 320 -rw-r--r-- _getsw.s
# 544 -rw-r--r-- acos.s
# 398 -rw-r--r-- acosh.s
# 401 -rw-r--r-- asin.s
# 397 -rw-r--r-- asinh.s
# 331 -rw-r--r-- atan.s
# 931 -rw-r--r-- atan2.s
# 438 -rw-r--r-- atanh.s
# 682 -rw-r--r-- ceil.s
# 555 -rw-r--r-- copysign.s
# 320 -rw-r--r-- cos.s
# 486 -rw-r--r-- cosh.s
# 3424 -rw-r--r-- d2dcomb.summ
# 381 -rw-r--r-- drem.s
# 2681 -rw-r--r-- erf.c
# 400 -rw-r--r-- exp.s
# 404 -rw-r--r-- exp10.s
# 387 -rw-r--r-- exp2.s
# 418 -rw-r--r-- expm1.s
# 322 -rw-r--r-- fabs.s
# 447 -rw-r--r-- finite.s
# 628 -rw-r--r-- floor.s
# 380 -rw-r--r-- fmod.s
# 2314 -rw-r--r-- fpumath.h
# 604 -rw-r--r-- gamma.c
# 377 -rw-r--r-- hypot.s
# 1977 -rw-r--r-- ieee_ext.s
# 853 -rw-r--r-- ieee_retro.c
# 1295 -rw-r--r-- ieee_values.s
# 394 -rw-r--r-- infinity.s
# 5099 -rw-r--r-- j0.c
# 5169 -rw-r--r-- j1.c
# 2310 -rw-r--r-- jn.c
# 3382 -rw-r--r-- lgamma.c
# 329 -rw-r--r-- log.s
# 333 -rw-r--r-- log10.s
# 346 -rw-r--r-- log1p.s
# 329 -rw-r--r-- log2.s
# 339 -rw-r--r-- logb.s
# 2996 -rw-r--r-- mathimpl.h
# 1392 -rw-r--r-- nextafter.c
# 57545 -rw-r--r-- paranoia.c
# 2128 -rw-r--r-- pow.s
# 326 -rw-r--r-- rint.s
# 343 -rw-r--r-- scalb.s
# 377 -rw-r--r-- setcont.s
# 449 -rw-r--r-- setinternal.s
# 320 -rw-r--r-- sin.s
# 487 -rw-r--r-- sinh.s
# 359 -rw-r--r-- sqrt.s
# 605 -rw-r--r-- sqrtp.s
# 334 -rw-r--r-- tan.s
# 493 -rw-r--r-- tanh.s
#
# ============= CHANGELOG ==============
if test -f 'CHANGELOG' -a X"$1" != X"-c"; then
echo 'x - skipping CHANGELOG (File already exists)'
else
echo 'x - extracting CHANGELOG (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'CHANGELOG' &&
Thu Nov 15 07:05:39 EDT 1990
Shar'd beta version for release to testers.
X
Sat Nov 17 17:06:17 EDT 1990
Fixed floor and ceil so that they address their local variables at -ve
offsets relative to the frame pointer.
X
Fixed hypot.s so that unnecessary register loads are avoided. I wasn't aware
of the correct assembler mnemonic to multiply %st(0) by the contents of a memory
location.
X
Sat Nov 17 18:49:31 EDT 1990
Copysign rewritten in assembler.
X
Sun Nov 18 11:21:24 EDT 1990
Sinh and cosh now in assembler.
Tanh now in assembler.
Fixed nasty bug in acos.s.
exp__E no longer needed - source deleted.
X
Tue Nov 20 21:24:26 EDT 1990
Sinh and cosh now multiply by 0.5 instead of dividing by 2. It's faster.
Added setinternal for convenience in some cases. Using this makes the
output of paranoia closer to that obtained on a Sun 4.
Added setcont.
Modified the Makefile so that test and paranoia are synonymous.
Modified atan2.c to avoid a floating point register load.
Fixed bug in atan2 - was giving incorrect result when y/x = +/-inf.
X
Thu Nov 22 06:52:47 EDT 1990
Fixed local variable handling in assembler files. Stack needs to be
allocated for local variables.
All assembler routines are now .align 4.
X
Thu Nov 22 19:24:09 EDT 1990
Added MASK_ALL to fpumath.h
X
Fri Nov 23 21:11:33 EDT 1990
Inverse hyperbolics missing from fpumath.h - fixed.
Some definitions missing from atanh.c - didn't work on some systems.
Coded inverse hyperbolics in assembler.
Coded atan2 in assembler.
X
Sat Nov 24 17:12:34 EDT 1990
Modified hypot.s to avoid use of ffree.
atan2 was not working correctly - fixed
Ready for netwide release.
X
Mon Nov 26 19:25:56 EDT 1990
Put in explicit call to cc for compiling lgamma. This means that the
library can link to code not compiled with gcc (not quite true pow.s has to be
changed).
X
Wed Nov 28 22:57:43 EDT 1990
Coded the is??? functions required for IEEE conformance
Coded ieee_retrospective() which prints a summary of IEEE exceptions that are on
when the function is called.
X
Fri Nov 30 17:35:29 EDT 1990
Removed the last vestige of GNU specifics. You can now use the library with 'cc'
although you need gcc/gas to create it.
X
Sat Dec 1 15:48:19 EDT 1990
Found a bug in pow (pow(-x,x) shoudn't work for x non-integral) fixed.
X
Sat Dec 1 21:00:48 EST 1990
Rewrote copysign in a more efficient way.
X
Sun Dec 2 06:51:59 EDT 1990
Added infinity().
Amended README to show the current state.
Changed fpumath.h so as to define HUGE = infinity() if you compile (user
code) with IEEE defined.
X
Sun Dec 2 08:26:20 EDT 1990
Modified log functions and inverse hyperbolics for greater accuracy.
X
Sun Dec 2 10:08:47 EDT 1990
Added MASK_ALL1 to fpumath.h. This enables 64 bit precision significands.
Coded sqrtp to ensure that sqrt is rounded to 53 bits - used in
paranoia.c.
Put sqrtp in fpumath.h and in README.
Made very minor changes to pow.s, floor.s and ceil.s.
ieee_retrospective was printing its leading info for a LOS exception -
fixed.
Rewrote finite so that it does not use the floating point unit and so does
not raise any exceptions.
Produced d2dcomb.summ.
X
Sun Dec 2 18:40:47 EDT 1990
Fixed exp, exp10, expm1, exp2, sinh, cosh, tanh to gain more speed.
X
Wed Dec 5 17:25:50 EDT 1990
Coded max_..., min_..., etc.
X
Wed Dec 5 19:38:56 EDT 1990
Fixed some local parameter handling errors in a few of the assembler
files.
X
Wed Dec 5 20:54:37 EDT 1990
Define MAX(MIN)DOUBLE in terms of max(min)_(sub)normal in fpumath.h.
X
Sat Dec 8 06:09:01 EDT 1990
Fixed bug in isnan 8(%ebp) not $8(%ebp) - typo!
Fixed bug in isinf - not doing correct test.
X
Sat Dec 8 06:36:07 EDT 1990
Coded nextafter.
X
Sat Dec 8 06:46:44 EDT 1990
MAXDOUBLE is wrong in math.h it is correct in fpumath.h
SHAR_EOF
chmod 0644 CHANGELOG ||
echo 'restore of CHANGELOG failed'
Wc_c="`wc -c < 'CHANGELOG'`"
test 3702 -eq "$Wc_c" ||
echo 'CHANGELOG: original size 3702, current size' "$Wc_c"
fi
# ============= COPYING ==============
if test -f 'COPYING' -a X"$1" != X"-c"; then
echo 'x - skipping COPYING (File already exists)'
else
echo 'x - extracting COPYING (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'COPYING' &&
X
X GNU GENERAL PUBLIC LICENSE
X Version 1, February 1989
X
X Copyright (C) 1989 Free Software Foundation, Inc.
X 675 Mass Ave, Cambridge, MA 02139, USA
X Everyone is permitted to copy and distribute verbatim copies
X of this license document, but changing it is not allowed.
X
X Preamble
X
X The license agreements of most software companies try to keep users
at the mercy of those companies. By contrast, our General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users. The
General Public License applies to the Free Software Foundation's
software and to any other program whose authors commit to using it.
You can use it for your programs, too.
X
X When we speak of free software, we are referring to freedom, not
price. Specifically, the General Public License is designed to make
sure that you have the freedom to give away or sell copies of free
software, that you receive source code or can get it if you want it,
that you can change the software or use pieces of it in new free
programs; and that you know you can do these things.
X
X To protect your rights, we need to make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.
X
X For example, if you distribute copies of a such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have. You must make sure that they, too, receive or can get the
source code. And you must tell them their rights.
X
X We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.
X
X Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software. If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.
X
X The precise terms and conditions for copying, distribution and
modification follow.
X
X GNU GENERAL PUBLIC LICENSE
X TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
X
X 0. This License Agreement applies to any program or other work which
contains a notice placed by the copyright holder saying it may be
distributed under the terms of this General Public License. The
"Program", below, refers to any such program or work, and a "work based
on the Program" means either the Program or any work containing the
Program or a portion of it, either verbatim or with modifications. Each
licensee is addressed as "you".
X
X 1. You may copy and distribute verbatim copies of the Program's source
code as you receive it, in any medium, provided that you conspicuously and
appropriately publish on each copy an appropriate copyright notice and
disclaimer of warranty; keep intact all the notices that refer to this
General Public License and to the absence of any warranty; and give any
other recipients of the Program a copy of this General Public License
along with the Program. You may charge a fee for the physical act of
transferring a copy.
X
X 2. You may modify your copy or copies of the Program or any portion of
it, and copy and distribute such modifications under the terms of Paragraph
1 above, provided that you also do the following:
X
X a) cause the modified files to carry prominent notices stating that
X you changed the files and the date of any change; and
X
X b) cause the whole of any work that you distribute or publish, that
X in whole or in part contains the Program or any part thereof, either
X with or without modifications, to be licensed at no charge to all
X third parties under the terms of this General Public License (except
X that you may choose to grant warranty protection to some or all
X third parties, at your option).
X
X c) If the modified program normally reads commands interactively when
X run, you must cause it, when started running for such interactive use
X in the simplest and most usual way, to print or display an
X announcement including an appropriate copyright notice and a notice
X that there is no warranty (or else, saying that you provide a
X warranty) and that users may redistribute the program under these
X conditions, and telling the user how to view a copy of this General
X Public License.
X
X d) You may charge a fee for the physical act of transferring a
X copy, and you may at your option offer warranty protection in
X exchange for a fee.
X
Mere aggregation of another independent work with the Program (or its
derivative) on a volume of a storage or distribution medium does not bring
the other work under the scope of these terms.
X
X 3. You may copy and distribute the Program (or a portion or derivative of
it, under Paragraph 2) in object code or executable form under the terms of
Paragraphs 1 and 2 above provided that you also do one of the following:
X
X a) accompany it with the complete corresponding machine-readable
X source code, which must be distributed under the terms of
X Paragraphs 1 and 2 above; or,
X
X b) accompany it with a written offer, valid for at least three
X years, to give any third party free (except for a nominal charge
X for the cost of distribution) a complete machine-readable copy of the
X corresponding source code, to be distributed under the terms of
X Paragraphs 1 and 2 above; or,
X
X c) accompany it with the information you received as to where the
X corresponding source code may be obtained. (This alternative is
X allowed only for noncommercial distribution and only if you
X received the program in object code or executable form alone.)
X
Source code for a work means the preferred form of the work for making
modifications to it. For an executable file, complete source code means
all the source code for all modules it contains; but, as a special
exception, it need not include source code for modules which are standard
libraries that accompany the operating system on which the executable
file runs, or for standard header files or definitions files that
accompany that operating system.
X
X 4. You may not copy, modify, sublicense, distribute or transfer the
Program except as expressly provided under this General Public License.
Any attempt otherwise to copy, modify, sublicense, distribute or transfer
the Program is void, and will automatically terminate your rights to use
the Program under this License. However, parties who have received
copies, or rights to use copies, from you under this General Public
License will not have their licenses terminated so long as such parties
remain in full compliance.
X
X 5. By copying, distributing or modifying the Program (or any work based
on the Program) you indicate your acceptance of this license to do so,
and all its terms and conditions.
X
X 6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the original
licensor to copy, distribute or modify the Program subject to these
terms and conditions. You may not impose any further restrictions on the
recipients' exercise of the rights granted herein.
X
X 7. The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time. Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.
X
Each version is given a distinguishing version number. If the Program
specifies a version number of the license which applies to it and "any
later version", you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation. If the Program does not specify a version number of
the license, you may choose any version ever published by the Free Software
Foundation.
X
X 8. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission. For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this. Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.
X
X NO WARRANTY
X
X 9. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.
X
X 10. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
X
X END OF TERMS AND CONDITIONS
X
X Appendix: How to Apply These Terms to Your New Programs
X
X If you develop a new program, and you want it to be of the greatest
possible use to humanity, the best way to achieve this is to make it
free software which everyone can redistribute and change under these
terms.
X
X To do so, attach the following notices to the program. It is safest to
attach them to the start of each source file to most effectively convey
the exclusion of warranty; and each file should have at least the
"copyright" line and a pointer to where the full notice is found.
X
X <one line to give the program's name and a brief idea of what it does.>
X Copyright (C) 19yy <name of author>
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
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"
fi
# ============= COPYRIGHT ==============
if test -f 'COPYRIGHT' -a X"$1" != X"-c"; then
echo 'x - skipping COPYRIGHT (File already exists)'
else
echo 'x - extracting COPYRIGHT (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'COPYRIGHT' &&
All the assembler source files and the C source files ieee_retro.c and gamma.c
are copyright 1990 to G. Geers.
X
Additional (library) C source and the header mathimpl.h are copyright to
the Regents of the University of California, Berkley.
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 316 -eq "$Wc_c" ||
echo 'COPYRIGHT: original size 316, current size' "$Wc_c"
fi
# ============= Makefile ==============
if test -f 'Makefile' -a X"$1" != X"-c"; then
echo 'x - skipping Makefile (File already exists)'
else
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
LIBDIR=/lib
INCDIR=/usr/include
# if you want to use setinternal in paranoia; define PTEST
#PTEST=
PTEST=-DTEST
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 jn.c ieee_retro.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 setinternal.s setcont.s \
X asinh.s acosh.s atanh.s atan2.s fmod.s ieee_ext.s _getsw.s \
X infinity.s sqrtp.s ieee_values.s nextafter.c
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 jn.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 setinternal.o setcont.o fmod.o ieee_ext.o \
X _getsw.o ieee_retro.o infinity.o sqrtp.o ieee_values.o \
X nextafter.o
X
TESTSRC=paranoia.c
TESTOBJ=paranoia.o
X
all: libfpu test
X
libfpu: $(LIBOBJ)
X ar vr libfpu.a $(LIBOBJ)
X
test: libfpu $(TESTOBJ)
X $(GCC) -o paranoia -DTEST $(TESTOBJ) libfpu.a
X
paranoia: test
X
install: libfpu
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 fpumath.h $(INCDIR)/fpumath.h
X
clean:
X rm -f $(LIBOBJ) $(TESTOBJ) core a.out
X
clobber:
X rm -f $(LIBOBJ) $(TESTOBJ) libfpu.a paranoia core a.out
X
# Dependencies generated using gcc -MM *.c
erf.o : erf.c
gamma.o : gamma.c
j0.o : j0.c mathimpl.h
j1.o : j1.c mathimpl.h
jn.o : jn.c mathimpl.h
lgamma.o : lgamma.c mathimpl.h
X $(CC) -O -c lgamma.c
ieee_retro.o : ieee_retro.c
paranoia.o : paranoia.c
X $(GCC) -c $(PTEST) paranoia.c
SHAR_EOF
chmod 0644 Makefile ||
echo 'restore of Makefile failed'
Wc_c="`wc -c < 'Makefile'`"
test 1984 -eq "$Wc_c" ||
echo 'Makefile: original size 1984, current size' "$Wc_c"
fi
# ============= PROBLEMS ==============
if test -f 'PROBLEMS' -a X"$1" != X"-c"; then
echo 'x - skipping PROBLEMS (File already exists)'
else
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.
SHAR_EOF
chmod 0644 PROBLEMS ||
echo 'restore of PROBLEMS failed'
Wc_c="`wc -c < 'PROBLEMS'`"
test 295 -eq "$Wc_c" ||
echo 'PROBLEMS: original size 295, current size' "$Wc_c"
fi
# ============= README ==============
if test -f 'README' -a X"$1" != X"-c"; then
echo 'x - skipping README (File already exists)'
else
echo 'x - extracting README (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'README' &&
VERSION: This is version 1.1
--------
This is version 1.1 of the alternative 386 math library. It supplants
all previous versions. Patches are available against version 1.0 (no
version numbers in 1.0 or beta releases) by anonymous ftp from
suphys.physics.su.oz.au.
X
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 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.
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
File List
---------
CHANGELOG atan.s exp2.s infinity.s paranoia.c
COPYING atan2.s expm1.s j0.c pow.s
COPYRIGHT atanh.s fabs.s j1.c rint.s
Makefile ceil.s finite.s jn.c scalb.s
PROBLEMS copysign.s floor.s lgamma.c setcont.s
README cos.s fmod.s log.s setinternal.s
TODO cosh.s fpumath.h log10.s sin.s
_getsw.s d2dcomb.summ gamma.c log1p.s sinh.s
acos.s drem.s hypot.s log2.s sqrt.s
acosh.s erf.c ieee_ext.s logb.s sqrtp.s
asin.s exp.s ieee_retro.c mathimpl.h tan.s
asinh.s exp10.s ieee_values.s nextafter.c tanh.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
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 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,
X Berkeley and its contributors''
********************************************************************************
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. 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.
X
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.
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
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@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 9579 -eq "$Wc_c" ||
echo 'README: original size 9579, current size' "$Wc_c"
fi
# ============= TODO ==============
if test -f 'TODO' -a X"$1" != X"-c"; then
echo 'x - skipping TODO (File already exists)'
else
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. Write a float version.
4. Get rid of Berkley code if possible.
SHAR_EOF
chmod 0644 TODO ||
echo 'restore of TODO failed'
Wc_c="`wc -c < 'TODO'`"
test 244 -eq "$Wc_c" ||
echo 'TODO: original size 244, current size' "$Wc_c"
fi
# ============= _getsw.s ==============
if test -f '_getsw.s' -a X"$1" != X"-c"; then
echo 'x - skipping _getsw.s (File already exists)'
else
echo 'x - extracting _getsw.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > '_getsw.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl _getsw
X
_getsw:
X pushl %ebp
X movl %esp,%ebp
X
X fstsw %ax
X
X leave
X ret
SHAR_EOF
chmod 0644 _getsw.s ||
echo 'restore of _getsw.s failed'
Wc_c="`wc -c < '_getsw.s'`"
test 320 -eq "$Wc_c" ||
echo '_getsw.s: original size 320, current size' "$Wc_c"
fi
# ============= acos.s ==============
if test -f 'acos.s' -a X"$1" != X"-c"; then
echo 'x - skipping acos.s (File already exists)'
else
echo 'x - extracting acos.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'acos.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.Lhalfpi:
X .double 1.57079632679489661923
X
X .align 4
X .globl acos
acos:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X
X ftst
X fstsw %ax
X sahf
X jz .Lzero
X
X fst %st(1)
X fmull 8(%ebp)
X fld1
X fsubp
X fsqrt
X fdivp
X fld1
X fpatan
X jnc .Ldone
X
X fldpi
X faddp
X leave
X ret
X
.Lzero:
X fldl .Lhalfpi
X
.Ldone:
X leave
X ret
SHAR_EOF
chmod 0644 acos.s ||
echo 'restore of acos.s failed'
Wc_c="`wc -c < 'acos.s'`"
test 544 -eq "$Wc_c" ||
echo 'acos.s: original size 544, current size' "$Wc_c"
fi
# ============= acosh.s ==============
if test -f 'acosh.s' -a X"$1" != X"-c"; then
echo 'x - skipping acosh.s (File already exists)'
else
echo 'x - extracting acosh.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'acosh.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl acosh
acosh:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X
X fmull 8(%ebp)
X fld1
X fsubrp
X fsqrt
X faddl 8(%ebp)
X fldln2
X fxch %st(1)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 acosh.s ||
echo 'restore of acosh.s failed'
Wc_c="`wc -c < 'acosh.s'`"
test 398 -eq "$Wc_c" ||
echo 'acosh.s: original size 398, current size' "$Wc_c"
fi
# ============= asin.s ==============
if test -f 'asin.s' -a X"$1" != X"-c"; then
echo 'x - skipping asin.s (File already exists)'
else
echo 'x - extracting asin.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'asin.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl asin
asin:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fst %st(1)
X fst %st(2)
X fmulp
X fld1
X fsubp
X fsqrt
X fld1
X fdivp
X fmulp
X fld1
X fpatan
X
X leave
X ret
SHAR_EOF
chmod 0644 asin.s ||
echo 'restore of asin.s failed'
Wc_c="`wc -c < 'asin.s'`"
test 401 -eq "$Wc_c" ||
echo 'asin.s: original size 401, current size' "$Wc_c"
fi
# ============= asinh.s ==============
if test -f 'asinh.s' -a X"$1" != X"-c"; then
echo 'x - skipping asinh.s (File already exists)'
else
echo 'x - extracting asinh.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'asinh.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl asinh
asinh:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X
X fmull 8(%ebp)
X fld1
X faddp
X fsqrt
X faddl 8(%ebp)
X fldln2
X fxch %st(1)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 asinh.s ||
echo 'restore of asinh.s failed'
Wc_c="`wc -c < 'asinh.s'`"
test 397 -eq "$Wc_c" ||
echo 'asinh.s: original size 397, current size' "$Wc_c"
fi
# ============= atan.s ==============
if test -f 'atan.s' -a X"$1" != X"-c"; then
echo 'x - skipping atan.s (File already exists)'
else
echo 'x - extracting atan.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atan.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl atan
atan:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fld1
X fpatan
X
X leave
X ret
SHAR_EOF
chmod 0644 atan.s ||
echo 'restore of atan.s failed'
Wc_c="`wc -c < 'atan.s'`"
test 331 -eq "$Wc_c" ||
echo 'atan.s: original size 331, current size' "$Wc_c"
fi
# ============= atan2.s ==============
if test -f 'atan2.s' -a X"$1" != X"-c"; then
echo 'x - skipping atan2.s (File already exists)'
else
echo 'x - extracting atan2.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atan2.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.Lpi:
X .double 3.14159265358979323846
X
X .align 4
.Lmpi:
X .double -3.14159265358979323846
X
X .align 4
.Lhalfpi:
X .double 1.57079632679489661923
X
X .align 4
.Lmhalfpi:
X .double -1.57079632679489661923
X
X .align 4
X .globl atan2
atan2:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 16(%ebp)
X ftst
X fnstsw %ax
X sahf
X fldl 8(%ebp)
X jz .Lgotzero
X jc .Lgotneg
X
X fdivp
X fld1
X fpatan
X
X leave
X ret
X
.Lgotneg:
X ftst
X fnstsw %ax
X sahf
X jc .Lneg1
X
X fdivp
X fld1
X fpatan
X fldl .Lmpi
X fsubrp
X
X leave
X ret
X
.Lneg1:
X fdivp
X fld1
X fpatan
X fldl .Lpi
X fsubrp
X
X leave
X ret
X
.Lgotzero:
X ftst
X fnstsw %ax
X sahf
X jz .Lzero
X jc .Lneg
X
X fldl .Lhalfpi
X
X leave
X ret
X
.Lzero:
X fldz
X
X leave
X ret
X
.Lneg:
X fldl .Lmhalfpi
X
X leave
X ret
SHAR_EOF
chmod 0644 atan2.s ||
echo 'restore of atan2.s failed'
Wc_c="`wc -c < 'atan2.s'`"
test 931 -eq "$Wc_c" ||
echo 'atan2.s: original size 931, current size' "$Wc_c"
fi
# ============= atanh.s ==============
if test -f 'atanh.s' -a X"$1" != X"-c"; then
echo 'x - skipping atanh.s (File already exists)'
else
echo 'x - extracting atanh.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atanh.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.Lhalf:
X .double 0.500
X
X .align 4
X .globl atanh
atanh:
X pushl %ebp
X movl %esp,%ebp
X
X fld1
X faddl 8(%ebp)
X fld1
X fsubl 8(%ebp)
X fdivrp
X
X fldln2
X fxch %st(1)
X fyl2x
X
X fldl .Lhalf
X fmulp
X
X leave
X ret
SHAR_EOF
chmod 0644 atanh.s ||
echo 'restore of atanh.s failed'
Wc_c="`wc -c < 'atanh.s'`"
test 438 -eq "$Wc_c" ||
echo 'atanh.s: original size 438, current size' "$Wc_c"
fi
# ============= ceil.s ==============
if test -f 'ceil.s' -a X"$1" != X"-c"; then
echo 'x - skipping ceil.s (File already exists)'
else
echo 'x - extracting ceil.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ceil.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Could use ceil(x) = -floor(-x) but this is quicker.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl ceil
ceil:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X fldl 8(%ebp) /* load data */
X
X fstcw -12(%ebp) /* store control word */
X fstcw -16(%ebp) /* store it again */
X orw $0x0800, -16(%ebp) /* round toward +inf */
X fldcw -16(%ebp) /* store new control word */
X frndint /* rounding gives ceil(x) */
X fldcw -12(%ebp) /* restore original control word */
X
X leave
X ret
SHAR_EOF
chmod 0644 ceil.s ||
echo 'restore of ceil.s failed'
Wc_c="`wc -c < 'ceil.s'`"
test 682 -eq "$Wc_c" ||
echo 'ceil.s: original size 682, current size' "$Wc_c"
fi
# ============= copysign.s ==============
if test -f 'copysign.s' -a X"$1" != X"-c"; then
echo 'x - skipping copysign.s (File already exists)'
else
echo 'x - extracting copysign.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'copysign.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.globl copysign
copysign:
X pushl %ebp
X movl %esp,%ebp
X
X movl 20(%ebp), %eax
X andl $0x80000000, %eax
X cmpl $0x80000000, %eax
X je .Lneg
X andl $0x7fffffff, 12(%ebp)
X fldl 8(%ebp) /* Store argument for return */
X leave
X ret
.Lneg:
X orl $0x80000000, 12(%ebp)
X fldl 8(%ebp) /* Store argument for return */
X leave
X ret
SHAR_EOF
chmod 0644 copysign.s ||
echo 'restore of copysign.s failed'
Wc_c="`wc -c < 'copysign.s'`"
test 555 -eq "$Wc_c" ||
echo 'copysign.s: original size 555, current size' "$Wc_c"
fi
# ============= cos.s ==============
if test -f 'cos.s' -a X"$1" != X"-c"; then
echo 'x - skipping cos.s (File already exists)'
else
echo 'x - extracting cos.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cos.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl cos
cos:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fcos
X
X leave
X ret
SHAR_EOF
chmod 0644 cos.s ||
echo 'restore of cos.s failed'
Wc_c="`wc -c < 'cos.s'`"
test 320 -eq "$Wc_c" ||
echo 'cos.s: original size 320, current size' "$Wc_c"
fi
# ============= cosh.s ==============
if test -f 'cosh.s' -a X"$1" != X"-c"; then
echo 'x - skipping cosh.s (File already exists)'
else
echo 'x - extracting cosh.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'cosh.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.LC2:
X .double 0.500
X
X .align 4
X .globl cosh
cosh:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X fst %st(1)
X
X fld1
X fdivp
X faddp
X
X fldl .LC2
X fmulp
X
X leave
X ret
SHAR_EOF
chmod 0644 cosh.s ||
echo 'restore of cosh.s failed'
Wc_c="`wc -c < 'cosh.s'`"
test 486 -eq "$Wc_c" ||
echo 'cosh.s: original size 486, current size' "$Wc_c"
fi
# ============= d2dcomb.summ ==============
if test -f 'd2dcomb.summ' -a X"$1" != X"-c"; then
echo 'x - skipping d2dcomb.summ (File already exists)'
else
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"
fi
# ============= drem.s ==============
if test -f 'drem.s' -a X"$1" != X"-c"; then
echo 'x - skipping drem.s (File already exists)'
else
echo 'x - extracting drem.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'drem.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl drem
drem:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 16(%ebp)
X fldl 8(%ebp)
.Lnotred:
X fprem1
X
X fstsw %ax
X sahf
X jp .Lnotred
X
X leave
X ret
SHAR_EOF
chmod 0644 drem.s ||
echo 'restore of drem.s failed'
Wc_c="`wc -c < 'drem.s'`"
test 381 -eq "$Wc_c" ||
echo 'drem.s: original size 381, current size' "$Wc_c"
fi
# ============= erf.c ==============
if test -f 'erf.c' -a X"$1" != X"-c"; then
echo 'x - skipping erf.c (File already exists)'
else
echo 'x - extracting erf.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'erf.c' &&
/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved. The Berkeley software License Agreement
X * specifies the terms and conditions for redistribution.
X */
X
#ifndef lint
static char sccsid[] = "@(#)erf.c 5.2 (Berkeley) 4/29/88";
#endif /* not lint */
X
/*
X C program for floating point error function
X
X erf(x) returns the error function of its argument
X erfc(x) returns 1.0-erf(x)
X
X erf(x) is defined by
X ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
X
X the entry for erfc is provided because of the
X extreme loss of relative accuracy if erf(x) is
X called for large x and the result subtracted
X from 1. (e.g. for x= 10, 12 places are lost).
X
X There are no error returns.
X
X Calls exp.
X
X Coefficients for large x are #5667 from Hart & Cheney (18.72D).
*/
X
#define M 7
#define N 9
X
extern double sqrt();
extern double exp();
X
static double torp = 1.1283791670955125738961589031;
static double p1[] = {
X 0.804373630960840172832162e5,
X 0.740407142710151470082064e4,
X 0.301782788536507577809226e4,
X 0.380140318123903008244444e2,
X 0.143383842191748205576712e2,
X -.288805137207594084924010e0,
X 0.007547728033418631287834e0,
};
static double q1[] = {
X 0.804373630960840172826266e5,
X 0.342165257924628539769006e5,
X 0.637960017324428279487120e4,
X 0.658070155459240506326937e3,
X 0.380190713951939403753468e2,
X 0.100000000000000000000000e1,
X 0.0,
};
static double p2[] = {
X 0.18263348842295112592168999e4,
X 0.28980293292167655611275846e4,
X 0.2320439590251635247384768711e4,
X 0.1143262070703886173606073338e4,
X 0.3685196154710010637133875746e3,
X 0.7708161730368428609781633646e2,
X 0.9675807882987265400604202961e1,
X 0.5641877825507397413087057563e0,
X 0.0,
};
static double q2[] = {
X 0.18263348842295112595576438e4,
X 0.495882756472114071495438422e4,
X 0.60895424232724435504633068e4,
X 0.4429612803883682726711528526e4,
X 0.2094384367789539593790281779e4,
X 0.6617361207107653469211984771e3,
X 0.1371255960500622202878443578e3,
X 0.1714980943627607849376131193e2,
X 1.0,
};
X
double
erf(arg) double arg;{
X double erfc();
X int sign;
X double argsq;
X double d, n;
X int i;
X
X sign = 1;
X if(arg < 0.){
X arg = -arg;
X sign = -1;
X }
X if(arg < 0.5){
X argsq = arg*arg;
X for(n=0,d=0,i=M-1; i>=0; i--){
X n = n*argsq + p1[i];
X d = d*argsq + q1[i];
X }
X return(sign*torp*arg*n/d);
X }
X if(arg >= 10.)
X return(sign*1.);
X return(sign*(1. - erfc(arg)));
}
X
double
erfc(arg) double arg;{
X double erf();
X double exp();
X double n, d;
X int i;
X
X if(arg < 0.)
X return(2. - erfc(-arg));
/*
X if(arg < 0.5)
X return(1. - erf(arg));
*/
X if(arg >= 10.)
X return(0.);
X
X for(n=0,d=0,i=N-1; i>=0; i--){
X n = n*arg + p2[i];
X d = d*arg + q2[i];
X }
X return(exp(-arg*arg)*n/d);
}
SHAR_EOF
chmod 0644 erf.c ||
echo 'restore of erf.c failed'
Wc_c="`wc -c < 'erf.c'`"
test 2681 -eq "$Wc_c" ||
echo 'erf.c: original size 2681, current size' "$Wc_c"
fi
# ============= exp.s ==============
if test -f 'exp.s' -a X"$1" != X"-c"; then
echo 'x - skipping exp.s (File already exists)'
else
echo 'x - extracting exp.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'exp.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl exp
exp:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 exp.s ||
echo 'restore of exp.s failed'
Wc_c="`wc -c < 'exp.s'`"
test 400 -eq "$Wc_c" ||
echo 'exp.s: original size 400, current size' "$Wc_c"
fi
# ============= exp10.s ==============
if test -f 'exp10.s' -a X"$1" != X"-c"; then
echo 'x - skipping exp10.s (File already exists)'
else
echo 'x - extracting exp10.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'exp10.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl exp10
exp10:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fldl2t
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 exp10.s ||
echo 'restore of exp10.s failed'
Wc_c="`wc -c < 'exp10.s'`"
test 404 -eq "$Wc_c" ||
echo 'exp10.s: original size 404, current size' "$Wc_c"
fi
# ============= exp2.s ==============
if test -f 'exp2.s' -a X"$1" != X"-c"; then
echo 'x - skipping exp2.s (File already exists)'
else
echo 'x - extracting exp2.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'exp2.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl exp2
exp2:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 exp2.s ||
echo 'restore of exp2.s failed'
Wc_c="`wc -c < 'exp2.s'`"
test 387 -eq "$Wc_c" ||
echo 'exp2.s: original size 387, current size' "$Wc_c"
fi
# ============= expm1.s ==============
if test -f 'expm1.s' -a X"$1" != X"-c"; then
echo 'x - skipping expm1.s (File already exists)'
else
echo 'x - extracting expm1.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'expm1.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl expm1
expm1:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X fld1
X fsubrp
X
X leave
X ret
SHAR_EOF
chmod 0644 expm1.s ||
echo 'restore of expm1.s failed'
Wc_c="`wc -c < 'expm1.s'`"
test 418 -eq "$Wc_c" ||
echo 'expm1.s: original size 418, current size' "$Wc_c"
fi
true || echo 'restore of fabs.s failed'
echo End of part 1, continue with part 2
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'