fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 1 (of 1)."
# Contents: README doc funcs funcs/doc funcs/environ funcs/src
# funcs/src/unused funcs/unused tests tests/unused
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:17 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f README -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"README\"
else
echo shar: Extracting \"README\" \(3902 characters\)
sed "s/^X//" >README <<'END_OF_README'
X
XThe original version of pml (Portable Math Library) was written several
Xyears ago, when I was just learning C and trying to switch most of my work
Xfor engineering classes from Fortran to C. I quickly ran into problems
Xwith portability, missing functions, compiler bugs, etc, etc, etc.
XConsidering I was using some pretty primitive stuff at the time, it's a
Xwonder anything worked! Anyway, pml came about from pressures to get
Xsomething that would work in a number of environments and also as a learning
Xtool as to how floating point libraries worked. I don't claim to be a
Xnumerical analyst, then or now, but I did learn a lot writing this stuff
X(including how to write in C).
X
XSeveral times I have started to go back and redo the library to use
Xnewer C features, make it compatible with the Unix version, etc. For
Xexample, the original C compilers I had access to could not pass structures,
Xjust pointers to structures, which made for some pretty ugly code. Changing
Xthe code to pass a 'complex' structure cleaned it up a lot. Given
Xseveral aborted attempts to redo the whole thing (aborted due to lack
Xof motivation and/or time), it's now in somewhat of a mess, and probably
Xinconsistent in several areas. Still, someone might find it useful, or
Xeven be motivated enough to carry on with it, given this code as a base.
X
XSo, in the hopes that someone will 'adopt' pml, and make it their pet,
XI've decided to post my most current code, and relinquish all rights to
Xit. I.E., make it truly public domain, rather than simply freely
Xredistributable.
X
XHints for adoptive parents:
X
X 1. Much of the documentation and internal comments
X is now either wrong, misleading, or both. Beware.
X
X 2. Most C functions use pieces from various iterations
X of my macro based debugging package. For example,
X the LEAVE() macro hasn't been used directly in years...
X You should either chuck this stuff wholesale, or
X convert it to use my latest version (posted on the
X net some months ago).
X
X 3. The implementation of functions like 'inverse complex
X hyperbolic arccosine' may very well be numerically
X naive.
X
X 4. This is your baby now, don't ask me how to change its
X diapers. I wouldn't mind getting pictures now and then
X though. :-) :-)
X
XBelow is the original README file, for reference.
X
X Fred Fish
X 10-Apr-87
X
X======================================================================
X
XPML --- Portable Math Library for C programs.
X
XThis directory contains the PML math library distribution.
XSince it is intended to be a highly portable library, useful
Xon a wide variety of machines, no installation command files
Xare provided. It is assumed that the installer is sufficiently
Xknowledgeable to successfully install the library given the
Xfollowing general guidelines:
X
X(1) The constants in "pml.h" must be suitably defined for
X the environment in which the library is to be used.
X
X(2) The appropriate environment #define variable is
X defined somewhere, either within the preprocessor
X itself, in <stdio.h>, is added to "pml.h", or
X is included on the compiler invocation command
X line.
X
X For example, "#define PDP10" or "#define pdp11".
X
X(3) The file "pmluser.h" is moved to the global header
X file directory, so that an "#include <pmluser.h>"
X will be properly processed.
X
X(4) All of the library source files are compiled.
X
X(5) The furnished test routines are compiled, linked
X with the library routines, and executed.
X
XThis library currently runs essentially unchanged on a wide
Xvariety of machines, under a wide variety of operating systems.
XKnown installations include PDP-11s running RSX-11M and using
Xthe DECUS C system, DECSYSTEM-20 running TOPS-20 using an
XMIT C compiler, Callan Data Systems 68000 system running
XUniplus (Unix port from Unisoft), and an IBM Personal
XComputer running PC-DOS using a CII C compiler. If you
Xknow of any others, please let me know.
X
X Fred Fish
X
END_OF_README
if test 3902 -ne `wc -c <README`; then
echo shar: \"README\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test ! -d doc ; then
echo shar: Creating directory \"doc\"
mkdir doc
fi
if test ! -d funcs ; then
echo shar: Creating directory \"funcs\"
mkdir funcs
fi
if test ! -d funcs/doc ; then
echo shar: Creating directory \"funcs/doc\"
mkdir funcs/doc
fi
if test ! -d funcs/environ ; then
echo shar: Creating directory \"funcs/environ\"
mkdir funcs/environ
fi
if test ! -d funcs/src ; then
echo shar: Creating directory \"funcs/src\"
mkdir funcs/src
fi
if test ! -d funcs/src/unused ; then
echo shar: Creating directory \"funcs/src/unused\"
mkdir funcs/src/unused
fi
if test ! -d funcs/unused ; then
echo shar: Creating directory \"funcs/unused\"
mkdir funcs/unused
fi
if test ! -d tests ; then
echo shar: Creating directory \"tests\"
mkdir tests
fi
if test ! -d tests/unused ; then
echo shar: Creating directory \"tests/unused\"
mkdir tests/unused
fi
echo shar: End of archive 1 \(of 1\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 1 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 2 (of 6)."
# Contents: funcs/environ/testfrexp.c funcs/environ/testldexp.c
# funcs/environ/testmodf.c funcs/src/acos.c funcs/src/acosh.c
# funcs/src/asin.c funcs/src/asinh.c funcs/src/atan2.c
# funcs/src/atanh.c funcs/src/cacos.c funcs/src/cdiv.c
# funcs/src/cmult.c funcs/src/cosh.c funcs/src/csqrt.c
# funcs/src/csubt.c funcs/src/poly.c funcs/src/sinh.c
# funcs/src/tan.c funcs/src/tanh.c funcs/unused/dint.c
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:41 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f funcs/environ/testfrexp.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/environ/testfrexp.c\"
else
echo shar: Extracting \"funcs/environ/testfrexp.c\" \(2017 characters\)
sed "s/^X//" >funcs/environ/testfrexp.c <<'END_OF_funcs/environ/testfrexp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FILE
X *
X * testfrexp.c test the runtime environment function frexp
X *
X * DESCRIPTION
X *
X * This simple minded program is provided to aid in testing
X * the "frexp" function assumed to be provided in the runtime
X * environment. If not provided, a suitable substitute can
X * be coded in C, however the necessary code is very machine
X * dependent and is generally almost trivial to code in assembly
X * language for the specific host machine.
X *
X * The frexp() function takes two arguments, the first is a double
X * and the second is a pointer to an int. Frexp() returns the
X * mantissa of the first arg, and stores the exponent indirectly
X * in the location pointed to by the second arg.
X *
X * See "frexp(3C)" in the Unix System V User's Manual for more
X * information.
X *
X * This program is typically used as:
X *
X * testfrexp <testfrexp.in >junkfile
X * diff testfrexp.out junkfile
X *
X */
X
X#include <stdio.h>
X
Xextern double frexp ();
X
Xmain ()
X{
X double input; /* Input value */
X double dmant; /* Mantissa, 0.5 LE |dmant| LT 1.0 */
X int intexp; /* Exponent as power of 2 */
X
X while (scanf ("%le", &input) == 1) {
X dmant = frexp (input, &intexp);
X printf ("%le %d\n", dmant, intexp);
X }
X}
END_OF_funcs/environ/testfrexp.c
if test 2017 -ne `wc -c <funcs/environ/testfrexp.c`; then
echo shar: \"funcs/environ/testfrexp.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/environ/testldexp.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/environ/testldexp.c\"
else
echo shar: Extracting \"funcs/environ/testldexp.c\" \(2027 characters\)
sed "s/^X//" >funcs/environ/testldexp.c <<'END_OF_funcs/environ/testldexp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FILE
X *
X * testldexp.c test the runtime environment function ldexp
X *
X * DESCRIPTION
X *
X * This simple minded program is provided to aid in testing
X * the "ldexp" function assumed to be provided in the runtime
X * environment. If not provided, a suitable substitute can
X * be coded in C, however the necessary code is very machine
X * dependent and is generally almost trivial to code in assembly
X * language for the specific host machine.
X *
X * The ldexp() function takes two arguments, the first is a double
X * which is the mantissa of the returned double value. The second
X * is an int which is the exponent of the returned double. The
X * result is (mantissa * (2 ** exp)).
X *
X * See "frexp(3C)" in the Unix System V User's Manual for more
X * information.
X *
X * This program is typically used as:
X *
X * testldexp <testldexp.in >junkfile
X * diff testldexp.out junkfile
X *
X */
X
X#include <stdio.h>
X
Xextern double ldexp ();
X
Xmain ()
X{
X double dmant; /* Mantissa, 0.5 LE |dmant| LT 1.0 */
X int intexp; /* Exponent as power of 2 */
X double result; /* dmant times 2 to the intexp */
X
X while (scanf ("%le %d", &dmant, &intexp) == 2) {
X result = ldexp (dmant, intexp);
X printf ("%le\n", result);
X }
X}
END_OF_funcs/environ/testldexp.c
if test 2027 -ne `wc -c <funcs/environ/testldexp.c`; then
echo shar: \"funcs/environ/testldexp.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/environ/testmodf.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/environ/testmodf.c\"
else
echo shar: Extracting \"funcs/environ/testmodf.c\" \(1977 characters\)
sed "s/^X//" >funcs/environ/testmodf.c <<'END_OF_funcs/environ/testmodf.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FILE
X *
X * testmodf.c test the runtime environment function modf
X *
X * This simple minded program is provided to aid in testing
X * the "modf" function assumed to be provided in the runtime
X * environment. If not provided, a suitable substitute can
X * be coded in C, however the necessary code is very machine
X * dependent and is generally almost trivial to code in assembly
X * language for the specific host machine.
X *
X * The modf() function takes two arguments. The first is a double value
X * and the second is a pointer to a double. Modf() returns the
X * signed fraction part of the first argument, and stores the integral
X * part (as a double) indirectly in the location pointed to by the
X * second argument. Note that both the direct and indirect result will
X * have the same sign, and:
X *
X * <integral part> + <fractional part> = <original value>
X *
X * See "frexp(3C)" in the Unix System V User's Manual for more
X * information.
X *
X */
X
Xextern double modf ();
X
Xmain ()
X{
X double input; /* Input value */
X double frac;
X double ipart;
X
X while (scanf ("%le", &input) == 1) {
X frac = modf (input, &ipart);
X printf ("%le %le\n", frac, ipart);
X }
X}
END_OF_funcs/environ/testmodf.c
if test 1977 -ne `wc -c <funcs/environ/testmodf.c`; then
echo shar: \"funcs/environ/testmodf.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/acos.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/acos.c\"
else
echo shar: Extracting \"funcs/src/acos.c\" \(2942 characters\)
sed "s/^X//" >funcs/src/acos.c <<'END_OF_funcs/src/acos.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * acos double precision arc cosine
X *
X * KEY WORDS
X *
X * acos
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision arc cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double acos (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-1.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * in atan is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes arccosine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then call
X * matherr and return 0.0 by default.
X *
X * (2) If x = 0.0 then acos(x) = PI/2.
X *
X * (3) If x = 1.0 then acos(x) = 0.0
X *
X * (4) If x = -1.0 then acos(x) = PI.
X *
X * (5) If 0.0 < x < 1.0 then
X * acos(x) = atan(Y)
X * Y = sqrt[1-(x**2)] / x
X *
X * (4) If -1.0 < x < 0.0 then
X * acos(x) = atan(Y) + PI
X * Y = sqrt[1-(x**2)] / x
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "acos";
X
X
Xdouble acos (x)
Xdouble x;
X{
X double y;
X extern double atan();
X extern double sqrt();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("acosin", "arg %le", x);
X if ( x > 1.0 || x < -1.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X } else if (x == 0.0) {
X xcpt.retval = HALFPI;
X } else if (x == 1.0) {
X xcpt.retval = 0.0;
X } else if (x == -1.0) {
X xcpt.retval = PI;
X } else {
X y = atan ( sqrt (1.0 - (x * x)) / x );
X if (x > 0.0) {
X xcpt.retval = y;
X } else {
X xcpt.retval = y + PI;
X }
X }
X DBUG_3 ("acosout", "result %le", x);
X DBUG_RETURN (x);
X}
END_OF_funcs/src/acos.c
if test 2942 -ne `wc -c <funcs/src/acos.c`; then
echo shar: \"funcs/src/acos.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/acosh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/acosh.c\"
else
echo shar: Extracting \"funcs/src/acosh.c\" \(2642 characters\)
sed "s/^X//" >funcs/src/acosh.c <<'END_OF_funcs/src/acosh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * acosh double precision hyperbolic arc cosine
X *
X * KEY WORDS
X *
X * acosh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double acosh (x)
X * double x;
X *
X * RESTRICTIONS
X *
X * The range of the ACOSH function is all real numbers greater
X * than or equal to 1.0 however large arguments may cause
X * overflow in the x squared portion of the function evaluation.
X *
X * For precision information refer to documentation of the
X * floating point library primatives called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes acosh(x) from:
X *
X * 1. If x < 1.0 then report illegal
X * argument and return zero.
X *
X * 2. If x > sqrt(MAXDOUBLE) then
X * set x = sqrt(MAXDOUBLE and
X * continue after reporting overflow.
X *
X * 3. acosh(x) = log [x+sqrt(x**2 - 1)]
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "acosh";
X
X
Xdouble acosh (x)
Xdouble x;
X{
X auto struct exception xcpt;
X extern double log ();
X extern double sqrt ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("acoshin", "arg %le", x);
X if (x < 1.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = ERANGE;
X xcpt.retval = 0.0;
X }
X } else if (x > SQRT_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X x = SQRT_MAXDOUBLE;
X xcpt.retval = log (2* SQRT_MAXDOUBLE);
X }
X } else {
X xcpt.retval = log (x + sqrt (x*x - 1.0));
X }
X DBUG_3 ("acoshout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
X
END_OF_funcs/src/acosh.c
if test 2642 -ne `wc -c <funcs/src/acosh.c`; then
echo shar: \"funcs/src/acosh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/asin.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/asin.c\"
else
echo shar: Extracting \"funcs/src/asin.c\" \(2675 characters\)
sed "s/^X//" >funcs/src/asin.c <<'END_OF_funcs/src/asin.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * asin double precision arc sine
X *
X * KEY WORDS
X *
X * asin
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision arc sine of double precision
X * floating point argument.
X *
X * If argument is less than -1.0 or greater than +1.0, calls
X * matherr with a DOMAIN error. If matherr does not handle
X * the error then prints error message and returns 0.
X *
X * USAGE
X *
X * double asin (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-2.
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the floating
X * point library primatives called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes arcsine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then
X * call matherr and return 0.0 by default.
X *
X * (2) If x = 0.0 then asin(x) = 0.0
X *
X * (3) If x = 1.0 then asin(x) = PI/2.
X *
X * (4) If x = -1.0 then asin(x) = -PI/2
X *
X * (5) If -1.0 < x < 1.0 then
X * asin(x) = atan(y)
X * y = x / sqrt[1-(x**2)]
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "asin";
X
X
Xdouble asin (x)
Xdouble x;
X{
X extern double atan ();
X extern double sqrt ();
X struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("asinin", "arg %le", x);
X if ( x > 1.0 || x < -1.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X } else if (x == 0.0) {
X xcpt.retval = 0.0;
X } else if (x == 1.0) {
X xcpt.retval = HALFPI;
X } else if (x == -1.0) {
X xcpt.retval = -HALFPI;
X } else {
X xcpt.retval = atan ( x / sqrt (1.0 - (x * x)) );
X }
X DBUG_3 ("asinout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/asin.c
if test 2675 -ne `wc -c <funcs/src/asin.c`; then
echo shar: \"funcs/src/asin.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/asinh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/asinh.c\"
else
echo shar: Extracting \"funcs/src/asinh.c\" \(2315 characters\)
sed "s/^X//" >funcs/src/asinh.c <<'END_OF_funcs/src/asinh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * asinh double precision hyperbolic arc sine
X *
X * KEY WORDS
X *
X * asinh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc sine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double asinh (x)
X * double x;
X *
X * RESTRICTIONS
X *
X * The domain of the ASINH function is the entire real axis
X * however the evaluation of x squared may cause overflow
X * for large magnitudes.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes asinh(x) from:
X *
X * 1. Let xmax = sqrt(MAXDOUBLE - 1)
X *
X * 2. If x < -xmax or xmax < x then
X * let x = xmax and flag overflow.
X *
X * 3. asinh(x) = log [x+sqrt(x**2 + 1)]
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "asinh";
X
X
Xdouble asinh (x)
Xdouble x;
X{
X auto struct exception xcpt;
X extern double log ();
X extern double sqrt ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("asinhin", "arg %le", x);
X if (x < -SQRT_MAXDOUBLE || x > SQRT_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = log (2 * SQRT_MAXDOUBLE);
X }
X } else {
X xcpt.retval = log (x + sqrt(x*x + 1.0));
X }
X DBUG_3 ("asinhout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/asinh.c
if test 2315 -ne `wc -c <funcs/src/asinh.c`; then
echo shar: \"funcs/src/asinh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/atan2.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/atan2.c\"
else
echo shar: Extracting \"funcs/src/atan2.c\" \(2363 characters\)
sed "s/^X//" >funcs/src/atan2.c <<'END_OF_funcs/src/atan2.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * atan2 double precision arc tangent of two arguments
X *
X * KEY WORDS
X *
X * atan2
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision arc tangent of two
X * double precision floating point arguments ( atan(Y/X) ).
X *
X * USAGE
X *
X * double atan2(x,y)
X * double x;
X * double y;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-4.
X *
X * RESTRICTIONS
X *
X * Note that the argument usage is exactly the reverse of the
X * common FORTRAN usage where atan2(x,y) computes atan(x/y).
X * The usage here is less confusing than remembering that x is
X * really y and y is really x.
X *
X * For precision information refer to documentation of the
X * other floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X * INTERNALS
X *
X * Computes atan(y/x) from:
X *
X * 1. If x = 0 then
X * atan(x,y) = PI/2 * (sign(y))
X *
X * 2. If x > 0 then
X * atan(x,y) = atan(y/x)
X *
X * 3. If x < 0 then atan2(x,y) =
X * PI*(sign(y)) + atan(y/x)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble atan2 (x, y)
Xdouble x;
Xdouble y;
X{
X double result;
X extern double sign();
X extern double atan();
X
X ENTER ("atan2");
X DEBUG4 ("atan2in", "x = %le y = %le", x, y);
X if (x == 0.0) {
X result = sign (HALFPI, y);
X } else if (x > 0.0) {
X result = atan (y/x);
X } else {
X result = atan (y/x) + sign (PI, y);
X }
X DEBUG3 ("atan2out", "result %le", result);
X LEAVE ();
X return (result);
X}
END_OF_funcs/src/atan2.c
if test 2363 -ne `wc -c <funcs/src/atan2.c`; then
echo shar: \"funcs/src/atan2.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/atanh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/atanh.c\"
else
echo shar: Extracting \"funcs/src/atanh.c\" \(2301 characters\)
sed "s/^X//" >funcs/src/atanh.c <<'END_OF_funcs/src/atanh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * atanh double precision hyperbolic arc tangent
X *
X * KEY WORDS
X *
X * atanh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc tangent of double precision
X * floating point number.
X *
X * USAGE
X *
X * double atanh (x)
X * double x;
X *
X * RESTRICTIONS
X *
X * The range of the atanh function is -1.0 to +1.0 exclusive.
X * Certain pathological cases near 1.0 may cause overflow
X * in evaluation of 1+x / 1-x, depending upon machine exponent
X * range and mantissa precision.
X *
X * For precision information refer to documentation of the
X * other floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes atanh(x) from:
X *
X * 1. If x <= -1.0 or x >= 1.0
X * then report argument out of
X * range and return 0.0
X *
X * 2. atanh(x) = 0.5 * log((1+x)/(1-x))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "atanh";
X
X
Xdouble atanh (x)
Xdouble x;
X{
X auto struct exception xcpt;
X extern double log ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("atanhin", "arg %le", x);
X if (x <= -1.0 || x >= 1.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = ERANGE;
X xcpt.retval = 0.0;
X }
X } else {
X xcpt.retval = 0.5 * log ((1+x) / (1-x));
X }
X DBUG_3 ("atanhout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/atanh.c
if test 2301 -ne `wc -c <funcs/src/atanh.c`; then
echo shar: \"funcs/src/atanh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cacos.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cacos.c\"
else
echo shar: Extracting \"funcs/src/cacos.c\" \(1983 characters\)
sed "s/^X//" >funcs/src/cacos.c <<'END_OF_funcs/src/cacos.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cacos complex double precision arc cosine
X *
X * KEY WORDS
X *
X * cacos
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex arc cosine of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX cacos (z)
X * COMPLEX z;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe,Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex arc cosine of z = x + jy from:
X *
X * cacos(z) = -j * clog(z + j * csqrt(1-z*z))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX cacos (z)
XCOMPLEX z;
X{
X COMPLEX temp;
X double swaptemp;
X extern COMPLEX cmult (), csqrt (), clog ();
X
X ENTER ("cacos");
X DEBUG4 ("cacosin", "arg %le %le", z.real, z.imag);
X temp = cmult(z, z);
X temp.real = 1.0 - temp.real;
X temp.imag = -temp.imag;
X temp = csqrt (temp);
X swaptemp = temp.real;
X temp.real = -temp.imag;
X temp.imag = swaptemp;
X temp.real += z.real;
X temp.imag += z.imag;
X temp = clog (temp);
X z.real = temp.imag;
X z.imag = -temp.real;
X DEBUG4 ("cacosout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
X
END_OF_funcs/src/cacos.c
if test 1983 -ne `wc -c <funcs/src/cacos.c`; then
echo shar: \"funcs/src/cacos.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cdiv.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cdiv.c\"
else
echo shar: Extracting \"funcs/src/cdiv.c\" \(2460 characters\)
sed "s/^X//" >funcs/src/cdiv.c <<'END_OF_funcs/src/cdiv.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cdiv double precision complex division
X *
X * KEY WORDS
X *
X * cdiv
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex result of division of
X * first double precision complex argument by second double
X * precision complex argument.
X *
X * USAGE
X *
X * COMPLEX cdiv (numerator, denominator)
X * COMPLEX numerator;
X * COMPLEX denominator;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes cdiv(znum,zden) from:
X *
X * 1. Let znum = a + j b
X * Let zden = c + j d
X *
X * 2. denom = c*c + d*d
X *
X * 3. If denom is zero then log error,
X * set r_cdiv = maximum floating value,
X * i_cdiv = 0, and go to step 5.
X *
X * 4. r_cdiv = (a*c + b*d) / denom
X * i_cdiv = (c*b - a*d) / denom
X *
X * 5. Then cdiv(znum,zden) = r_cdiv + j i_cdiv
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX cdiv (znum, zden)
XCOMPLEX znum;
XCOMPLEX zden;
X{
X COMPLEX result;
X double denom;
X
X ENTER ("cdiv");
X DEBUG4 ("cdivin", "arg1 %le %le", znum.real, znum.imag);
X DEBUG4 ("cdivin", "arg2 %le %le", zden.real, zden.imag);
X denom = (zden.real * zden.real) + (zden.imag * zden.imag);
X if (denom == 0.0) {
X/****
X pmlerr(C_DIV_ZERO);
X result.real = MAX_POS_DBLF;
X******/
X result.real = 0.0; /* TERRIBLY WRONG! */
X result.imag = 0.0;
X } else {
X result.real = ((znum.real*zden.real)+(znum.imag*zden.imag)) / denom;
X result.imag = ((zden.real*znum.imag)-(znum.real*zden.imag)) / denom;
X }
X DEBUG4 ("cdivout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
END_OF_funcs/src/cdiv.c
if test 2460 -ne `wc -c <funcs/src/cdiv.c`; then
echo shar: \"funcs/src/cdiv.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cmult.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cmult.c\"
else
echo shar: Extracting \"funcs/src/cmult.c\" \(1994 characters\)
sed "s/^X//" >funcs/src/cmult.c <<'END_OF_funcs/src/cmult.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cmult double precision complex multiplication
X *
X * KEY WORDS
X *
X * cmult
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex result of multiplication of
X * first double precision complex argument by second double
X * precision complex argument.
X *
X * USAGE
X *
X * COMPLEX cmult (z1, z2)
X * COMPLEX z1;
X * COMPLEX z2;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes cmult(z1,z2) from:
X *
X * 1. Let z1 = a + j b
X * Let z2 = c + j d
X *
X * 2. r_cmult = (a*c - b*d)
X * i_cmult = (a*d + c*b)
X *
X * 3. Then cmult(z1,z2) = r_cmult + j i_cmult
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX cmult (z1, z2)
XCOMPLEX z1;
XCOMPLEX z2;
X{
X COMPLEX result;
X
X ENTER ("cmult");
X DEBUG4 ("cmultin", "arg1 %le %le", z1.real, z1.imag);
X DEBUG4 ("cmultin", "arg2 %le %le", z2.real, z2.imag);
X result.real = (z1.real * z2.real) - (z1.imag * z2.imag);
X result.imag = (z1.real * z2.imag) + (z2.real * z1.imag);
X DEBUG4 ("cmultout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
END_OF_funcs/src/cmult.c
if test 1994 -ne `wc -c <funcs/src/cmult.c`; then
echo shar: \"funcs/src/cmult.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cosh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cosh.c\"
else
echo shar: Extracting \"funcs/src/cosh.c\" \(2424 characters\)
sed "s/^X//" >funcs/src/cosh.c <<'END_OF_funcs/src/cosh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cosh double precision hyperbolic cosine
X *
X * KEY WORDS
X *
X * cosh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double cosh (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-4
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes hyperbolic cosine from:
X *
X * cosh(X) = 0.5 * (exp(X) + exp(-X))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "cosh";
X
X
Xdouble cosh (x)
Xdouble x;
X{
X auto struct exception xcpt;
X extern double exp ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("coshin", "arg %le", x);
X if (x > LOGE_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MAXDOUBLE;
X }
X } else if (x < LOGE_MINDOUBLE) {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: UNDERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MINDOUBLE;
X }
X } else {
X x = exp (x);
X xcpt.retval = 0.5 * (x + 1.0/x);
X }
X DBUG_3 ("coshout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/cosh.c
if test 2424 -ne `wc -c <funcs/src/cosh.c`; then
echo shar: \"funcs/src/cosh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/csqrt.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/csqrt.c\"
else
echo shar: Extracting \"funcs/src/csqrt.c\" \(2679 characters\)
sed "s/^X//" >funcs/src/csqrt.c <<'END_OF_funcs/src/csqrt.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * csqrt complex double precision square root
X *
X * KEY WORDS
X *
X * csqrt
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex square root of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX csqrt (z)
X * COMPLEX z;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-12
X *
X * RESTRICTIONS
X *
X * The relative error in the double precision square root
X * computation is 10**(-30.1) after three applications
X * of Heron's iteration for the square root.
X *
X * However, this assumes exact arithmetic in the iterations
X * and initial approximation. Additional errors may occur
X * due to truncation, rounding, or machine precision limits.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex square root of z = x + j y from:
X *
X * 1. If z = 0 + j 0 then return z.
X *
X * 2. root = sqrt((dabs(x) + cabs(z)) / 2)
X *
X * 3. q = y / (2 * root)
X *
X * 4. If x >= 0 then
X * csqrt(z) = (root,q)
X *
X * 5. If x < 0 and y >= 0 then
X * csqrt(z) = (q,root)
X *
X * 6. If x < 0 and y < 0 then
X * csqrt(z) = (-q,root)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX csqrt (z)
XCOMPLEX z;
X{
X double root, q;
X extern double dabs(), sqrt(), cabs ();
X
X ENTER ("csqrt");
X DEBUG4 ("csqrtin", "arg %le %le", z.real, z.imag);
X if (z.real != 0.0 || z.imag != 0.0) {
X root = sqrt (0.5 * (dabs (z.real) + cabs (z)));
X q = z.imag / (2.0 * root);
X if (z.real >= 0.0) {
X z.real = root;
X z.imag = q;
X } else if (z.imag < 0.0) {
X z.real = -q;
X z.imag = -root;
X } else {
X z.real = q;
X z.imag = root;
X }
X }
X DEBUG4 ("csqrtout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
END_OF_funcs/src/csqrt.c
if test 2679 -ne `wc -c <funcs/src/csqrt.c`; then
echo shar: \"funcs/src/csqrt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/csubt.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/csubt.c\"
else
echo shar: Extracting \"funcs/src/csubt.c\" \(2005 characters\)
sed "s/^X//" >funcs/src/csubt.c <<'END_OF_funcs/src/csubt.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * csubt double precision complex subtraction
X *
X * KEY WORDS
X *
X * csubt
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex result of subtraction of
X * second double precision complex argument from first double
X * precision complex argument.
X *
X * Note that the complex subtraction function is so simple that
X * it would not normally be called as a function but simply
X * done "inline". It is supplied mostly for completeness.
X *
X * USAGE
X *
X * COMPLEX csubt (z1, z2)
X * COMPLEX z1;
X * COMPLEX z2;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes csubt (z1, z2) from:
X *
X * 1. Let z1 = a + j b
X * Let z2 = c + j d
X *
X * 2. Then csubt(z1,z2) = (a - c) + j (b - d)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX csubt (z1, z2)
XCOMPLEX z1;
XCOMPLEX z2;
X{
X ENTER ("csubt");
X DEBUG4 ("csubtin", "arg %le %le", z1.real, z1.imag);
X DEBUG4 ("csubtin", "arg2 %le %le", z2.real, z2.imag);
X z1.real -= z2.real;
X z1.imag -= z2.imag;
X DEBUG4 ("csubtout", "result %le %le", z1.real, z1.imag);
X LEAVE ();
X return (z1);
X}
END_OF_funcs/src/csubt.c
if test 2005 -ne `wc -c <funcs/src/csubt.c`; then
echo shar: \"funcs/src/csubt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/poly.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/poly.c\"
else
echo shar: Extracting \"funcs/src/poly.c\" \(2049 characters\)
sed "s/^X//" >funcs/src/poly.c <<'END_OF_funcs/src/poly.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * poly double precision polynomial evaluation
X *
X * KEY WORDS
X *
X * poly
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Evaluates a polynomial and returns double precision
X * result. Is passed a the order of the polynomial,
X * a pointer to an array of double precision polynomial
X * coefficients (in ascending order), and the independent
X * variable.
X *
X * USAGE
X *
X * double poly (order, coeffs, x)
X * int order;
X * double *coeffs;
X * double x;
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Evalates the polynomial using recursion and the form:
X *
X * P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble poly (order, coeffs, x)
Xregister int order;
Xdouble *coeffs;
Xdouble x;
X{
X auto double curr_coeff;
X auto double rtn_value;
X
X DBUG_ENTER ("poly");
X DBUG_5 ("polyin", "args %d %#x %le", order, coeffs, x);
X if (order <= 0) {
X rtn_value = *coeffs;
X } else {
X curr_coeff = *coeffs; /* Bug in Unisoft's compiler. Does not */
X coeffs++; /* generate good code for *coeffs++ */
X rtn_value = curr_coeff + x * poly (--order, coeffs, x);
X }
X DBUG_3 ("polyout", "result %le", rtn_value);
X DBUG_RETURN (rtn_value);
X}
END_OF_funcs/src/poly.c
if test 2049 -ne `wc -c <funcs/src/poly.c`; then
echo shar: \"funcs/src/poly.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/sinh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/sinh.c\"
else
echo shar: Extracting \"funcs/src/sinh.c\" \(2423 characters\)
sed "s/^X//" >funcs/src/sinh.c <<'END_OF_funcs/src/sinh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sinh double precision hyperbolic sine
X *
X * KEY WORDS
X *
X * sinh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic sine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double sinh (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X * RESTRICTIONS
X *
X * Inputs greater than ln(MAXDOUBLE) result in overflow.
X * Inputs less than ln(MINDOUBLE) result in underflow.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes hyperbolic sine from:
X *
X * sinh(x) = 0.5 * (exp(x) - 1.0/exp(x))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "sinh";
X
X
Xdouble sinh (x)
Xdouble x;
X{
X extern double exp ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("sinhin", "arg %le", x);
X if (x > LOGE_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MAXDOUBLE;
X }
X } else if (x < LOGE_MINDOUBLE) {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: UNDERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MINDOUBLE;
X }
X } else {
X x = exp (x);
X xcpt.retval = 0.5 * (x - (1.0 / x));
X }
X DBUG_3 ("sinhout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/sinh.c
if test 2423 -ne `wc -c <funcs/src/sinh.c`; then
echo shar: \"funcs/src/sinh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/tan.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/tan.c\"
else
echo shar: Extracting \"funcs/src/tan.c\" \(2208 characters\)
sed "s/^X//" >funcs/src/tan.c <<'END_OF_funcs/src/tan.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * tan Double precision tangent
X *
X * KEY WORDS
X *
X * tan
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns tangent of double precision floating point number.
X *
X * USAGE
X *
X * double tan (x)
X * double x;
X *
X * INTERNALS
X *
X * Computes the tangent from tan(x) = sin(x) / cos(x).
X *
X * If cos(x) = 0 and sin(x) >= 0, then returns largest
X * floating point number on host machine.
X *
X * If cos(x) = 0 and sin(x) < 0, then returns smallest
X * floating point number on host machine.
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp. B-8
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "tan";
X
Xdouble tan (x)
Xdouble x;
X{
X double sinx;
X double cosx;
X auto struct exception xcpt;
X extern double sin ();
X extern double cos();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("tanin", "arg %le", x);
X sinx = sin (x);
X cosx = cos (x);
X if (cosx == 0.0) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X if (sinx >= 0.0) {
X xcpt.retval = MAXDOUBLE;
X } else {
X xcpt.retval = -MAXDOUBLE;
X }
X }
X } else {
X xcpt.retval = sinx / cosx;
X }
X DBUG_3 ("tanout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/tan.c
if test 2208 -ne `wc -c <funcs/src/tan.c`; then
echo shar: \"funcs/src/tan.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/tanh.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/tanh.c\"
else
echo shar: Extracting \"funcs/src/tanh.c\" \(2312 characters\)
sed "s/^X//" >funcs/src/tanh.c <<'END_OF_funcs/src/tanh.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * tanh double precision hyperbolic tangent
X *
X * KEY WORDS
X *
X * tanh
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic tangent of double precision
X * floating point number.
X *
X * USAGE
X *
X * double tanh (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes hyperbolic tangent from:
X *
X * tanh(x) = 1.0 for x > TANH_MAXARG
X * = -1.0 for x < -TANH_MAXARG
X * = sinh(x) / cosh(x) otherwise
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "tanh";
X
X
Xdouble tanh (x)
Xdouble x;
X{
X auto struct exception xcpt;
X register int positive;
X extern double sinh ();
X extern double cosh ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("tanhin", "arg %le", x);
X if (x > TANH_MAXARG || x < -TANH_MAXARG) {
X if (x > 0.0) {
X positive = 1;
X } else {
X positive = 0;
X }
X xcpt.type = PLOSS;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: PLOSS error\n", funcname);
X errno = ERANGE;
X if (positive) {
X xcpt.retval = 1.0;
X } else {
X xcpt.retval = -1.0;
X }
X }
X } else {
X xcpt.retval = sinh (x) / cosh (x);
X }
X DBUG_3 ("tanhout", "result %le", xcpt.retval);
X return (xcpt.retval);
X}
END_OF_funcs/src/tanh.c
if test 2312 -ne `wc -c <funcs/src/tanh.c`; then
echo shar: \"funcs/src/tanh.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/dint.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/dint.c\"
else
echo shar: Extracting \"funcs/unused/dint.c\" \(2630 characters\)
sed "s/^X//" >funcs/unused/dint.c <<'END_OF_funcs/unused/dint.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * dint double precision integer portion
X *
X * KEY WORDS
X *
X * dint
X * machine dependent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns integer portion of double precision number as double
X * precision number.
X *
X * USAGE
X *
X * double dint(x)
X * double x;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * RESTRICTIONS
X *
X * The current DEC-20 C system treats double as float. This
X * routine will need to be modified when true double precision
X * is implemented.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
X#ifdef PDP10
X
X#define W1_FBITS 27 /* Number of fractional bits in word 1 */
X#define WORD_MASK 0777777777777 /* Mask for all 36 bits of first word */
X
Xdouble dint(x)
Xdouble x;
X{
X register int *vpntr, exponent;
X int xexp();
X
X vpntr = &x;
X if ((exponent=xexp(x)) <= 0) {
X x = 0.0;
X } else if (exponent <= W1_FBITS) {
X *vpntr &= (WORD_MASK << (W1_FBITS - exponent));
X } else {
X pmlerr(DINT_2BIG);
X x = 0.0;
X }
X if (x < 0.0) {
X x += 1.0;
X }
X return (x);
X}
X
X#else
X
X#define W1_FBITS 24 /* (NOTE HIDDEN BIT NORMALIZATION) */
X#define W2_FBITS 32 /* Number of fractional bits in word 2 */
X#define WORD_MASK 0xFFFFFFFF /* Mask for each long word of double */
X
Xstatic union {
X double dval;
X long lval[2];
X} share;
X
Xdouble dint(x)
Xdouble x;
X{
X int exponent, xexp(), fbitdown;
X register long *lpntr;
X
X lpntr = &share.lval[0];
X share.dval = x;
X if ((exponent=xexp(x)) <= 0) {
X share.dval = 0.0;
X } else if (exponent <= W1_FBITS) {
X *lpntr++ &= (WORD_MASK << (W1_FBITS - exponent));
X *lpntr++ = 0;
X } else if (exponent <= (fbitdown = W1_FBITS+W2_FBITS)) {
X lpntr++;
X *lpntr++ &= (WORD_MASK << (fbitdown - exponent));
X } else {
X pmlerr(DINT_2BIG);
X }
X return (share.dval);
X}
X#endif
END_OF_funcs/unused/dint.c
if test 2630 -ne `wc -c <funcs/unused/dint.c`; then
echo shar: \"funcs/unused/dint.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 2 \(of 6\).
cp /dev/null ark2isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 3 (of 6)."
# Contents: funcs/src/atan.c funcs/src/cos.c funcs/src/log.c
# funcs/src/pml.h funcs/src/sin.c funcs/src/sqrt.c
# funcs/unused/scale.c funcs/unused/xexp.c funcs/unused/xmant.c
# tests/c2c.dat tests/cc2c.dat
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:43 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f funcs/src/atan.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/atan.c\"
else
echo shar: Extracting \"funcs/src/atan.c\" \(5187 characters\)
sed "s/^X//" >funcs/src/atan.c <<'END_OF_funcs/src/atan.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * atan double precision arc tangent
X *
X * KEY WORDS
X *
X * atan
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision arc tangent of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double atan (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 120-130.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes arctangent(x) from:
X *
X * (1) If x < 0 then negate x, perform steps
X * 2, 3, and 4, and negate the returned
X * result. This makes use of the identity
X * atan(-x) = -atan(x).
X *
X * (2) If argument x > 1 at this point then
X * test to be sure that x can be inverted
X * without underflowing. If not, reduce
X * x to largest possible number that can
X * be inverted, issue warning, and continue.
X * Perform steps 3 and 4 with arg = 1/x
X * and subtract returned result from
X * pi/2. This makes use of the identity
X * atan(x) = pi/2 - atan(1/x) for x>0.
X *
X * (3) At this point 0 <= x <= 1. If
X * x > tan(pi/12) then perform step 4
X * with arg = (x*sqrt(3)-1)/(sqrt(3)+x)
X * and add pi/6 to returned result. This
X * last transformation maps arguments
X * greater than tan(pi/12) to arguments
X * in range 0...tan(pi/12).
X *
X * (4) At this point the argument has been
X * transformed so that it lies in the
X * range -tan(pi/12)...tan(pi/12).
X * Since very small arguments may cause
X * underflow in the polynomial evaluation,
X * a final check is performed. If the
X * argument is less than the underflow
X * bound, the function returns x, since
X * atan(x) approaches asin(x) which
X * approaches x, as x goes to zero.
X *
X * (5) atan(x) is now computed by one of the
X * approximations given in the cited
X * reference (Hart). That is:
X *
X * atan(x) = x * SUM [ C[i] * x**(2*i) ]
X * over i = {0,1,...8}.
X *
X * Where:
X *
X * C[0] = .9999999999999999849899
X * C[1] = -.333333333333299308717
X * C[2] = .1999999999872944792
X * C[3] = -.142857141028255452
X * C[4] = .11111097898051048
X * C[5] = -.0909037114191074
X * C[6] = .0767936869066
X * C[7] = -.06483193510303
X * C[8] = .0443895157187
X * (coefficients from HART table #4945 pg 267)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double atan_coeffs[] = {
X .9999999999999999849899, /* P0 must be first */
X -.333333333333299308717,
X .1999999999872944792,
X -.142857141028255452,
X .11111097898051048,
X -.0909037114191074,
X .0767936869066,
X -.06483193510303,
X .0443895157187 /* Pn must be last */
X};
X
Xstatic char funcname[] = "atan";
X
X#define LAST_BOUND 0.2679491924311227074725 /* tan (PI/12) */
X
X
Xdouble atan (x)
Xdouble x;
X{
X register int order;
X double xt2;
X double t1;
X double t2;
X extern double poly ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("atanin", "arg %le", x);
X if (x < 0.0) {
X xcpt.retval = -(atan (-x));
X } else if (x > 1.0) {
X if (x < MAXDOUBLE && x > -MAXDOUBLE) {
X xcpt.retval = HALFPI - atan (1.0 / x);
X } else {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: UNDERFLOW error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X }
X } else if (x > LAST_BOUND) {
X t1 = x * SQRT3 - 1.0;
X t2 = SQRT3 + x;
X xcpt.retval = SIXTHPI + atan (t1 / t2);
X } else if (x < X16_UNDERFLOWS) {
X xcpt.type = PLOSS;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: PLOSS error\n", funcname);
X errno = EDOM;
X xcpt.retval = x;
X }
X } else {
X xt2 = x * x;
X order = sizeof (atan_coeffs) / sizeof(double);
X order -= 1;
X xcpt.retval = x * poly (order, atan_coeffs, xt2);
X }
X DBUG_3 ("atanout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/atan.c
if test 5187 -ne `wc -c <funcs/src/atan.c`; then
echo shar: \"funcs/src/atan.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cos.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cos.c\"
else
echo shar: Extracting \"funcs/src/cos.c\" \(5011 characters\)
sed "s/^X//" >funcs/src/cos.c <<'END_OF_funcs/src/cos.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cos double precision cosine
X *
X * KEY WORDS
X *
X * cos
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double cos (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes cos(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call cos recursively
X * using relation cos(x) = -cos(x - PI).
X *
X * (3) If x < -PI/2 then call cos recursively
X * using relation cos(x) = -cos(x + PI).
X *
X * (4) If x > PI/4 then call sin using
X * relation cos(x) = sin(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation cos(x) = sin(PI/2 + x).
X *
X * (6) If x would cause underflow in approx
X * evaluation arithmetic then return
X * sqrt(1.0 - x**2).
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 119 can be used:
X *
X * cos(x) = ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.12905394659037374438571854e+7
X * P1 = -0.3745670391572320471032359e+6
X * P2 = 0.134323009865390842853673e+5
X * P3 = -0.112314508233409330923e+3
X * Q0 = 0.12905394659037373590295914e+7
X * Q1 = 0.234677731072458350524124e+5
X * Q2 = 0.2096951819672630628621e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3843 pg 244)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X *
X */
X
X# include <stdio.h>
X# include <pmluser.h>
X# include "pml.h"
X
Xstatic double cos_pcoeffs[] = {
X 0.12905394659037374438e7,
X -0.37456703915723204710e6,
X 0.13432300986539084285e5,
X -0.11231450823340933092e3
X};
X
Xstatic double cos_qcoeffs[] = {
X 0.12905394659037373590e7,
X 0.23467773107245835052e5,
X 0.20969518196726306286e3,
X 1.0
X};
X
Xstatic char funcname[] = "cos";
X
X
Xdouble cos (x)
Xdouble x;
X{
X auto double y;
X auto double yt2;
X auto double rtnval;
X extern double mod ();
X extern double sin ();
X extern double sqrt ();
X extern double poly ();
X struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("cosin", "arg %le", x);
X if (x < -PI || x > PI) {
X x = mod (x, TWOPI);
X if (x > PI) {
X x = x - TWOPI;
X } else if (x < -PI) {
X x = x + TWOPI;
X }
X }
X if (x > HALFPI) {
X xcpt.retval = -(cos (x - PI));
X } else if (x < -HALFPI) {
X xcpt.retval = -(cos (x + PI));
X } else if (x > FOURTHPI) {
X xcpt.retval = sin (HALFPI - x);
X } else if (x < -FOURTHPI) {
X xcpt.retval = sin (HALFPI + x);
X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
X xcpt.retval = sqrt (1.0 - (x * x));
X } else {
X y = x / FOURTHPI;
X yt2 = y * y;
X xcpt.retval = poly (3, cos_pcoeffs, yt2) / poly (3, cos_qcoeffs, yt2);
X }
X DBUG_3 ("cosout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/cos.c
if test 5011 -ne `wc -c <funcs/src/cos.c`; then
echo shar: \"funcs/src/cos.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/log.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/log.c\"
else
echo shar: Extracting \"funcs/src/log.c\" \(4649 characters\)
sed "s/^X//" >funcs/src/log.c <<'END_OF_funcs/src/log.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * log double precision natural log
X *
X * KEY WORDS
X *
X * log
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision natural log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 105-111.
X *
X * RESTRICTIONS
X *
X * The absolute error in the approximating polynomial is
X * 10**(-19.38). Note that contrary to DEC's assertion
X * in the F4P user's guide, the error is absolute and not
X * relative.
X *
X * This error bound assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes log(X) from:
X *
X * (1) If argument is zero then flag an error
X * and return minus infinity (or rather its
X * machine representation).
X *
X * (2) If argument is negative then flag an
X * error and return minus infinity.
X *
X * (3) Given that x = m * 2**k then extract
X * mantissa m and exponent k.
X *
X * (4) Transform m which is in range [0.5,1.0]
X * to range [1/sqrt(2), sqrt(2)] by:
X *
X * s = m * sqrt(2)
X *
X * (5) Compute z = (s - 1) / (s + 1)
X *
X * (6) Now use the approximation from HART
X * page 111 to find log(s):
X *
X * log(s) = z * ( P(z**2) / Q(z**2) )
X *
X * Where:
X *
X * P(z**2) = SUM [ Pj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * Q(z**2) = SUM [ Qj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = -0.240139179559210509868484e2
X * P1 = 0.30957292821537650062264e2
X * P2 = -0.96376909336868659324e1
X * P3 = 0.4210873712179797145
X * Q0 = -0.120069589779605254717525e2
X * Q1 = 0.19480966070088973051623e2
X * Q2 = -0.89111090279378312337e1
X * Q3 = 1.0000
X *
X * (coefficients from HART table #2705 pg 229)
X *
X * (7) Finally, compute log(x) from:
X *
X * log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double log_pcoeffs[] = {
X -0.24013917955921050986e2,
X 0.30957292821537650062e2,
X -0.96376909336868659324e1,
X 0.4210873712179797145
X};
X
Xstatic double log_qcoeffs[] = {
X -0.12006958977960525471e2,
X 0.19480966070088973051e2,
X -0.89111090279378312337e1,
X 1.0000
X};
X
Xstatic char funcname[] = "log";
X
X
Xdouble log (x)
Xdouble x;
X{
X auto int k;
X auto double s;
X auto double z;
X auto double zt2;
X auto double pqofz;
X auto struct exception xcpt;
X extern double frexp ();
X extern double poly ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("login", "arg %le", x);
X if (x == 0.0) {
X xcpt.type = SING;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: SINGULARITY error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else {
X s = SQRT2 * frexp (x, &k);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "s = %le", s);
X z = (s - 1.0) / (s + 1.0);
X DBUG_3 ("log", "z = %le", z);
X zt2 = z * z;
X DBUG_3 ("log", "zt2 = %le", zt2);
X pqofz = z * (poly (3, log_pcoeffs, zt2) / poly (3, log_qcoeffs, zt2));
X DBUG_3 ("pqofz", "pqofz = %le", pqofz);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "LN2 = %le", LN2);
X x = k * LN2;
X DBUG_3 ("log", "x = %le", x);
X x -= LNSQRT2;
X DBUG_3 ("log", "x = %le", x);
X x += pqofz;
X DBUG_3 ("log", "x = %le", x);
X xcpt.retval = x;
X }
X DBUG_3 ("logout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/log.c
if test 4649 -ne `wc -c <funcs/src/log.c`; then
echo shar: \"funcs/src/log.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/pml.h -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/pml.h\"
else
echo shar: Extracting \"funcs/src/pml.h\" \(3819 characters\)
sed "s/^X//" >funcs/src/pml.h <<'END_OF_funcs/src/pml.h'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * This file gets included with all of the floating point math
X * library routines when they are compiled. Note that
X * this is the proper place to put machine dependencies
X * whenever possible.
X *
X * It should be pointed out that for simplicity's sake, the
X * environment parameters are defined as floating point constants,
X * rather than octal or hexadecimal initializations of allocated
X * storage areas. This means that the range of allowed numbers
X * may not exactly match the hardware's capabilities. For example,
X * if the maximum positive double precision floating point number
X * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
X * defined to be 1.11E100 then the numbers between 1.11E100 and
X * 1.11...E100 are considered to be undefined. For most
X * applications, this will cause no problems.
X *
X * An alternate method is to allocate a global static "double" variable,
X * say "maxdouble", and use a union declaration and initialization
X * to initialize it with the proper bits for the EXACT maximum value.
X * This was not done because the only compilers available to the
X * author did not fully support union initialization features.
X *
X */
X
X#ifndef NO_DBUG
X# include <dbug.h>
X#else
X# define DBUG_ENTER(a1)
X# define DBUG_RETURN(a1) return(a1)
X# define DBUG_VOID_RETURN return
X# define DBUG_EXECUTE(keyword,a1)
X# define DBUG_2(keyword,format)
X# define DBUG_3(keyword,format,a1)
X# define DBUG_4(keyword,format,a1,a2)
X# define DBUG_5(keyword,format,a1,a2,a3)
X# define DBUG_PUSH(a1)
X# define DBUG_POP()
X# define DBUG_PROCESS(a1)
X# define DBUG_FILE (stderr)
X#endif
X
X#include <values.h>
X#include <math.h>
X#include <errno.h>
X
X/*
X * MAXDOUBLE => Maximum double precision number
X * MINDOUBLE => Minimum double precision number
X * DMAXEXP => Maximum exponent of a double
X * DMINEXP => Minimum exponent of a double
X */
X
X#define MINDOUBLE (1.0/MAXDOUBLE)
X#define LOG2_MAXDOUBLE (DMAXEXP + 1)
X#define LOG2_MINDOUBLE (DMINEXP - 1)
X#define LOGE_MAXDOUBLE (LOG2_MAXDOUBLE / LOG2E)
X#define LOGE_MINDOUBLE (LOG2_MINDOUBLE / LOG2E)
X
X/*
X * The following are hacks which should be fixed when I understand all
X * the issues a little better. |tanh(TANH_MAXARG)| = 1.0
X */
X#define TANH_MAXARG 16
X#define SQRT_MAXDOUBLE 1.304380e19
X
X#define PI 3.14159265358979323846
X#define TWOPI (2.0 * PI)
X#define HALFPI (PI / 2.0)
X#define FOURTHPI (PI / 4.0)
X#define SIXTHPI (PI / 6.0)
X#define LOG2E 1.4426950408889634074 /* Log to base 2 of e */
X#define LOG10E 0.4342944819032518276
X#define SQRT2 1.4142135623730950488
X#define SQRT3 1.7320508075688772935
X#define LN2 0.6931471805599453094
X#define LNSQRT2 0.3465735902799726547
X
X
X/*
X * MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X */
X
X#ifdef IEEE
X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define VOID void
END_OF_funcs/src/pml.h
if test 3819 -ne `wc -c <funcs/src/pml.h`; then
echo shar: \"funcs/src/pml.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/sin.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/sin.c\"
else
echo shar: Extracting \"funcs/src/sin.c\" \(4983 characters\)
sed "s/^X//" >funcs/src/sin.c <<'END_OF_funcs/src/sin.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sin double precision sine
X *
X * KEY WORDS
X *
X * sin
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision sine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sin (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes sin(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call sin recursively
X * using relation sin(x) = -sin(x - PI).
X *
X * (3) If x < -PI/2 then call sin recursively
X * using relation sin(x) = -sin(x + PI).
X *
X * (4) If x > PI/4 then call cos using
X * relation sin(x) = cos(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation sin(x) = -cos(PI/2 + x).
X *
X * (6) If x is small enough that polynomial
X * evaluation would cause underflow
X * then return x, since sin(x)
X * approaches x as x approaches zero.
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 118 can be used:
X *
X * sin(x) = y * ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.206643433369958582409167054e+7
X * P1 = -0.18160398797407332550219213e+6
X * P2 = 0.359993069496361883172836e+4
X * P3 = -0.2010748329458861571949e+2
X * Q0 = 0.263106591026476989637710307e+7
X * Q1 = 0.3927024277464900030883986e+5
X * Q2 = 0.27811919481083844087953e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3063 pg 234)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double sin_pcoeffs[] = {
X 0.20664343336995858240e7,
X -0.18160398797407332550e6,
X 0.35999306949636188317e4,
X -0.20107483294588615719e2
X};
X
Xstatic double sin_qcoeffs[] = {
X 0.26310659102647698963e7,
X 0.39270242774649000308e5,
X 0.27811919481083844087e3,
X 1.0
X};
X
Xstatic char funcname[] = "sin";
X
X
Xdouble sin (x)
Xdouble x;
X{
X double y;
X double yt2;
X double rtnval;
X extern double mod ();
X extern double cos ();
X extern double poly ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("sinin", "arg %le", x);
X if (x < -PI || x > PI) {
X x = mod (x, TWOPI);
X if (x > PI) {
X x = x - TWOPI;
X } else if (x < -PI) {
X x = x + TWOPI;
X }
X }
X if (x > HALFPI) {
X xcpt.retval = -(sin (x - PI));
X } else if (x < -HALFPI) {
X xcpt.retval = -(sin (x + PI));
X } else if (x > FOURTHPI) {
X xcpt.retval = cos (HALFPI - x);
X } else if (x < -FOURTHPI) {
X xcpt.retval = -(cos (HALFPI + x));
X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
X xcpt.retval = x;
X } else {
X y = x / FOURTHPI;
X yt2 = y * y;
X xcpt.retval = y * (poly (3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2));
X }
X DBUG_3 ("sinout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/sin.c
if test 4983 -ne `wc -c <funcs/src/sin.c`; then
echo shar: \"funcs/src/sin.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/sqrt.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/sqrt.c\"
else
echo shar: Extracting \"funcs/src/sqrt.c\" \(4591 characters\)
sed "s/^X//" >funcs/src/sqrt.c <<'END_OF_funcs/src/sqrt.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sqrt double precision square root
X *
X * KEY WORDS
X *
X * sqrt
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision square root of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sqrt (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 89-96.
X *
X * RESTRICTIONS
X *
X * The relative error is 10**(-30.1) after three applications
X * of Heron's iteration for the square root.
X *
X * However, this assumes exact arithmetic in the iterations
X * and initial approximation. Additional errors may occur
X * due to truncation, rounding, or machine precision limits.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes square root by:
X *
X * (1) Range reduction of argument to [0.5,1.0]
X * by application of identity:
X *
X * sqrt(x) = 2**(k/2) * sqrt(x * 2**(-k))
X *
X * k is the exponent when x is written as
X * a mantissa times a power of 2 (m * 2**k).
X * It is assumed that the mantissa is
X * already normalized (0.5 =< m < 1.0).
X *
X * (2) An approximation to sqrt(m) is obtained
X * from:
X *
X * u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
X *
X * P0 = 0.594604482
X * P1 = 2.54164041
X * Q0 = 2.13725758
X * Q1 = 1.0
X *
X * (coefficients from HART table #350 pg 193)
X *
X * (3) Three applications of Heron's iteration are
X * performed using:
X *
X * y[n+1] = 0.5 * (y[n] + (m/y[n]))
X *
X * where y[0] = u = approx sqrt(m)
X *
X * (4) If the value of k was odd then y is either
X * multiplied by the square root of two or
X * divided by the square root of two for k positive
X * or negative respectively. This rescales y
X * by multiplying by 2**frac(k/2).
X *
X * (5) Finally, y is rescaled by int(k/2) which
X * is equivalent to multiplication by 2**int(k/2).
X *
X * The result of steps 4 and 5 is that the value
X * of y between 0.5 and 1.0 has been rescaled by
X * 2**(k/2) which removes the original rescaling
X * done prior to finding the mantissa square root.
X *
X * NOTES
X *
X * The Convergent Technologies compiler optimizes division
X * by powers of two to "arithmetic shift right" instructions.
X * This is ok for positive integers but gives different
X * results than other C compilers for negative integers.
X * For example, (-1)/2 is -1 on the CT box, 0 on every other
X * machine I tried.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#define P0 0.594604482 /* Approximation coeff */
X#define P1 2.54164041 /* Approximation coeff */
X#define Q0 2.13725758 /* Approximation coeff */
X#define Q1 1.0 /* Approximation coeff */
X
X#define ITERATIONS 3 /* Number of iterations */
X
Xstatic char funcname[] = "sqrt";
X
X
Xdouble sqrt (x)
Xdouble x;
X{
X auto int k;
X register int bugfix;
X register int kmod2;
X register int count;
X auto int exponent;
X auto double m;
X auto double u;
X auto double y;
X auto double rtnval;
X auto struct exception xcpt;
X extern double frexp ();
X extern double ldexp ();
X
X DBUG_ENTER ("sqrt");
X DBUG_3 ("sqrtin", "arg %le", x);
X if (x == 0.0) {
X rtnval = 0.0;
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X } else {
X m = frexp (x, &k);
X u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
X for (count = 0, y = u; count < ITERATIONS; count++) {
X y = 0.5 * (y + (m / y));
X }
X if ((kmod2 = (k % 2)) < 0) {
X y /= SQRT2;
X } else if (kmod2 > 0) {
X y *= SQRT2;
X }
X bugfix = 2;
X xcpt.retval = ldexp (y, k/bugfix);
X }
X DBUG_3 ("sqrtout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/sqrt.c
if test 4591 -ne `wc -c <funcs/src/sqrt.c`; then
echo shar: \"funcs/src/sqrt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/scale.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/scale.c\"
else
echo shar: Extracting \"funcs/unused/scale.c\" \(4092 characters\)
sed "s/^X//" >funcs/unused/scale.c <<'END_OF_funcs/unused/scale.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * scale scale a double precision number by power of 2
X *
X * KEY WORDS
X *
X * scale
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * double scale (value, scale)
X * double value;
X * integer scale;
X *
X * DESCRIPTION
X *
X * Adds a specified integer to a double precision number's
X * exponent, effectively multiplying by a power of 2 for positive
X * scale values and divided by a power of 2 for negative
X * scale values.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl
X * Tempe, Az. 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and the 68000 the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X * Scale extracts the exponent, shifts it into the lower
X * order bits, adds the scale value to it, checks for
X * exponent overflow or underflow, shifts it back into the
X * high bits, and inserts it back into the value.
X *
X * NOTE: On the DECSYSTEM-20, the exponent for negative
X * numbers is complemented. This is not true for the PDP-11.
X *
X * NOTE: On the 68000, the mantissa is the same for both
X * positive and negative numbers, only the sign is different.
X * On the PDP-11, the mantissa is two's complement. Both
X * use hidden bit normalization.
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#ifdef PDP10
X#define EXP_MASK 0377000000000 /* Mask for exponent */
X#define MANT_MASK 0400777777777 /* Mask for mantissa */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#define EXP_SHIFTS 27 /* Shifts for exp in LSBs */
X#endif
X
X#ifdef pdp11
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#endif
X
X#ifdef mc68000
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#endif
X
Xstatic union {
X double dval;
X long lval[2];
X} share;
X
X
X
Xdouble scale(value,scale)
Xdouble value;
Xregister int scale;
X{
X register long temp1, temp2, *lpntr;
X
X lpntr = &share.lval[0];
X share.dval = value;
X temp1 = *lpntr;
X temp2 = (temp1 & EXP_MASK) >> EXP_SHIFTS;
X#ifdef PDP10
X if (value < 0.0) {
X temp2 = ~temp2 & LEXP_MASK;
X }
X#endif
X temp2 += scale;
X if (temp2 > MAX_EXPONENT+128) {
X pmlerr(SCALE_OVERFLOW);
X } else if (temp2 < 0) {
X pmlerr(SCALE_UNDERFLOW);
X } else {
X#ifdef PDP10
X if (value < 0.0) {
X temp2 = ~temp2 & LEXP_MASK;
X }
X#endif
X temp1 &= MANT_MASK;
X temp2 = temp2 << EXP_SHIFTS;
X *lpntr = temp1 | temp2;
X }
X return (share.dval);
X}
END_OF_funcs/unused/scale.c
if test 4092 -ne `wc -c <funcs/unused/scale.c`; then
echo shar: \"funcs/unused/scale.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/xexp.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/xexp.c\"
else
echo shar: Extracting \"funcs/unused/xexp.c\" \(3294 characters\)
sed "s/^X//" >funcs/unused/xexp.c <<'END_OF_funcs/unused/xexp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * xexp extract double precision number's exponent
X *
X * KEY WORDS
X *
X * xexp
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * int xexp(value)
X * double value;
X *
X * DESCRIPTION
X *
X * Extracts exponent from a double precision number and
X * returns it as a normal length integer.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl.
X * Tempe, Az 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and mc68000 (non IEEE format) the double
X * precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * For the mc68000 using IEEE format the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEEEEMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
X#ifdef PDP10
X#define EXP_MASK 0377000000000 /* Mask for exponent */
X#define EXP_SHIFTS 27 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X
X#ifdef mc68000
X#ifdef IEEE
X#define EXP_MASK 0x7FF00000 /* Mask for exponent */
X#define EXP_SHIFTS 20 /* Shifts to get into LSB's */
X#define EXP_BIAS 1023 /* Exponent bias */
X#else
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X#endif
X
X#ifdef pdp11
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X
Xunion dtol {
X double dval;
X int ival[2];
X};
X
X
Xint xexp (value)
Xunion dtol value;
X{
X register int *ipntr;
X
X if (value.dval == 0.0) {
X return (0);
X } else {
X ipntr = &value.ival[0];
X#ifdef PDP10
X if (value.dval < 0.0) { /* Exponent is complemented */
X *ipntr = ~*ipntr;
X }
X#endif
X *ipntr &= EXP_MASK;
X *ipntr >>= EXP_SHIFTS;
X *ipntr -= EXP_BIAS;
X return (*ipntr);
X }
X}
END_OF_funcs/unused/xexp.c
if test 3294 -ne `wc -c <funcs/unused/xexp.c`; then
echo shar: \"funcs/unused/xexp.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/xmant.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/xmant.c\"
else
echo shar: Extracting \"funcs/unused/xmant.c\" \(3950 characters\)
sed "s/^X//" >funcs/unused/xmant.c <<'END_OF_funcs/unused/xmant.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * xmant extract double precision number's mantissa
X *
X * KEY WORDS
X *
X * xmant
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * double xmant (value)
X * double value;
X *
X * DESCRIPTION
X *
X * Extracts a double precision number's mantissa and returns it
X * as a double precision number with a normalized mantissa and
X * a zero exponent.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl
X * Tempe, Az 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and the mc68000 the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * For the mc68000 using IEEE format the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEEEEMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X * NOTE: Beware of 0.0; on some machines which use excess 128
X * notation for the exponent, if the mantissa is zero the exponent
X * is also.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#ifdef PDP10
X#define MANT_MASK 0400777777777 /* Mantissa extraction mask */
X#define ZPOS_MASK 0200000000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0177000000000 /* Negative # mask for exp = 0 */
X#endif
X
X#ifdef pdp11
X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */
X#endif
X
X#ifdef mc68000
X#ifdef IEEE
X#define MANT_MASK 0x800FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x3FF00000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x3FF00000 /* Negative # mask for exp = 0 */
X#else
X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */
X#endif
X#endif
X
Xunion dtol {
X double dval;
X int ival[2];
X};
X
X
Xdouble xmant (value)
Xunion dtol value;
X{
X register int *ipntr;
X
X ipntr = &value.ival[0];
X *ipntr &= MANT_MASK;
X *ipntr |= ZPOS_MASK;
X return (value.dval);
X}
X
X
X/****************** OLD WAY ***********
Xdouble xmant (value)
Xunion dtol value;
X{
X register long *ipntr;
X register int neg;
X
X if (value == 0.0) {
X return (value);
X } else {
X if (value < 0.0) {
X neg = TRUE;
X } else {
X neg = FALSE;
X }
X ipntr = &value.ival[0];
X *ipntr &= MANT_MASK;
X if (neg) {
X *ipntr |= ZNEG_MASK;
X } else {
X *ipntr |= ZPOS_MASK;
X }
X return (value.dval);
X }
X}
X**********************/
END_OF_funcs/unused/xmant.c
if test 3950 -ne `wc -c <funcs/unused/xmant.c`; then
echo shar: \"funcs/unused/xmant.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/c2c.dat -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/c2c.dat\"
else
echo shar: Extracting \"tests/c2c.dat\" \(5354 characters\)
sed "s/^X//" >tests/c2c.dat <<'END_OF_tests/c2c.dat'
Xcsin 1.0e00 2.0e00 3.165778547525405883e00 1.959601044654846191e00
Xcsin 0.0e00 2.0e00 0.0e00 3.626860409975051879e00
Xcsin 1.0e00 0.0e00 8.414709866046905517e-01 0.0e00
Xcsin -2.0e00 2.0e00 -3.420954853296279907e00 -1.509306490421295166e00
Xcsin -1.0e00 -2.0e00 -3.165778547525405883e00 -1.959601044654846191e00
Xcsin 1.0e00 -2.0e00 3.165778547525405883e00 -1.959601044654846191e00
Xccos 1.0e00 2.0e00 2.032723009586334228e00 -3.051897794008255004e00
Xccos 0.0e00 2.0e00 3.762195706367492675e00 0.0e00
Xccos 1.0e00 0.0e00 5.403023064136505127e-01 0.0e00
Xccos -2.0e00 2.0e00 -1.565625846385955810e00 3.297894805669784545e00
Xccos -1.0e00 -2.0e00 2.032723009586334228e00 -3.051897794008255004e00
Xccos 1.0e00 -2.0e00 2.032723009586334228e00 3.051897794008255004e00
Xcexp 1.0e00 2.0e00 -1.131204381585121154e00 2.471726655960083007e00
Xcexp 0.0e00 2.0e00 -4.161468371748924255e-01 9.092974215745925903e-01
Xcexp 1.0e00 0.0e00 2.718281835317611694e00 0.0e00
Xcexp -2.0e00 2.0e00 -5.631934991106390953e-02 1.230600243434309959e-01
Xcexp -1.0e00 -2.0e00 -1.530918665230274200e-01 -3.345118276774883270e-01
Xcexp 1.0e00 -2.0e00 -1.131204381585121154e00 -2.471726655960083007e00
Xctanh 1.0e00 2.0e00 1.166736245155334472e00 -2.434581927955150604e-01
Xctanh 0.0e00 2.0e00 0.0e00 -2.185039848089218139e00
Xctanh 1.0e00 0.0e00 7.615941539406776428e-01 0.0e00
Xctanh -2.0e00 2.0e00 -1.023835599422454834e00 -2.839295566082000732e-02
Xctanh -1.0e00 -2.0e00 -1.166736245155334472e00 2.434581927955150604e-01
Xctanh 1.0e00 -2.0e00 1.166736245155334472e00 2.434581927955150604e-01
Xctan 1.0e00 2.0e00 3.381283208727836608e-02 1.014793619513511657e00
Xctan 0.0e00 2.0e00 0.0e00 9.640275761485099792e-01
Xctan 1.0e00 0.0e00 1.557407721877098083e00 0.0e00
Xctan -2.0e00 2.0e00 2.839295566082000732e-02 1.023835599422454834e00
Xctan -1.0e00 -2.0e00 -3.381283208727836608e-02 -1.014793619513511657e00
Xctan 1.0e00 -2.0e00 3.381283208727836608e-02 -1.014793619513511657e00
Xcsqrt 0.0 0.0 0.000000000000000000e00 0.000000000000000000e00
Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01
Xcsqrt 0.0 2.0 1.000000000000000000e00 1.000000000000000000e00
Xcsqrt 1.0 0.0 1.000000000000000000e00 0.000000000000000000e00
Xcsqrt -2.0 2.0 6.435942575335502624e-01 1.553773969411849975e00
Xcsqrt -1.0 2.0 7.861513718962669372e-01 1.272019654512405395e00
Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01
Xcsinh 1.0e00 2.0e00 -4.890562593936920166e-01 1.403119236230850219e00
Xcsinh 0.0e00 2.0e00 0.0e00 9.092974215745925903e-01
Xcsinh 1.0e00 0.0e00 1.175201192498207092e00 0.0e00
Xcsinh -2.0e00 2.0e00 1.509306475520133972e00 3.420954853296279907e00
Xcsinh -1.0e00 -2.0e00 4.890562593936920166e-01 -1.403119236230850219e00
Xcsinh 1.0e00 -2.0e00 -4.890562593936920166e-01 -1.403119236230850219e00
Xcrcp 1.0e00 2.0e00 1.999999992549419403e-01 -3.999999985098838806e-01
Xcrcp 0.0e00 2.0e00 0.0e00 -5.0e-01
Xcrcp 1.0e00 0.0e00 1.0e00 0.0e00
Xcrcp -2.0e00 2.0e00 -2.5e-01 -2.5e-01
Xcrcp -1.0e00 -2.0e00 -1.999999992549419403e-01 3.999999985098838806e-01
Xcrcp 1.0e00 -2.0e00 1.999999992549419403e-01 3.999999985098838806e-01
Xclog 1.0e00 2.0e00 8.047189563512802124e-01 1.107148721814155578e00
Xclog 0.0e00 2.0e00 6.931471824645996093e-01 1.570796325802803039e00
Xclog 1.0e00 0.0e00 0.0e00 0.0e00
Xclog -2.0e00 2.0e00 1.039720773696899414e00 2.356194496154785156e00
Xclog -1.0e00 -2.0e00 8.047189563512802124e-01 -2.034443944692611694e00
Xclog 1.0e00 -2.0e00 8.047189563512802124e-01 -1.107148721814155578e00
Xccosh 1.0e00 2.0e0 -6.421481221914291381e-01 1.068607419729232788e00
Xccosh 0.0e00 2.0e00 -4.161468371748924255e-01 0.0e00
Xccosh 1.0e00 0.0e00 1.543080642819404602e00 0.0e00
Xccosh -2.0e00 2.0e00 -1.565625831484794616e00 -3.297894805669784545e00
Xccosh -1.0e00 -2.0e00 -6.421481221914291381e-01 1.068607419729232788e00
Xccosh 1.0e00 -2.0e00 -6.421481221914291381e-01 -1.068607419729232788e00
Xcatan 1.0e00 2.0e00 1.338972523808479309e00 4.023594781756401062e-01
Xcatan 0.0e00 2.0e00 1.570796325802803039e00 5.493061468005180358e-01 /* -real part */
Xcatan 1.0e00 0.0e00 7.853981629014015197e-01 0.0e00
Xcatan -2.0e00 2.0e00 -1.311223268508911132e00 2.388778626918792724e-01
Xcatan -1.0e00 -2.0e00 -1.338972523808479309e00 -4.023594819009304046e-01
Xcatan 1.0e00 -2.0e00 1.338972523808479309e00 -4.023594819009304046e-01
Xcasin 1.0 2.0 4.270785786211490631e-01 1.528570890426635742e00
Xcasin 0.0 2.0 0.000000000000000000e00 1.443635478615760803e00
Xcasin 1.0 0.0 1.570796325802803039e00 0.000000000000000000e00
Xcasin -2.0 2.0 -7.542491331696510314e-01 1.734324455261230468e00
Xcasin -1.0 -2.0 -4.270785823464393615e-01 -1.528570920228958129e00
Xcasin 1.0 -2.0 4.270785823464393615e-01 -1.528570920228958129e00
Xcacos 1.0 2.0 1.143717750906944274e00 -1.528570920228958129e00
Xcacos 0.0 2.0 1.570796325802803039e00 -1.443635478615760803e00
Xcacos 1.0 0.0 0.000000000000000000e00 0.000000000000000000e00
Xcacos -2.0 2.0 2.325045466423034668e00 -1.734324529767036438e00
Xcacos -1.0 -2.0 1.997874900698661804e00 1.528570890426635742e00
Xcacos 1.0 -2.0 1.143717750906944274e00 1.528570890426635742e00
END_OF_tests/c2c.dat
if test 5354 -ne `wc -c <tests/c2c.dat`; then
echo shar: \"tests/c2c.dat\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/cc2c.dat -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/cc2c.dat\"
else
echo shar: Extracting \"tests/cc2c.dat\" \(3825 characters\)
sed "s/^X//" >tests/cc2c.dat <<'END_OF_tests/cc2c.dat'
Xcadd 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.683000028133392334e00 6.440000057220458984e00
Xcadd -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 2.003005027770996093e00 -1.887652993202209472e00
Xcadd 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 3.719695717096328735e00 4.449983000755310058e00
Xcadd 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.954332172870635986e00 0.000000000000000000e00
Xcadd -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -4.234000027179718017e00 -6.334455549716949462e00
Xcadd 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.300034515380859375e03
Xcdiv 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.609979800879955291e-01 1.263787318021059036e-01
Xcdiv -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.863302707672119140e-01 6.862613558769226074e-02
Xcdiv 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 5.729839205741882324e-02 -4.476501746103167533e-02
Xcdiv 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.663971476256847381e-01 0.000000000000000000e00
Xcdiv -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 5.215928852558135986e-01 8.269466646015644073e-02
Xcdiv 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 3.765932763671875000e04 -3.554461171875000000e05
Xcmult 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -5.596120119094848632e00 1.293470001220703125e01
Xcmult -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 4.705165147781372070e00 1.302998638153076171e01
Xcmult 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 8.510770574212074279e-01 2.011436134576797485e00
Xcmult 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 5.323814511299133300e00 0.000000000000000000e00
Xcmult -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -5.635822236537933349e00 1.193936669826507568e01
Xcmult 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 -4.487600040435791015e01 4.235604019165039062e02
Xcsubt 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -2.437000006437301635e00 -1.760000020265579223e00
Xcsubt -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.694344997406005859e00 6.356786966323852539e00
Xcsubt 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 -2.948183119297027587e00 -4.249976992607116699e00
Xcsubt 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 -1.802817046642303466e00 0.000000000000000000e00
Xcsubt -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 1.766000002622604370e00 1.665544450283050537e00
Xcsubt 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.299965484619140625e03
END_OF_tests/cc2c.dat
if test 3825 -ne `wc -c <tests/cc2c.dat`; then
echo shar: \"tests/cc2c.dat\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 3 \(of 6\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 4 (of 6)."
# Contents: doc/header.r funcs/src/exp.c funcs/src/pmlerr.c
# funcs/src/unused/OLDpml.h tests/d2d.dat tests/errors.c
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:45 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f doc/header.r -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"doc/header.r\"
else
echo shar: Extracting \"doc/header.r\" \(7828 characters\)
sed "s/^X//" >doc/header.r <<'END_OF_doc/header.r'
X.tp
X.(1C
XP M L U S E R S G U I D E
X.sp 2
XP o r t a b l e M a t h L i b r a r y
X.sp 2
XVersion 2.0
X.sp 2
X\*(td
X.sp 4
Xby
X.sp
XFred Fish
X.)1
X.bp
X.sh 1
XINTRODUCTION
X.(x
Xdocument topics
X.)x
X.sh 2
XDocument structure
X.pp
XThis document describes the PML library, a
Xmathematical function library for C programs.
XIt contains the following major sections:
X.(1
XIntroduction to the document (this section).
XA module description section.
XAn index.
X.)1
X.sh 2
XDesign Considerations
X
X.pp
XIn writing this library many tradeoffs had to be considered. The
Xprimary design goal was transportability. It was desired that
Xthe final library be easily transportable among a wide class
Xof machines. As of this release, the library has been used
Xwith only minor modifications on
Xa DECSYSTEM-20 (a 36 bit machine), a VAX-11/780 under compatibility
Xmode, and various PDP-11's.
X
X.pp
XThis portability was achieved by careful isolation of machine
Xdependencies and parameterization of the environment (see references).
XThe only assumption made is that the C compiler can generate proper code for
Xthe four basic operations (+,-,/,*).
X
X.pp
XEven though efficiency
Xwas considered to be of only secondary importance, the final routines
Xcompared favorably with an informal test "race" against the
XDECSYSTEM-20 FORTRAN, which has optimized assembly language
Xlibrary routines. The PML library routines seldom took
Xmore than twice as long as the FORTRAN library, and many were
Xclose enough to call a draw.
X
X.pp
XThere are currently only four highly machine dependent routines in
Xthe Portable Math Library. When transporting the library to a new machine,
Xthese should be the only ones in which recoding is necessary.
XThese routines, written in machine targeted C, are:
X
X.pp
Xdscale --- scale a double precision floating point number by a specified
Xpower of two. Equivalent to multiplication or division
Xby a power of two, depending upon sign of the scale value.
X
X.pp
Xdxexp --- extract exponent of double precision floating point number and
Xreturn it as an integer.
X
X.pp
Xdxmant --- extract mantissa of double precision floating point number and
Xreturn it as a double precision floating point number.
X
X.pp
Xdint --- discard fractional part of double precision number, returning
Xinteger part as a double precision floating point
Xnumber.
X
X.pp
XThe entire Portable Math Library is built upon six "primitives"
Xwhich compute their values from polynomial approximations. All
Xothers can be defined in terms of the primitives using various
Xidentities. The primitives are (1) datan, (2) dexp, (3) dln,
X(4) dsin, (5) dcos, and (6) dsqrt. Strictly speaking, only
Xone of dsin and dcos could be chosen as a primitive and the
Xother defined with an appropriate identity. In this implementation
Xhowever, dsin and dcos call themselves recursively, and each other,
Xto perform range reduction.
X
X.sh 2
XError Handling
X
X.pp
XNo assumptions are made about whether the four basic operations are done by
Xhardware or software. Any overflows or underflows in the basic operations
Xare assumed to be handled by the environment, if at all. Pathological
Xcases in the library routines are trapped internally and control
Xis passed to an error handler routine "pmlerr"
X(which may be replaced with
Xone of the user's choosing) for error recovery.
X
X.pp
XThe default error handler is conceptually similar to the one used by DEC for
Xthe FORTRAN compilers. It contains an internal table which
Xallows various actions to be taken for each error recognized.
XCurrently each error has a corresponding flag word with three
Xbits, each bit assigned as follows:
X
X.pp
XCONTINUE --- If the bit is set control is returned to the calling
Xroutine after completion of error processing. Otherwise
Xthe task exits with an error status.
X
X.pp
XLOG --- If the log bit is set then an error warning message is
Xwritten to the standard error output channel prior to exiting
Xor continuing. If reset, no message is given.
X
X.pp
XCOUNT --- If the count bit is set then the task's PML
Xerror count (internal to the error handler) is incremented.
XIf the total error count exceeds the maximum allowed
Xthen the task exits with
Xerror status. If the count bit is reset then the error
Xis ignored with respect to the error count and exit on limit.
X
X.pp
XThe default conditions in the error handler for all errors
Xis that the CONTINUE, LOG, and COUNT bits are all set. The
Xerror limit is set at 10. These values can be changed by
Xsuitably editing the header file "pml.h"
Xand <pmluser.h>.
X
X.pp
XThe error handler responses can also be changed dynamically
Xvia the following routines:
X
X.pp
XPMLSFS --- Sets the specified bits in the specified error's
Xflag word.
XFor example, "pmlsfs(NEG__DSQRT,CONTINUE | LOG)" sets
Xthe CONTINUE and LOG bits for the "double precision square root
Xof a negative number" error. The COUNT bit is not affected.
XThe manifest constant values are defined in <pmluser.h>.
X
X.pp
XPMLCFS --- Clears the specified bits in the specified error's
Xflag word.
XFor example, "pmlcfs(NEG__DSQRT,CONTINUE | LOG)" clears
Xthe CONTINUE and LOG bits for the "double precision square root
Xof a negative number" error. The COUNT bit is not affected.
XThe manifest constant values are defined in <pmluser.h>.
X
X.pp
XPMLLIM --- Changes the task's PML error limit to
Xthe specified value and returns the previous value.
X
X.pp
XPMLCNT --- Returns the current value of the PML error
Xcount and resets it back to zero.
X
X.sh 2
XFunction Names
X
X.pp
XIn general, FORTRAN naming conventions were followed since no
Xother more obvious convention was apparent. There is one
Xstrong exception however, and that is the natural
Xlogarithm functions use the generic name "ln" while the
Xlogarithm to the base 10 functions use the name "log".
XThis is consistent with the normal usage in virtually all
Xmodern mathematical and engineering texts. How FORTRAN came
Xto use "log" and "log10" respectively is beyond me.
XThis usage has bitten many a starting FORTRAN programmer.
X
X.sh 2
XInstallation
X
X.pp
XAs part of the installation kit, some simple minded testing programs
Xare provided. They are by no means exhaustive tests or ones that
Xcarefully check possible trouble areas. They are intended to
Xprovide a quick and dirty way of verifying that no gross errors
Xare inadvertently incorporated in the library as a result of
Ximprovements or bug fixes, and that the installation is successful.
XFuture releases may incorporate more extensive test routines
Xand suggestions are solicited.
X
X.sh 2
XBugs
X
X.pp
XOn the subject of bugs, all exterminators are encouraged
Xto notify the author of any problems or fixes so that
Xthey can be incorporated into the next release and renegade
Xversions of the library can be minimized.
XContact:
X.sp 2
X.nf
X Fred Fish
X 1346 W 10th Place
X Tempe, Arizona 85281
X.sp 1
X (602) 894-6881 (home)
X (602) 932-7000 (work @ GACA)
X.fi
X.sp 2
X.pp
XThose users with strong numerical analysis backgrounds or experience
Xmay be able to suggest better methods for some of the library routines.
XMany of the higher level routines are simple minded implementations
Xof identities, and may not be nearly as stable as some more obscure
Xmethods.
X
X.sh 2
XTransporting To Other Machines
X
X.pp
XTo transport the Portable Math Library to a processor other than
Xthe PDP11 or DECSYSTEM-20, do the following:
X.(1
XDefine the machine dependent parameters in "pml.h".
XImplement the four machine dependent modules listed previously.
XCompile the rest of the library modules and the testing routines.
XLink and run the testing routines to verify successful installation.
XRepeat as required.
X.)1
X
X.bp
X.sh 1
XRUNTIME MODULES
X
X.pp
XThe PML modules are documented on the following pages.
XAll the rest of this section was produced by the DEX (documentation
Xextractor) utility directly from the source files. Thus
Xthe information reflects the current state of the runtime
Xmodules.
END_OF_doc/header.r
if test 7828 -ne `wc -c <doc/header.r`; then
echo shar: \"doc/header.r\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/exp.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/exp.c\"
else
echo shar: Extracting \"funcs/src/exp.c\" \(5510 characters\)
sed "s/^X//" >funcs/src/exp.c <<'END_OF_funcs/src/exp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * exp double precision exponential
X *
X * KEY WORDS
X *
X * exp
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision exponential of double precision
X * floating point number.
X *
X * USAGE
X *
X * double exp (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 96-104.
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.4). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes exponential from:
X *
X * exp(x) = 2**y * 2**z * 2**w
X *
X * Where:
X *
X * y = int ( x * log2(e) )
X *
X * v = 16 * frac ( x * log2(e))
X *
X * z = (1/16) * int (v)
X *
X * w = (1/16) * frac (v)
X *
X * Note that:
X *
X * 0 =< v < 16
X *
X * z = {0, 1/16, 2/16, ...15/16}
X *
X * 0 =< w < 1/16
X *
X * Then:
X *
X * 2**z is looked up in a table of 2**0, 2**1/16, ...
X *
X * 2**w is computed from an approximation:
X *
X * 2**w = (A + B) / (A - B)
X *
X * A = C + (D * w * w)
X *
X * B = w * (E + (F * w * w))
X *
X * C = 20.8137711965230361973
X *
X * D = 1.0
X *
X * E = 7.2135034108448192083
X *
X * F = 0.057761135831801928
X *
X * Coefficients are from HART, table #1121, pg 206.
X *
X * Effective multiplication by 2**y is done by a
X * floating point scale with y as scale argument.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X# define C 20.8137711965230361973 /* Polynomial approx coeff. */
X# define D 1.0 /* Polynomial approx coeff. */
X# define E 7.2135034108448192083 /* Polynomial approx coeff. */
X# define F 0.057761135831801928 /* Polynomial approx coeff. */
X
X
X/************************************************************************
X * *
X * This table is fixed in size and reasonably hardware *
X * independent. The given constants were generated on a *
X * DECSYSTEM 20 using double precision FORTRAN. *
X * *
X ************************************************************************
X */
X
Xstatic double fpof2tbl[] = {
X 1.00000000000000000000, /* 2 ** 0/16 */
X 1.04427378242741384020, /* 2 ** 1/16 */
X 1.09050773266525765930, /* 2 ** 2/16 */
X 1.13878863475669165390, /* 2 ** 3/16 */
X 1.18920711500272106640, /* 2 ** 4/16 */
X 1.24185781207348404890, /* 2 ** 5/16 */
X 1.29683955465100966610, /* 2 ** 6/16 */
X 1.35425554693689272850, /* 2 ** 7/16 */
X 1.41421356237309504880, /* 2 ** 8/16 */
X 1.47682614593949931110, /* 2 ** 9/16 */
X 1.54221082540794082350, /* 2 ** 10/16 */
X 1.61049033194925430820, /* 2 ** 11/16 */
X 1.68179283050742908600, /* 2 ** 12/16 */
X 1.75625216037329948340, /* 2 ** 13/16 */
X 1.83400808640934246360, /* 2 ** 14/16 */
X 1.91520656139714729380 /* 2 ** 15/16 */
X};
X
Xstatic char funcname[] = "exp";
X
X
Xdouble exp (x)
Xdouble x;
X{
X register int y;
X register int index;
X auto double w;
X auto double v;
X auto double a;
X auto double b;
X auto double t;
X auto double temp;
X auto double wpof2;
X auto double zpof2;
X auto double rtnval;
X extern double dabs ();
X extern double ldexp ();
X extern double modf ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("expin", "arg %le", x);
X if (x > LOGE_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MAXDOUBLE;
X }
X } else if (x <= LOGE_MINDOUBLE) {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = 0.0;
X }
X } else {
X t = x * LOG2E;
X (void) modf (t, &temp);
X y = temp;
X v = 16.0 * modf (t, &temp);
X (void) modf (dabs (v), &temp);
X index = temp;
X if (x < 0.0) {
X zpof2 = 1.0 / fpof2tbl[index];
X } else {
X zpof2 = fpof2tbl[index];
X }
X w = modf (v, &temp) / 16.0;
X a = C + (D * w * w);
X b = w * (E + (F * w * w));
X wpof2 = (a + b) / (a - b);
X xcpt.retval = ldexp ((wpof2 * zpof2), y);
X }
X DBUG_3 ("expout", "result %le", rtnval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/exp.c
if test 5510 -ne `wc -c <funcs/src/exp.c`; then
echo shar: \"funcs/src/exp.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/pmlerr.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/pmlerr.c\"
else
echo shar: Extracting \"funcs/src/pmlerr.c\" \(7682 characters\)
sed "s/^X//" >funcs/src/pmlerr.c <<'END_OF_funcs/src/pmlerr.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/************************************************************************
X * *
X * PORTABLE MATH LIBRARY -- ERROR HANDLER *
X * *
X * This is a sample PML library error handler. *
X * It may be used as is, or another of the user's choice *
X * substituted. *
X * *
X * In any case, the global "pmlerr" must be defined somewhere *
X * in the user task, since many Portable Math Library routines *
X * reference it. The other routines in this file are not called *
X * by any library routines and may be omitted. *
X * *
X ************************************************************************
X */
X
X# include <stdio.h>
X# include <pmluser.h>
X# include "pml.h"
X
Xstatic struct pml_err {
X int flags; /* Flag word; bits defined in pml.h */
X char *msg; /* Error message */
X char *func; /* Function in which error occured */
X};
X
Xstatic struct pml_err pml_errs[] = {
X CONTINUE | COUNT | LOG, "overflow", "exp",
X CONTINUE | COUNT | LOG, "underflow", "exp",
X CONTINUE | COUNT | LOG, "exponent overflow", "scale",
X CONTINUE | COUNT | LOG, "negative argument", "sqrt",
X CONTINUE | COUNT | LOG, "zero argument", "log",
X CONTINUE | COUNT | LOG, "negative argument", "log",
X CONTINUE | COUNT | LOG, "argument magnitude greater than 1.0", "acos",
X CONTINUE | COUNT | LOG, "argument magnitude greater than 1.0", "asin",
X CONTINUE | COUNT | LOG, "overflow", "tan",
X CONTINUE | COUNT | LOG, "overflow", "cosh",
X CONTINUE | COUNT | LOG, "underflow", "cosh",
X CONTINUE | COUNT | LOG, "overflow", "sinh",
X CONTINUE | COUNT | LOG, "underflow", "sinh",
X CONTINUE | COUNT | LOG, "overflow", "asinh",
X CONTINUE | COUNT | LOG, "argument less than 1.0", "acosh",
X CONTINUE | COUNT | LOG, "overflow", "acosh",
X CONTINUE | COUNT | LOG, "argument magnitude not < 1.0", "atanh",
X CONTINUE | COUNT | LOG, "underflow", "atan",
X CONTINUE | COUNT | LOG, "complex division by zero", "cdiv",
X CONTINUE | COUNT | LOG, "complex reciprocal of zero", "crcp",
X CONTINUE | COUNT | LOG, "exponent underflow", "scale",
X CONTINUE | COUNT | LOG, "argument has no fractional part", "dint",
X};
X
Xstatic int err_count = 0; /* Counter for PML errors */
Xstatic int err_limit = MAX_ERRORS; /* PML error limit */
X
X/*
X * FUNCTION
X *
X * pmlcfs Clear specified PML error handler flags
X *
X * KEY WORDS
X *
X * pmlcfs
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Clear the specified PML error handler flags for the
X * specified error. Two or more flags may be cleared simultaneously
X * by "or-ing" them in the call, for example "LOG | CONTINUE".
X * The manifest constants for the flags and error codes are
X * defined in <pmluser.h>.
X *
X * USAGE
X *
X * pmlcfs(err_code,flags)
X * int err_code;
X * int flags;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xpmlcfs(err_code,flag_word)
Xregister int err_code;
Xregister int flag_word;
X{
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"pmlcfs: invalid error code %d\n",err_code);
X } else {
X pml_errs[err_code].flags &= ~flag_word;
X }
X}
X
X
X/*
X * FUNCTION
X *
X * pmlcnt get PML error count and reset it to zero
X *
X * KEY WORDS
X *
X * pmlcnt
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns the total number of PML errors seen
X * prior to the call, and resets the error count to zero.
X *
X * USAGE
X *
X * int pmlcnt()
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xint pmlcnt()
X{
X register int rtn_val;
X
X rtn_val = err_count;
X err_count = 0;
X return(rtn_val);
X}
X
X
X/*
X * FUNCTION
X *
X * pmlerr Portable Math Library error handler
X *
X * KEY WORDS
X *
X * pmlerr
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Provides a sample PML error handler. Does
X * not use any available hardware "traps" so is machine
X * independent. Generally only called internally by the
X * other PML routines.
X *
X * There are currently three flags which control the
X * response for specific errors:
X *
X * (1) LOG When set an error message is sent
X * to the user terminal.
X *
X * (2) COUNT When set the error is counted
X * against the PML error limit.
X *
X * (3) CONTINUE When set the task continues
X * providing the error count has not
X * exceeded the PML error limit.
X *
X * Each of these flags can be set or reset independently
X * by "pmlsfs" or "pmlcfs" respectively.
X *
X * USAGE
X *
X * pmlerr(err_code)
X * int err_code;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xpmlerr(err_code)
Xregister int err_code;
X{
X register struct pml_err *err;
X
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"pmlerr: invalid error code %d\n",err_code);
X } else {
X err = &pml_errs[err_code];
X if (err->flags & LOG) {
X fprintf(stderr,"pml: %s in function \"%s\"\n",err->msg,err->func);
X }
X if (err->flags & COUNT) {
X err_count++;
X }
X if ((err->flags & CONTINUE) && (err_count <= err_limit)) {
X return;
X } else {
X fprintf(stderr,"pml: error limit exceeded\n");
X fprintf(stderr,"pml: task aborted with %d error(s)\n",
X err_count);
X exit(-1);
X }
X }
X}
X
X
X/*
X * FUNCTION
X *
X * pmllim Set Portable Math Library error limit
X *
X * KEY WORDS
X *
X * pmllim
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Sets the PML error limit to the specified value
X * and returns it previous value.
X * Does not affect the current error count (which may be reset
X * to zero by a call to "pmlcnt"). Note that the default error
X * limit is set at compile time by the value in "pml.h".
X *
X * USAGE
X *
X * int pmllim(limit)
X * int limit;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xint pmllim(limit)
Xregister int limit;
X{
X register int rtn_val;
X
X rtn_val = err_limit;
X err_limit = limit;
X return(rtn_val);
X}
X
X
X/*
X * FUNCTION
X *
X * pmlsfs Set specified PML error handler flags
X *
X * KEY WORDS
X *
X * pmlsfs
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Set the specified PML error handler flags for the
X * specified error. Two or more flags may be set simultaneously
X * by "or-ing" them in the call, for example "LOG | CONTINUE".
X * The manifest constants for the flags and error codes are
X * defined in <pmluser.h>.
X *
X * USAGE
X *
X * pmlsfs(err_code,flags)
X * int err_code;
X * int flags;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
Xpmlsfs(err_code,flag_word)
Xregister int err_code;
Xregister int flag_word;
X{
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"? pmlsfs --- invalid error code %d.\n",err_code);
X } else {
X pml_errs[err_code].flags |= flag_word;
X }
X}
END_OF_funcs/src/pmlerr.c
if test 7682 -ne `wc -c <funcs/src/pmlerr.c`; then
echo shar: \"funcs/src/pmlerr.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/unused/OLDpml.h -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/unused/OLDpml.h\"
else
echo shar: Extracting \"funcs/src/unused/OLDpml.h\" \(6798 characters\)
sed "s/^X//" >funcs/src/unused/OLDpml.h <<'END_OF_funcs/src/unused/OLDpml.h'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * This file gets included with all of the floating point math
X * library routines when they are compiled. Note that
X * this is the proper place to put machine dependencies
X * whenever possible.
X *
X * It should be pointed out that for simplicity's sake, the
X * environment parameters are defined as floating point constants,
X * rather than octal or hexadecimal initializations of allocated
X * storage areas. This means that the range of allowed numbers
X * may not exactly match the hardware's capabilities. For example,
X * if the maximum positive double precision floating point number
X * is EXACTLY 1.11...E100 and the constant "MAX_POS_DBLF is
X * defined to be 1.11E100 then the numbers between 1.11E100 and
X * 1.11...E100 are considered to be undefined. For most
X * applications, this will cause no problems.
X *
X * An alternate method is to allocate a global static "double" variable,
X * say "max_pos_dblf", and use a union declaration and initialization
X * to initialize it with the proper bits for the EXACT maximum value.
X * This was not done because the only compilers available to the
X * author did not fully support union initialization features.
X *
X */
X
X
X/*
X * PDP10 (DECSYSTEM 20) HARDWARE DEPENDENCIES
X *
X * ******** W A R N I N G W A R N I N G ********
X * The current DEC20 implementation treats double precision
X * floats as single precision floats. When true double
X * precision is implemented, the flowing floats may have to
X * change.
X *
X * Since the PDP10 compiler has no static variables, and only
X * first 5 letters are significant, must use preprocessor
X * to avoid name clashes.
X */
X
X#ifndef NO_DBUG
X# include <dbug.h>
X#else
X# define DBUG_ENTER(a1)
X# define DBUG_RETURN(a1) return(a1)
X# define DBUG_VOID_RETURN return
X# define DBUG_EXECUTE(keyword,a1)
X# define DBUG_2(keyword,format)
X# define DBUG_3(keyword,format,a1)
X# define DBUG_4(keyword,format,a1,a2)
X# define DBUG_5(keyword,format,a1,a2,a3)
X# define DBUG_PUSH(a1)
X# define DBUG_POP()
X# define DBUG_PROCESS(a1)
X# define DBUG_FILE (stderr)
X#endif
X
X#ifdef PDP10
X#define MAX_POS_DBLF 1.7014118e38 /* Max positive double */
X#define MIN_POS_DBLF 1.4693680e-39 /* Min positive double */
X#define MAX_NEG_DBLF -1.7014118e38 /* Max negative double */
X#define MIN_NEG_DBLF -1.4693680e-39 /* Min negative double */
X#define MAX_EXPONENT 127 /* Max exponent allowed */
X#define RECIP_MIN 5.877471e-39 /* MAX_POS_DBLF >= 1/RECIP_MIN */
X#define RECIP_MAX 1.7014118e38 /* MIN_POS_DBLF <= 1/RECIP_MAX */
X#define LN_MAXPOSDBL 88.0 /* LN(MAX_POS_DBLF) */
X#define LN_MINPOSDBL -89.4 /* LN(MIN_POS_DBLF) */
X#define TANH_MAXARG 16 /* |TANH(maxarg)| = 1.0 */
X#define SQRT_MPDF 1.304380e19 /* SQRT(MAX_POS_DBLF) */
X#define X6_UNDERFLOWS 3.37174e-7 /* X**6 almost underflows */
X#define X16_UNDERFLOWS 3.74063e-3 /* X**16 almost underflows */
X
X#define atan_coeffs qzzz1 /* rename to avoid name clash */
X#define cos_pcoeffs qzzz2 /* rename to avoid name clash */
X#define cos_qcoeffs qzzz3 /* rename to avoid name clash */
X#define ln_pcoeffs qzzz4 /* rename to avoid name clash */
X#define ln_qcoeffs qzzz5 /* rename to avoid name clash */
X#define sin_pcoeffs qzzz6 /* rename to avoid name clash */
X#define sin_qcoeffs qzzz7 /* rename to avoid name clash */
X#endif
X
X
X/*
X * PDP11 HARDWARE DEPENDENCIES
X *
X */
X
X#ifdef pdp11
X#define MAX_POS_DBLF 1.7014118e38 /* Max positive double */
X#define MIN_POS_DBLF 1.4693680e-39 /* Min positive double */
X#define MAX_NEG_DBLF -1.7014118e38 /* Max negative double */
X#define MIN_NEG_DBLF -1.4693680e-39 /* Min negative double */
X#define MAX_EXPONENT 127 /* Max exponent allowed */
X#define RECIP_MIN 5.877471e-39 /* MAX_POS_DBLF >= 1/RECIP_MIN */
X#define RECIP_MAX 1.7014118e38 /* MIN_POS_DBLF <= 1/RECIP_MAX */
X#define LN_MAXPOSDBL 88.0 /* LN(MAX_POS_DBLF) */
X#define LN_MINPOSDBL -89.4 /* LN(MIN_POS_DBLF) */
X#define TANH_MAXARG 16 /* |TANH(maxarg)| = 1.0 */
X#define SQRT_MPDF 1.304380e19 /* SQRT(MAX_POS_DBLF) */
X#define X6_UNDERFLOWS 3.37174e-7 /* X**6 almost underflows */
X#define X16_UNDERFLOWS 3.74063e-3 /* X**16 almost underflows */
X#endif
X
X
X/*
X * MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X */
X
X#if defined(mc68000) && defined(IEEE)
X#define MIN_POS_DBLF 2.22507385850720220000e-308; /* Min positive dbl */
X#define MAX_POS_DBLF 1.79769313486231440000e+308; /* Max positive dbl */
X#define MIN_NEG_DBLF -2.22507385850720220000e-308; /* Min negative dbl */
X#define MAX_NEG_DBLF -1.79769313486231440000e+308; /* Max negative dbl */
X#define MAX_PEXPONENT 1024 /* Max pos exponent */
X#define MAX_NEXPONENT -1022 /* Max neg exponent */
X#endif
X
X
X#if defined(mc68000) && !defined(IEEE)
X#define MAX_POS_DBLF 1.7014118e38 /* Max positive double */
X#define MIN_POS_DBLF 1.4693680e-39 /* Min positive double */
X#define MAX_NEG_DBLF -1.7014118e38 /* Max negative double */
X#define MIN_NEG_DBLF -1.4693680e-39 /* Min negative double */
X#define MAX_EXPONENT 127 /* Max exponent allowed */
X#define RECIP_MIN 5.877471e-39 /* MAX_POS_DBLF >= 1/RECIP_MIN */
X#define RECIP_MAX 1.7014118e38 /* MIN_POS_DBLF <= 1/RECIP_MAX */
X#define LN_MAXPOSDBL 88.0 /* LN(MAX_POS_DBLF) */
X#define LN_MINPOSDBL -89.4 /* LN(MIN_POS_DBLF) */
X#define TANH_MAXARG 16 /* |TANH(maxarg)| = 1.0 */
X#define SQRT_MPDF 1.304380e19 /* SQRT(MAX_POS_DBLF) */
X#define X6_UNDERFLOWS 3.37174e-7 /* X**6 almost underflows */
X#define X16_UNDERFLOWS 3.74063e-3 /* X**16 almost underflows */
X#endif
X
X
X/* Define some commonly used constants */
X
X#define PI 3.14159265358979323846
X#define TWOPI (2.0 * PI)
X#define HALFPI (PI / 2.0)
X#define FOURTHPI (PI / 4.0)
X#define SIXTHPI (PI / 6.0)
X#define LOG2E 1.4426950408889634073
X#define LOG10E 0.4342944819032518276
X#define SQRT2 1.4142135623730950488
X#define SQRT3 1.7320508075688772935
X#define LN2 0.6931471805599453094
X#define LNSQRT2 0.3465735902799726547
X
X
X/* Configuration options */
X
X#define MAX_ERRORS 10 /* Limit on number of errors before abort */
END_OF_funcs/src/unused/OLDpml.h
if test 6798 -ne `wc -c <funcs/src/unused/OLDpml.h`; then
echo shar: \"funcs/src/unused/OLDpml.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/d2d.dat -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/d2d.dat\"
else
echo shar: Extracting \"tests/d2d.dat\" \(8975 characters\)
sed "s/^X//" >tests/d2d.dat <<'END_OF_tests/d2d.dat'
Xdabs 1.2345 1.2345
Xdabs -1.2345 1.2345
Xdabs 1.2345e+20 1.2345e+20
Xdabs -9.987e+25 9.987e+25
Xdabs -7.89e-24 7.89e-24
Xdabs 5.67e-15 5.67e-15
Xexp -40.0 4.248354255291e-18
Xexp -30.0 9.357622968840e-14
Xexp -20.0 2.061153622438e-09
Xexp -10.0 4.539992976248e-05
Xexp -2.0 1.353352832366e-01
Xexp -1.0 0.36787944117144232159552377016146087
Xexp -0.25 0.77880078307140486824517026697832065
Xexp 0.0 1.000000000000e00
Xexp 0.25 1.28402541668774148407342056806243646
Xexp 1.0 2.71828182845904523536028747135266250
Xexp 2.0 7.389056098930e00
Xexp 10.0 2.202646579480e04
Xexp 20.0 4.851651954097e08
Xexp 30.0 1.068647458152e13
Xexp 40.0 2.353852668370e17
Xsqrt 0.0 0.0
Xsqrt 4.999999999999999999e-02 2.236067977499789696e-01
Xsqrt 2.500000000000000000e-01 5.000000000000000000e-01
Xsqrt 5.000000000000000000e-01 7.071067811865475243e-01
Xsqrt 1.000000000000000000e00 1.000000000000000000e00
Xsqrt 2.000000000000000000e00 1.414213562373095048e00
Xsqrt 1.000000000000000000e01 3.162277660168379332e00
Xsqrt 3.000000000000000000e01 5.477225575051661134e00
Xsqrt 7.000000000000000000e01 8.366600265340755481e00
Xsqrt 2.000000000000000000e02 1.414213562373095049e01
Xsqrt 1.999999999999999999e32 1.414213562373095048e16
Xsqrt 0.6366197723675813430755350534900574 0.79788456080286535587989211986876373
Xsqrt 3.1415926535897932384626433832795028 1.77245385090551602729816748334114518
Xlog 4.999999999999999999e-02 -2.995732273553990993e00
Xlog 2.500000000000000000e-01 -1.386294361119890618e00
Xlog 5.000000000000000000e-01 -6.931471805599453094e-01
Xlog 1.000000000000000000e00 0.000000000000000000e00
Xlog 2.0 0.69314718055994530941723212145817657
Xlog 4.0 1.38629436111989061883446424291635313
Xlog 10.0 2.30258509299404568401799145468436421
Xlog 1.000000000000000000e01 2.302585092994045684e00
Xlog 3.000000000000000000e01 3.401197381662155375e00
Xlog 5.000000000000000000e01 3.912023005428146058e00
Xlog 2.000000000000000000e02 5.298317366548036677e00
Xlog 1.999999999999999999e32 7.437587015636940721e01
Xatanh -9.999899999999999999e-01 -6.103033822758835533e00
Xatanh -5.000000000000000000e-01 -5.493061443340548457e-01
Xatanh -2.500000000000000000e-01 -2.554128118829953416e-01
Xatanh 0.000000000000000000e00 0.000000000000000000e00
Xatanh 2.500000000000000000e-01 2.554128118829953415e-01
Xatanh 5.000000000000000000e-01 5.493061443340548457e-01
Xatanh 9.999899999999999999e-01 6.103033822758835533e00
Xsinh -8.000000000000000000e01 -2.770311192196755026e34
Xsinh -1.000000000000000000e01 -1.101323287470339337e04
Xsinh -5.000000000000000000e-01 -5.210953054937473615e-01
Xsinh 0.000000000000000000e00 0.000000000000000000e00
Xsinh 1.500000000000000000e00 2.129279455094817497e00
Xsinh 1.000000000000000000e01 1.101323287470339337e04
Xsinh 8.000000000000000000e01 2.770311192196755026e34
Xcosh -8.000000000000000000e01 2.770311192196755026e34
Xcosh -1.000000000000000000e01 1.101323292010332313e04
Xcosh -5.000000000000000000e-01 1.127625965206380785e00
Xcosh 0.000000000000000000e00 1.000000000000000000e00
Xcosh 1.500000000000000000e00 2.352409615243247325e00
Xcosh 1.000000000000000000e01 1.101323292010332313e04
Xcosh 8.000000000000000000e01 2.770311192196755026e34
Xtanh -8.000000000000000000e01 -1.000000000000000000e00
Xtanh -1.000000000000000000e01 -9.999999958776927635e-01
Xtanh -5.000000000000000000e-01 -4.621171572600097586e-01
Xtanh 0.000000000000000000e00 0.000000000000000000e00
Xtanh 1.500000000000000000e00 9.051482536448664383e-01
Xtanh 1.000000000000000000e01 9.999999958776927635e-01
Xtanh 8.000000000000000000e01 1.000000000000000000e00
Xsin 5.00e-11 5.000000000000000000e-11
Xsin 0.0 0.0
Xsin 0.25 0.24740395925452292959684870484938920
Xsin 0.5 0.47942553860420300027328793521557139
Xsin 0.78539816339744830961556084581987572 0.70710678118654752440084436210484903
Xsin -0.7853981633974483096155608458198757 -0.7071067811865475244008443621048490
Xsin 1.0 0.84147098480789650665250232163029900
Xsin 2.00e00 9.092974268256816953e-01
Xsin 2.50e00 5.984721441039564939e-01
Xsin 3.50e00 -3.507832276896198480e-01
Xsin 4.00e00 -7.568024953079282514e-01
Xsin 5.00e00 -9.589242746631384689e-01
Xsin 6.00e00 -2.794154981989258727e-01
Xsin 7.50e00 9.379999767747388579e-01
Xsin 8.00e00 9.893582466233817778e-01
Xsin 9.00e00 4.121184852417565697e-01
Xsin 9.50e00 -7.515112046180930728e-02
Xsin -5.00e-01 -4.794255386042030000e-01
Xsin -1.00e00 -8.414709848078965066e-01
Xsin -2.00e00 -9.092974268256816953e-01
Xsin -2.50e00 -5.984721441039564939e-01
Xsin -3.50e00 3.507832276896198480e-01
Xsin -4.00e00 7.568024953079282514e-01
Xsin -5.00e00 9.589242746631384689e-01
Xsin -5.50e00 7.055403255703919063e-01
Xsin -6.50e00 -2.151199880878155242e-01
Xsin -7.50e00 -9.379999767747388579e-01
Xsin -8.00e00 -9.893582466233817778e-01
Xsin -9.00e00 -4.121184852417565697e-01
Xsin -9.50e00 7.515112046180930728e-02
Xcos 5.00e-11 9.999999999999999997e-01
Xcos 0.00e00 9.999999999999999997e-01
Xcos 5.00e-01 8.77582561890372716e-01
Xcos 1.00e00 5.403023058681397173e-01
Xcos 2.00e00 -4.161468365471423870e-01
Xcos 2.50e00 -8.011436155469337148e-01
Xcos 3.50e00 -9.364566872907963376e-01
Xcos 4.00e00 -6.536436208636119144e-01
Xcos 5.00e00 2.836621854632262644e-01
Xcos 5.50e00 7.086697742912599999e-01
Xcos 6.50e00 9.765876257280234999e-01
Xcos 7.50e00 3.466353178350258109e-01
Xcos 8.00e00 -1.455000338086135258e-01
Xcos 9.00e00 -9.111302618846769882e-01
Xcos 9.50e00 -9.971721561963784728e-01
Xcos -5.00e-01 8.775825618903727160e-01
Xcos -1.00e00 5.403023058681397173e-01
Xcos -2.00e00 -4.161468365471423870e-01
Xcos -2.50e00 -8.011436155469337148e-01
Xcos -3.50e00 -9.364566872907963376e-01
Xcos -4.00e00 -6.536436208636119144e-01
Xcos -5.00e00 2.836621854632262644e-01
Xcos -6.00e00 9.601702866503660205e-01
Xcos -6.50e00 9.765876257280234999e-01
Xcos -7.50e00 3.466353178350258109e-01
Xcos -8.00e00 -1.455000338086135258e-01
Xcos -9.00e00 -9.111302618846769882e-01
Xcos -9.50e00 -9.971721561963784728e-01
Xtan -2.000000000000000000e00 2.185039863261518991e00
Xtan -1.000000000000000000e00 -1.557407724654902230e00
Xtan -5.000000000000000000e-01 -5.463024898437905132e-01
Xtan 0.000000000000000000e00 0.000000000000000000e00
Xtan 0.39269908169872415480783042290993786 0.41421356237309504880168872420969807
Xtan 0.25 0.25534192122103626650448223649047368
Xtan 0.5 0.54630248984379051325517946578028538
Xtan 0.78539816339744830961556084581987572 1.0
Xtan -0.78539816339744830961556084581987572 -1.0
Xtan 1.0 1.55740772465490223050697480745836017
Xtan 2.000000000000000000e00 -2.185039863261518991e00
Xlog10 1.00e-10 -1.000000000000000000e01
Xlog10 1.00e-05 -5.000000000000000000e00
Xlog10 5.00e-01 -3.010299956639811952e-01
Xlog10 1.00e00 0.000000000000000000e00
Xlog10 1.50e00 1.760912590556812420e-01
Xlog10 2.0 0.30102999566398119521373889472449303
Xlog10 2.71828182845904523536028747135266250 0.43429448190325182765112891891660508
Xlog10 1.00e05 5.000000000000000000e00
Xlog10 1.00e20 2.000000000000000000e01
Xatan 0.000000000000e00 0.000000000000e00
Xatan 4.999999999999e-02 4.995839572194e-02
Xatan 2.500000000000e-01 2.449786631268e-01
Xatan 3.500000000000e-01 3.366748193867e-01
Xatan 7.000000000000e-01 6.107259643892e-01
Xatan 1.000000000000e00 7.85398163397e-01
Xatan 2.000000000000e00 1.10714871779e00
Xatan 1.400000000000e01 1.49948886200e00
Xatan 3.000000000000e01 1.53747533091e00
Xatan 5.000000000000e01 1.55079899282e00
Xatan 7.000000000000e01 1.55651158420e00
Xatan 9.000000000000e01 1.55968567289e00
Xatan 2.000000000000e02 1.56579636846e00
Xasinh -8.000000000000000000e01 -5.075212875445203927e00
Xasinh -1.000000000000000000e01 -2.998222950297969704e00
Xasinh -5.000000000000000000e-01 -4.812118250596034474e-01
Xasinh 0.000000000000000000e00 0.000000000000000000e00
Xasinh 1.500000000000000000e00 1.194763217287109304e00
Xasinh 1.000000000000000000e01 2.998222950297969738e00
Xasinh 8.000000000000000000e01 5.075212875445207223e00
Xasin -1.000000000000000000e00 -1.570796326794896619e00
Xasin -7.000000000000000000e-01 -7.753974966107530636e-01
Xasin -1.999999999999999999e-01 -2.013579207903307914e-01
Xasin 0.000000000000000000e00 0.000000000000000000e00
Xasin 1.999999999999999999e-01 2.013579207903307914e-01
Xasin 7.000000000000000000e-01 7.753974966107530636e-01
Xasin 1.000000000000000000e00 1.570796326794896619e00
Xacosh 1.000000000000000000e00 0.000000000000000000e00
Xacosh 2.000000000000000000e00 1.316957896924816708e00
Xacosh 1.000000000000000000e01 2.993222846126380897e00
Xacosh 2.000000000000000000e02 5.991458297049387423e00
Xacosh 1.000000000000000000e03 7.600902209541988611e00
Xacosh 1.000000000000000000e05 1.220607264550517372e01
Xacosh 1.000000000000000000e18 4.213967885445276762e01
Xacos -1.000000000000000000e00 3.141592653589793238e00
Xacos -7.000000000000000000e-01 2.346193823405649682e00
Xacos -1.999999999999999999e-01 1.772154247585227410e00
Xacos 0.000000000000000000e00 1.570796326794896619e00
Xacos 1.999999999999999999e-01 1.369438406004565827e00
Xacos 7.000000000000000000e-01 7.953988301841435554e-01
Xacos 1.000000000000000000e00 0.000000000000000000e00
END_OF_tests/d2d.dat
if test 8975 -ne `wc -c <tests/d2d.dat`; then
echo shar: \"tests/d2d.dat\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/errors.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/errors.c\"
else
echo shar: Extracting \"tests/errors.c\" \(7459 characters\)
sed "s/^X//" >tests/errors.c <<'END_OF_tests/errors.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * testerrors.c error handler test for portable math library
X *
X * KEY WORDS
X *
X * test routines
X * portable math library
X * error handler
X *
X * DESCRIPTION
X *
X * Tests functions of the optional error handler for the portable
X * math library. Tests the COUNT, LOG, and CONTINUE bits for
X * each function, along with the task error limit.
X *
X * USAGE
X *
X * testerrors
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Engineering Software Tools
X * P.O. Box 2035
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "../funcs/pml.h"
X
Xchar *test1[] = {
X "\n****** TEST #1 ******",
X "",
X " Generate all errors",
X "",
X " (1) LOG bits, COUNT bits, and CONTINUE bits set.",
X " (2) Error limit set high enough to avoid aborting.",
X "",
X " This test should generate an error message for each",
X " error tested, and the final error count should be printed",
X " when done.",
X "\n",
X 0
X};
X
Xchar *test2[] = {
X "\n****** TEST #2 ******",
X "",
X " Turn off logging.",
X "",
X " (1) LOG bits reset.",
X " (2) COUNT bits and CONTINUE bits set.",
X " (3) Error limit set high enough to avoid aborting.",
X "",
X " This test should not log any error messages but the error",
X " count should be the same as test 1.",
X "\n",
X 0
X};
X
Xchar *test3[] = {
X "\n****** TEST #3 ******",
X "",
X " Turn off error counting",
X "",
X " (1) LOG and COUNT bits reset.",
X " (2) CONTINUE bits set.",
X " (3) Error limit set to 0.",
X "",
X " This test should not log or count any errors. The final",
X " error count should be zero.",
X "\n",
X 0
X};
X
Xchar *test4[] = {
X "\n****** TEST #4 ******",
X "",
X " Enable error limit and abort.",
X "",
X " (1) LOG bits, COUNT bits, and CONTINUE bits set.",
X " (2) Error limit set to 5.",
X "",
X " This test should abort after logging the sixth error.",
X "\n",
X 0
X};
X
X
Xmain()
X{
X initialize();
X pmllim(40);
X prtdoc(test1);
X gen_errors();
X printf("\nTotal errors counted = %d\n",pmlcnt());
X
X log_off();
X pmllim(40);
X prtdoc(test2);
X gen_errors();
X printf("\nTotal errors counted = %d\n",pmlcnt());
X
X count_off();
X pmllim(0);
X prtdoc(test3);
X gen_errors();
X printf("\nTotal errors counted = %d\n",pmlcnt());
X
X initialize();
X pmllim(5);
X prtdoc(test4);
X gen_errors();
X}
X
X
Xprtdoc(dp)
Xregister char **dp;
X{
X while (*dp) {
X printf("%s\n",*dp++);
X }
X}
X
X
Xinitialize()
X{
X pmlsfs(EXP_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(EXP_UNDERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(SCALE_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(SCALE_UNDERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(NEG_SQRT,LOG|COUNT|CONTINUE);
X pmlsfs(LOG_OF_ZERO,LOG|COUNT|CONTINUE);
X pmlsfs(LOG_OF_NEGATIVE,LOG|COUNT|CONTINUE);
X pmlsfs(ACOS_BADARG,LOG|COUNT|CONTINUE);
X pmlsfs(ASIN_BADARG,LOG|COUNT|CONTINUE);
X pmlsfs(TAN_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(COSH_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(COSH_UNDERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(SINH_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(SINH_UNDERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(ASINH_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(ACOSH_BADARG,LOG|COUNT|CONTINUE);
X pmlsfs(ACOSH_OVERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(ATANH_BADARG,LOG|COUNT|CONTINUE);
X pmlsfs(ATAN_UNDERFLOW,LOG|COUNT|CONTINUE);
X pmlsfs(C_DIV_ZERO,LOG|COUNT|CONTINUE);
X pmlsfs(CRCP_OF_ZERO,LOG|COUNT|CONTINUE);
X pmlsfs(DINT_2BIG,LOG|COUNT|CONTINUE);
X}
X
Xlog_off()
X{
X pmlcfs(EXP_OVERFLOW,LOG);
X pmlcfs(EXP_UNDERFLOW,LOG);
X pmlcfs(SCALE_OVERFLOW,LOG);
X pmlcfs(SCALE_UNDERFLOW,LOG);
X pmlcfs(NEG_SQRT,LOG);
X pmlcfs(LOG_OF_ZERO,LOG);
X pmlcfs(LOG_OF_NEGATIVE,LOG);
X pmlcfs(ACOS_BADARG,LOG);
X pmlcfs(ASIN_BADARG,LOG);
X pmlcfs(TAN_OVERFLOW,LOG);
X pmlcfs(COSH_OVERFLOW,LOG);
X pmlcfs(COSH_UNDERFLOW,LOG);
X pmlcfs(SINH_OVERFLOW,LOG);
X pmlcfs(SINH_UNDERFLOW,LOG);
X pmlcfs(ASINH_OVERFLOW,LOG);
X pmlcfs(ACOSH_BADARG,LOG);
X pmlcfs(ACOSH_OVERFLOW,LOG);
X pmlcfs(ATANH_BADARG,LOG);
X pmlcfs(ATAN_UNDERFLOW,LOG);
X pmlcfs(C_DIV_ZERO,LOG);
X pmlcfs(CRCP_OF_ZERO,LOG);
X pmlcfs(DINT_2BIG,LOG);
X}
X
Xcount_off()
X{
X pmlcfs(EXP_OVERFLOW,COUNT);
X pmlcfs(EXP_UNDERFLOW,COUNT);
X pmlcfs(SCALE_OVERFLOW,COUNT);
X pmlcfs(SCALE_UNDERFLOW,COUNT);
X pmlcfs(NEG_SQRT,COUNT);
X pmlcfs(LOG_OF_ZERO,COUNT);
X pmlcfs(LOG_OF_NEGATIVE,COUNT);
X pmlcfs(ACOS_BADARG,COUNT);
X pmlcfs(ASIN_BADARG,COUNT);
X pmlcfs(TAN_OVERFLOW,COUNT);
X pmlcfs(COSH_OVERFLOW,COUNT);
X pmlcfs(COSH_UNDERFLOW,COUNT);
X pmlcfs(SINH_OVERFLOW,COUNT);
X pmlcfs(SINH_UNDERFLOW,COUNT);
X pmlcfs(ASINH_OVERFLOW,COUNT);
X pmlcfs(ACOSH_BADARG,COUNT);
X pmlcfs(ACOSH_OVERFLOW,COUNT);
X pmlcfs(ATANH_BADARG,COUNT);
X pmlcfs(ATAN_UNDERFLOW,COUNT);
X pmlcfs(C_DIV_ZERO,COUNT);
X pmlcfs(CRCP_OF_ZERO,COUNT);
X pmlcfs(DINT_2BIG,COUNT);
X}
X
Xgen_errors()
X{
X complex z1, z2;
X
X z1.r = 1.0;
X z1.i = 0.0;
X z2.r = 0.0;
X z2.i = 0.0;
X printf("Testing for exp overflow.\n");
X exp(100.0);
X printf("Testing for exp underflow.\n");
X exp(-100.0);
X printf("Testing for scale exponent overflow.\n");
X scale(1.0,500);
X printf("Testing for scale exponent underflow.\n");
X scale(1.0,-500);
X printf("Testing for sqrt of negative argument.\n");
X sqrt(-1.0);
X printf("Testing for ln of zero.\n");
X ln(0.0);
X printf("Testing for ln of negative argument.\n");
X ln(-1.0);
X printf("Testing for acos argument magnitude greater than 1.0\n");
X acos(2.0);
X printf("Testing for asin argument magnitude greater than 1.0\n");
X asin(-2.0);
X printf("Testing for tan overflow\n");
X tan(HALFPI);
X printf("Testing for cosh overflow\n");
X cosh(LN_MAXPOSDBL+1.0);
X printf("Testing for cosh underflow\n");
X cosh(LN_MINPOSDBL-1.0);
X printf("Testing for sinh overflow\n");
X sinh(LN_MAXPOSDBL+1.0);
X printf("Testing for sinh underflow\n");
X sinh(LN_MINPOSDBL-1.0);
X printf("Testing for asinh overflow\n");
X asinh(2.0 * SQRT_MPDF);
X printf("Testing for acosh argument less than 1.0\n");
X acosh(0.0);
X printf("Testing for acosh overflow\n");
X acosh(2.0 * SQRT_MPDF);
X printf("Testing for atanh argument magnitude >= 1.0\n");
X atanh(-1.0);
X printf("Testing for atan underflow\n");
X atan(RECIP_MAX);
X printf("Testing for complex division by zero\n");
X cdiv(&z1,&z2);
X printf("Testing for complex reciprocal of zero\n");
X crcp(&z2);
X printf("Testing for dint argument has no fractional part\n");
X dint(MAX_POS_DBLF);
X}
X
END_OF_tests/errors.c
if test 7459 -ne `wc -c <tests/errors.c`; then
echo shar: \"tests/errors.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 4 \(of 6\).
cp /dev/null ark4isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 5 (of 6)."
# Contents: tests/c2d.c tests/d2d.c tests/dd2d.c tests/unused/d2i.c
# tests/unused/di2d.c
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:47 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f tests/c2d.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/c2d.c\"
else
echo shar: Extracting \"tests/c2d.c\" \(10177 characters\)
sed "s/^X//" >tests/c2d.c <<'END_OF_tests/c2d.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * c2d.c test complex to double math functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect a single
X * double precision complex argument and return a double.
X *
X * Most of the test data in the current data file (c2d.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first. Naturally, an error in a lower level
X * function will cause propagation of errors up to higher
X * level functions.
X *
X * USAGE
X *
X * c2d [-esv] [-l limit]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -l => report errors greater than
X * specified limit (default 10**-6)
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (c2d.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit.
X * Each test directive has the form:
X *
X * <name> <real(arg)> <imag(arg)> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * The first field, "name", is the name of the function
X * to test. The second field is the real part of the argument.
X * The third field is the imaginary part of the argument.
X * The last field is the expected result.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X
X#ifndef NODBG
X#include <dbg.h>
X#else
X#include "dbg.h"
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xstatic double max_abs_err = MAX_ABS_ERR;
X
Xextern double cabs (); /* Complex absolute value */
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X double (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "cabs", cabs, 0.0, /* Complex magnitude */
X NULL, NULL, 0.0 /* Function list end marker */
X};
X
X
X/*
X * FUNCTION
X *
X * main entry point for c2d test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain (argc, argv)
Xint argc;
Xchar *argv[];
X{
X ENTER ("main");
X DEBUGWHO (argv[0]);
X options (argc, argv);
X dotests (argv);
X statistics ();
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default real part of argument to 0.0
X * Default imaginary part of argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
X
Xdotests (argv)
Xchar *argv[];
X{
X char buffer[256]; /* Directive buffer */
X char function[64]; /* Specified function name */
X COMPLEX argument; /* Specified function argument */
X double expected; /* Specified expected result */
X double result; /* Actual result */
X double error; /* Relative or absolute error */
X double abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup (); /* Returns function test pointer */
X register char *strp; /* Pointer to next token in string */
X extern char *strtok ();
X extern double atof ();
X
X ENTER ("dotests");
X while (fgets (buffer, sizeof(buffer), stdin) != NULL) {
X strcpy (function, "{null}");
X argument.real = 0.0;
X argument.imag = 0.0;
X expected = 0.0;
X sscanf (buffer, "%s %le %le %le",
X function, &argument.real, &argument.imag, &expected);
X testp = lookup (function);
X if (testp == NULL) {
X fprintf (stderr, "%s: unknown function \"%s\".\n",
X argv[0], function);
X } else {
X result = (*testp -> func)(argument);
X if (vflag) {
X printf ("%s(%le + j %le) = %30.23le.\n",
X function, argument.real, argument.imag, result);
X }
X if (expected != 0.0) {
X error = (result - expected) / expected;
X } else {
X error = result;
X }
X if (error < 0.0) {
X abs_err = -error;
X } else {
X abs_err = error;
X }
X if ((abs_err > max_abs_err) || eflag) {
X fprintf (stderr, "%s: error in \"%s\"\n", argv[0], function);
X fprintf (stderr, "\treal (arg)\t%25.20le\n", argument.real);
X fprintf (stderr, "\timag (arg)\t%25.20le\n", argument.imag);
X fprintf (stderr, "\tresult\t\t%25.20le\n", result);
X fprintf (stderr, "\texpected\t%25.20le\n", expected);
X }
X if (abs_err > testp -> max_err) {
X testp -> max_err = abs_err;
X }
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
X
Xoptions (argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int flag;
X extern int getopt ();
X extern char *optarg;
X
X ENTER ("options");
X eflag = sflag = vflag = FALSE;
X while ((flag = getopt (argc, argv, "#:el:sv")) != EOF) {
X switch (flag) {
X case '#':
X DEBUGPUSH (optarg);
X break;
X case 'e':
X eflag = TRUE;
X break;
X case 'l':
X sscanf (optarg, "%le", &max_abs_err);
X DEBUG3 ("args", "max_abs_err = %le", max_abs_err);
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup (funcname)
Xchar *funcname;
X{
X struct test *testp;
X struct test *rtnval;
X
X ENTER ("lookup");
X rtnval = (struct test *) NULL;
X for (testp = tests; testp -> name != NULL && rtnval == NULL; testp++) {
X if (!strcmp (testp -> name, funcname)) {
X rtnval = testp;
X }
X }
X LEAVE ();
X return (rtnval);
X}
X
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics ()
X{
X struct test *tp;
X
X ENTER ("statistics");
X if (sflag) {
X for (tp = tests; tp -> name != NULL; tp++) {
X printf ("%s:\tmaximum relative error %le\n",
X tp -> name, tp -> max_err);
X }
X }
X LEAVE ();
X}
END_OF_tests/c2d.c
if test 10177 -ne `wc -c <tests/c2d.c`; then
echo shar: \"tests/c2d.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/d2d.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/d2d.c\"
else
echo shar: Extracting \"tests/d2d.c\" \(11104 characters\)
sed "s/^X//" >tests/d2d.c <<'END_OF_tests/d2d.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FILE
X *
X * d2d.c test double to double math functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect a single
X * double precision argument and return a double.
X *
X * Most of the test data in the current data file (d2d.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first. Naturally, an error in a lower level
X * function will cause propagation of errors up to higher
X * level functions.
X *
X * USAGE
X *
X * d2d [-esv] [-l limit]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -l => report errors greater than
X * specified limit (default 10**-6)
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (d2d.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit.
X * Each test directive has the form:
X *
X * <name> <argument> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * The first field, "name", is the name of the function
X * to test (sqrt, ln, exp, etc). The second field
X * is the argument to use in calling the specified function.
X * The third field is the expected result.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X */
X
X
X#include <stdio.h>
X#include "pml.h"
X
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xstatic double max_abs_err = MAX_ABS_ERR;
X
Xextern double dabs ();
Xextern double acos ();
Xextern double acosh ();
Xextern double asin ();
Xextern double asinh ();
Xextern double atan ();
Xextern double atanh ();
Xextern double cos ();
Xextern double cosh ();
Xextern double exp ();
Xextern double log ();
Xextern double log10 ();
Xextern double sin ();
Xextern double sinh ();
Xextern double sqrt ();
Xextern double tan ();
Xextern double tanh ();
X
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X double (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "dabs", dabs, 0.0, /* Absolute value */
X "acos", acos, 0.0, /* Arc cosine (in radians) */
X "acosh", acosh, 0.0, /* Hyperbolic arc cosine (in radians) */
X "asin", asin, 0.0, /* Arc sine (in radians) */
X "asinh", asinh, 0.0, /* Hyperbolic arc sine (in radians) */
X "atan", atan, 0.0, /* Arc tangent (in radians) */
X "atanh", atanh, 0.0, /* Hyperbolic arc tangent (in radians) */
X "cos", cos, 0.0, /* Cosine (argument in radians) */
X "cosh", cosh, 0.0, /* Hyperbolic cosine (argument in radians) */
X "exp", exp, 0.0, /* Exponential */
X "log", log, 0.0, /* Natural logarithm */
X "log10", log10, 0.0, /* Log to base 10 */
X "sin", sin, 0.0, /* Sine (argument in radians) */
X "sinh", sinh, 0.0, /* Hyperbolic sine (argument in radians) */
X "sqrt", sqrt, 0.0, /* Garden variety square root */
X "tan", tan, 0.0, /* Tangent (argument in radians) */
X "tanh", tanh, 0.0, /* Hyperbolic tangent (argument in radians) */
X NULL, NULL, 0.0 /* Function list end marker */
X};
X
X
X/*
X * FUNCTION
X *
X * main entry point for d2d test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain (argc, argv)
Xint argc;
Xchar *argv[];
X{
X VOID dotests ();
X
X DBUG_ENTER ("main");
X DBUG_PROCESS (argv[0]);
X options (argc, argv);
X dotests (argv);
X statistics ();
X DBUG_RETURN (0);
X}
X
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
X
XVOID dotests (argv)
Xchar *argv[];
X{
X char buffer[256]; /* Directive buffer */
X char function[64]; /* Specified function name */
X double argument; /* Specified function argument */
X double expected; /* Specified expected result */
X double result; /* Actual result */
X double error; /* Relative or absolute error */
X double abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup (); /* Returns function test pointer */
X register char *strp; /* Pointer to next token in string */
X extern char *strtok ();
X extern double atof ();
X
X DBUG_ENTER ("dotests");
X while (fgets (buffer, sizeof(buffer), stdin) != NULL) {
X strcpy (function, "{null}");
X argument = 0.0;
X expected = 0.0;
X sscanf (buffer, "%s %le %le", function, &argument, &expected);
X testp = lookup (function);
X if (testp == NULL) {
X fprintf (stderr, "%s: unknown function \"%s\".\n",
X argv[0], function);
X } else {
X result = (*testp -> func)(argument);
X if (vflag) {
X printf ("%s(%le) = %30.23le.\n", function, argument, result);
X }
X if (expected != 0.0) {
X error = (result - expected) / expected;
X } else {
X error = result;
X }
X if (error < 0.0) {
X abs_err = -error;
X } else {
X abs_err = error;
X }
X if ((abs_err > max_abs_err) || eflag) {
X fprintf (stderr, "%s: error in \"%s\"\n", argv[0], function);
X fprintf (stderr, "\targument\t%25.20le\n", argument);
X fprintf (stderr, "\tresult\t\t%25.20le\n", result);
X fprintf (stderr, "\texpected\t%25.20le\n", expected);
X }
X if (abs_err > testp -> max_err) {
X testp -> max_err = abs_err;
X }
X }
X }
X DBUG_VOID_RETURN;
X}
X
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
X
Xoptions (argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int flag;
X extern int getopt ();
X extern char *optarg;
X
X DBUG_ENTER ("options");
X eflag = sflag = vflag = FALSE;
X while ((flag = getopt (argc, argv, "#:el:sv")) != EOF) {
X switch (flag) {
X case '#':
X DBUG_PUSH (optarg);
X break;
X case 'e':
X eflag = TRUE;
X break;
X case 'l':
X sscanf (optarg, "%le", &max_abs_err);
X DBUG_3 ("args", "max_abs_err = %le", max_abs_err);
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X }
X }
X DBUG_VOID_RETURN;
X}
X
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup (funcname)
Xchar *funcname;
X{
X struct test *testp;
X struct test *rtnval;
X
X DBUG_ENTER ("lookup");
X rtnval = (struct test *) NULL;
X for (testp = tests; testp -> name != NULL && rtnval == NULL; testp++) {
X if (!strcmp (testp -> name, funcname)) {
X rtnval = testp;
X }
X }
X DBUG_RETURN (rtnval);
X}
X
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics ()
X{
X struct test *tp;
X
X DBUG_ENTER ("statistics");
X if (sflag) {
X for (tp = tests; tp -> name != NULL; tp++) {
X printf ("%s:\tmaximum relative error %le\n",
X tp -> name, tp -> max_err);
X }
X }
X DBUG_VOID_RETURN;
X}
END_OF_tests/d2d.c
if test 11104 -ne `wc -c <tests/d2d.c`; then
echo shar: \"tests/d2d.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/dd2d.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/dd2d.c\"
else
echo shar: Extracting \"tests/dd2d.c\" \(10222 characters\)
sed "s/^X//" >tests/dd2d.c <<'END_OF_tests/dd2d.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * dd2d.c test dual double to double math functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect two
X * double precision arguments and return a double.
X *
X * Most of the test data in the current data file (dd2d.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first. Naturally, an error in a lower level
X * function will cause propagation of errors up to higher
X * level functions.
X *
X * USAGE
X *
X * dd2d [-esv] [-l limit]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -l => report errors greater than
X * specified limit (default 10**-6)
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (dd2d.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit.
X * Each test directive has the form:
X *
X * <name> <arg1> <arg2> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * The first field, "name", is the name of the function
X * to test (sqrt, ln, exp, etc). The second and third fields
X * are the arguments to use in calling the specified function.
X * The last field is the expected result.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
X
X#include <stdio.h>
X
X#ifndef NODBG
X#include <dbg.h>
X#else
X#include "dbg.h"
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xstatic double max_abs_err = MAX_ABS_ERR;
X
Xextern double atan2 ();
Xextern double max ();
Xextern double min ();
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X double (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "atan2", atan2, 0.0, /* Arc tangent with two args */
X "max", max, 0.0, /* Maximum value */
X "min", min, 0.0, /* Minimum value */
X NULL, NULL, 0.0 /* Function list end marker */
X};
X
X
X/*
X * FUNCTION
X *
X * main entry point for dd2d test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain (argc, argv)
Xint argc;
Xchar *argv[];
X{
X ENTER ("main");
X DEBUGWHO (argv[0]);
X options (argc, argv);
X dotests (argv);
X statistics ();
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default argument 1 to 0.0
X * Default argument 2 to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
X
Xdotests (argv)
Xchar *argv[];
X{
X char buffer[256]; /* Directive buffer */
X char function[64]; /* Specified function name */
X double arg1; /* Specified function argument 1 */
X double arg2; /* Specified function argument 2 */
X double expected; /* Specified expected result */
X double result; /* Actual result */
X double error; /* Relative or absolute error */
X double abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup (); /* Returns function test pointer */
X register char *strp; /* Pointer to next token in string */
X extern char *strtok ();
X extern double atof ();
X
X ENTER ("dotests");
X while (fgets (buffer, sizeof(buffer), stdin) != NULL) {
X strcpy (function, "{null}");
X arg1 = 0.0;
X arg2 = 0.0;
X expected = 0.0;
X sscanf (buffer, "%s %le %le %le", function, &arg1, &arg2, &expected);
X testp = lookup (function);
X if (testp == NULL) {
X fprintf (stderr, "%s: unknown function \"%s\".\n",
X argv[0], function);
X } else {
X result = (*testp -> func)(arg1, arg2);
X if (vflag) {
X printf ("%s(%le,%le) = %30.23le.\n",
X function, arg1, arg2, result);
X }
X if (expected != 0.0) {
X error = (result - expected) / expected;
X } else {
X error = result;
X }
X if (error < 0.0) {
X abs_err = -error;
X } else {
X abs_err = error;
X }
X if ((abs_err > max_abs_err) || eflag) {
X fprintf (stderr, "%s: error in \"%s\"\n", argv[0], function);
X fprintf (stderr, "\targument 1\t%25.20le\n", arg1);
X fprintf (stderr, "\targument 2\t%25.20le\n", arg2);
X fprintf (stderr, "\tresult\t\t%25.20le\n", result);
X fprintf (stderr, "\texpected\t%25.20le\n", expected);
X }
X if (abs_err > testp -> max_err) {
X testp -> max_err = abs_err;
X }
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
X
Xoptions (argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int flag;
X extern int getopt ();
X extern char *optarg;
X
X ENTER ("options");
X eflag = sflag = vflag = FALSE;
X while ((flag = getopt (argc, argv, "#:el:sv")) != EOF) {
X switch (flag) {
X case '#':
X DEBUGPUSH (optarg);
X break;
X case 'e':
X eflag = TRUE;
X break;
X case 'l':
X sscanf (optarg, "%le", &max_abs_err);
X DEBUG3 ("args", "max_abs_err = %le", max_abs_err);
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup (funcname)
Xchar *funcname;
X{
X struct test *testp;
X struct test *rtnval;
X
X ENTER ("lookup");
X rtnval = (struct test *) NULL;
X for (testp = tests; testp -> name != NULL && rtnval == NULL; testp++) {
X if (!strcmp (testp -> name, funcname)) {
X rtnval = testp;
X }
X }
X LEAVE ();
X return (rtnval);
X}
X
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics ()
X{
X struct test *tp;
X
X ENTER ("statistics");
X if (sflag) {
X for (tp = tests; tp -> name != NULL; tp++) {
X printf ("%s:\tmaximum relative error %le\n",
X tp -> name, tp -> max_err);
X }
X }
X LEAVE ();
X}
END_OF_tests/dd2d.c
if test 10222 -ne `wc -c <tests/dd2d.c`; then
echo shar: \"tests/dd2d.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/unused/d2i.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/unused/d2i.c\"
else
echo shar: Extracting \"tests/unused/d2i.c\" \(9030 characters\)
sed "s/^X//" >tests/unused/d2i.c <<'END_OF_tests/unused/d2i.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * d2i.c test portable math library functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect a single
X * double precision argument and return an integer.
X *
X * Most of the test data in the current data file (d2i.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first.
X *
X * USAGE
X *
X * d2i [-esv]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (d2i.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit (specified
X * at compile time). Each test directive has the form:
X *
X * <name> <argument> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * The first field, "name", is the name of the function
X * to test (sqrt, ln, exp, etc). The second field
X * is the argument to use in calling the specified function.
X * The third field is the expected result.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Engineering Software Tools
X * P.O. Box 2035
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
X#include <stdio.h>
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to check error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xextern int xexp();
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X int (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "xexp", xexp, 0.0, /* Extract exponent */
X NULL,NULL,0.0, /* Function list end marker */
X};
X
X/*
X * FUNCTION
X *
X * main entry point for d2i test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain(argc,argv)
Xint argc;
Xchar *argv[];
X{
X options(argv);
X dotests(argv);
X statistics();
X}
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
Xdotests(argv)
Xchar *argv[];
X{
X char buffer[128]; /* Directive buffer */
X char function[16]; /* Specified function name */
X double argument; /* Specified function argument */
X int expected; /* Specified expected result */
X int result; /* Actual result */
X int error; /* Relative or absolute error */
X int abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup(); /* Returns function test pointer */
X
X while (fgets(buffer,sizeof(buffer),stdin) != NULL) {
X strcpy(function,"{null}");
X argument = 0.0;
X expected = 0;
X sscanf(buffer,"%s%F%d",&function[0],&argument,&expected);
X testp = lookup(function);
X if (testp == NULL) {
X fprintf(stderr,"%s: unknown function \"%s\".\n",argv[0],function);
X } else {
X result = (*testp->func)(argument);
X if (vflag) {
X printf("%s(%e) = %d\n",function,argument,result);
X }
X if (expected != 0) {
X error = (result - expected) / expected;
X } else {
X error = result;
X }
X if (error < 0) {
X abs_err = -error;
X } else {
X abs_err = error;
X }
X if ((abs_err > MAX_ABS_ERR) || eflag) {
X fprintf(stderr,"%s: %s(%e) = %d; should be %d\n",
X argv[0],function,argument,result,expected);
X }
X if (abs_err > testp->max_err) {
X testp->max_err = abs_err;
X }
X }
X }
X}
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
Xstatic char *error1[] = {
X "%s: unknown flag '%c'\n"
X};
X
Xoptions(argv)
Xchar *argv[];
X{
X char *argpntr;
X
X eflag = sflag = vflag = FALSE;
X argpntr = argv[1];
X if (argpntr != NULL) {
X if (*argpntr = '-') {
X while (*++argpntr != NULL) {
X switch (*argpntr) {
X case 'e':
X eflag = TRUE;
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X default:
X fprintf(stderr,error1[0],argv[0],*argpntr);
X break;
X }
X }
X }
X }
X}
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup(funcname)
Xchar *funcname;
X{
X struct test *testp;
X
X for (testp = tests; testp->name != NULL; testp++) {
X if (!strcmp(testp->name,funcname)) {
X return(testp);
X }
X }
X return((struct test *)NULL);
X}
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics()
X{
X struct test *tp;
X
X if (sflag) {
X for (tp = tests; tp->name != NULL; tp++) {
X printf("%s:\tmaximum relative error %e\n",tp->name,tp->max_err);
X }
X }
X}
END_OF_tests/unused/d2i.c
if test 9030 -ne `wc -c <tests/unused/d2i.c`; then
echo shar: \"tests/unused/d2i.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/unused/di2d.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/unused/di2d.c\"
else
echo shar: Extracting \"tests/unused/di2d.c\" \(9150 characters\)
sed "s/^X//" >tests/unused/di2d.c <<'END_OF_tests/unused/di2d.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * di2d.c test some portable math library functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect a double
X * precision argument, followed by an integer argument, and
X * return a double.
X *
X * Most of the test data in the current data file (di2d.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first.
X *
X * USAGE
X *
X * di2d [-esv]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (di2d.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit (specified
X * at compile time). Each test directive has the form:
X *
X * <name> <arg1> <arg2> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * The first field, "name", is the name of the function
X * to test (sqrt, ln, exp, etc). The second field
X * is the argument to use in calling the specified function.
X * The third field is the expected result.
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Engineering Software Tools
X * P.O. Box 2035
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
X#include <stdio.h>
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to check error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xextern double scale();
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X double (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "scale", scale, 0.0, /* Scale a double by specified exponent */
X NULL,NULL,0.0, /* Function list end marker */
X};
X
X/*
X * FUNCTION
X *
X * main entry point for di2d test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain(argc,argv)
Xint argc;
Xchar *argv[];
X{
X options(argv);
X dotests(argv);
X statistics();
X}
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
Xdotests(argv)
Xchar *argv[];
X{
X char buffer[128]; /* Directive buffer */
X char function[16]; /* Specified function name */
X double arg1; /* Specified function argument */
X int arg2; /* Specified function argument #2 */
X double expected; /* Specified expected result */
X double result; /* Actual result */
X double error; /* Relative or absolute error */
X double abs_err, abs(); /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup(); /* Returns function test pointer */
X
X while (fgets(buffer,sizeof(buffer),stdin) != NULL) {
X strcpy(function,"{null}");
X arg1 = 0.0;
X arg2 = 0;
X expected = 0.0;
X sscanf(buffer,"%s%F%d%F",&function[0],&arg1,&arg2,&expected);
X testp = lookup(function);
X if (testp == NULL) {
X fprintf(stderr,"%s: unknown function \"%s\".\n",argv[0],function);
X } else {
X result = (*testp->func)(arg1,arg2);
X if (vflag) {
X printf("%s(%e,%d) = %25.18e.\n",function,arg1,arg2,result);
X }
X if (expected != 0.0) {
X error = (result - expected) / expected;
X } else {
X error = result;
X }
X abs_err = abs(error);
X if ((abs_err > MAX_ABS_ERR) || eflag) {
X fprintf(stderr,"%s: %s(%e,%d) = %20.15e; should be %20.15e.\n",
X argv[0],function,arg1,arg2,result,expected);
X }
X if (abs_err > testp->max_err) {
X testp->max_err = abs_err;
X }
X }
X }
X}
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
Xstatic char *error1[] = {
X "%s: unknown flag '%c'\n"
X};
X
Xoptions(argv)
Xchar *argv[];
X{
X char *argpntr;
X
X eflag = sflag = vflag = FALSE;
X argpntr = argv[1];
X if (argpntr != NULL) {
X if (*argpntr = '-') {
X while (*++argpntr != NULL) {
X switch (*argpntr) {
X case 'e':
X eflag = TRUE;
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X default:
X fprintf(stderr,error1[0],argv[0],*argpntr);
X break;
X }
X }
X }
X }
X}
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup(funcname)
Xchar *funcname;
X{
X struct test *testp;
X
X for (testp = tests; testp->name != NULL; testp++) {
X if (!strcmp(testp->name,funcname)) {
X return(testp);
X }
X }
X return((struct test *)NULL);
X}
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics()
X{
X struct test *tp;
X
X if (sflag) {
X for (tp = tests; tp->name != NULL; tp++) {
X printf("%s:\tmaximum relative error %e\n",tp->name,tp->max_err);
X }
X }
X}
END_OF_tests/unused/di2d.c
if test 9150 -ne `wc -c <tests/unused/di2d.c`; then
echo shar: \"tests/unused/di2d.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 5 \(of 6\).
cp /dev/null ark5isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =fnf@mcdsun.UUCP (04/11/87)
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 6 (of 6)."
# Contents: tests/c2c.c tests/cc2c.c
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:49 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f tests/c2c.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/c2c.c\"
else
echo shar: Extracting \"tests/c2c.c\" \(11638 characters\)
sed "s/^X//" >tests/c2c.c <<'END_OF_tests/c2c.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * c2c.c test complex to complex math functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect a single
X * double precision complex argument and return a double
X * precision complex result.
X *
X * Most of the test data in the current data file (c2c.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first. Naturally, an error in a lower level
X * function will cause propagation of errors up to higher
X * level functions.
X *
X * USAGE
X *
X * c2c [-esv] [-l limit]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -l => report errors greater than
X * specified limit (default 10**-6)
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (c2c.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit.
X * Each test directive has the form:
X *
X * <name> <argument> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * Both the argument and the expected result are given
X * in the form:
X *
X * <real> <imag>
X *
X * where <real> is the real part and <imag> is the
X * imaginary part, separated by one or more spaces.
X *
X *
X * NOTE
X *
X * This program uses both the csubt and the cdiv routines,
X * which are pml functions!. Thus if either of these screw
X * up, the results will be unpredictable. BEWARE!
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X
X#ifndef NODBG
X#include <dbg.h>
X#else
X#include "dbg.h"
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xstatic double max_abs_err = MAX_ABS_ERR;
X
X/*
X * External functions which are used internally.
X */
X
Xextern char *strtok ();
Xextern double atof ();
Xextern double cabs ();
Xextern COMPLEX csubt ();
Xextern COMPLEX cdiv ();
X
X/*
X * External functions to be tested.
X */
X
Xextern COMPLEX cacos();
Xextern COMPLEX casin();
Xextern COMPLEX catan();
Xextern COMPLEX ccos();
Xextern COMPLEX ccosh();
Xextern COMPLEX cexp();
Xextern COMPLEX clog();
Xextern COMPLEX crcp();
Xextern COMPLEX csin();
Xextern COMPLEX csinh();
Xextern COMPLEX csqrt();
Xextern COMPLEX ctan();
Xextern COMPLEX ctanh();
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X COMPLEX (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "cacos", cacos, 0.0, /* Complex arc cosine */
X "casin", casin, 0.0, /* Complex arc sine */
X "catan", catan, 0.0, /* Complex arc tangent */
X "ccos", ccos, 0.0, /* Complex cosine */
X "ccosh", ccosh, 0.0, /* Complex hyperbolic cosine */
X "cexp", cexp, 0.0, /* Complex exponential */
X "clog", clog, 0.0, /* Complex natural logarithm */
X "crcp", crcp, 0.0, /* Complex reciprocal */
X "csin", csin, 0.0, /* Complex sine */
X "csinh", csinh, 0.0, /* Complex hyperbolic sine */
X "csqrt", csqrt, 0.0, /* Complex square root */
X "ctan", ctan, 0.0, /* Complex tangent */
X "ctanh", ctanh, 0.0, /* Complex hyperbolic tangent */
X NULL, NULL, 0.0 /* Function list end marker */
X};
X
X
X/*
X * FUNCTION
X *
X * main entry point for c2c test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain (argc, argv)
Xint argc;
Xchar *argv[];
X{
X ENTER ("main");
X DEBUGWHO (argv[0]);
X options (argc, argv);
X dotests (argv);
X statistics ();
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default real part of argument to 0.0
X * Default imaginary part of argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
X
Xdotests (argv)
Xchar *argv[];
X{
X char buffer[256]; /* Directive buffer */
X char function[64]; /* Specified function name */
X COMPLEX argument; /* Specified function argument */
X COMPLEX expected; /* Specified expected result */
X COMPLEX result; /* Actual result */
X COMPLEX error; /* Relative or absolute error */
X double abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup (); /* Returns function test pointer */
X register char *strp; /* Pointer to next token in string */
X
X ENTER ("dotests");
X while (fgets (buffer, sizeof(buffer), stdin) != NULL) {
X strcpy (function, "{null}");
X argument.real = argument.imag = 0.0;
X expected.real = expected.imag = 0.0;
X sscanf (buffer, "%s %le %le %le %le", function,
X &argument.real, &argument.imag, &expected.real, &expected.imag);
X testp = lookup (function);
X if (testp == NULL) {
X fprintf (stderr, "%s: unknown function \"%s\".\n",
X argv[0], function);
X } else {
X result = (*testp -> func)(argument);
X if (vflag) {
X printf ("%s(%le + j %le) \n = %30.23le + j %30.23le.\n",
X function, argument.real, argument.imag, result.real,
X result.imag);
X }
X if (expected.real != 0.0 || expected.imag != 0.0) {
X error = csubt (result, expected);
X error = cdiv (error, expected);
X } else {
X error = result;
X }
X abs_err = cabs (error);
X if ((abs_err > max_abs_err) || eflag) {
X fprintf (stderr, "%s: error in \"%s\"\n", argv[0], function);
X fprintf (stderr, "\treal (arg)\t\t%25.20le\n", argument.real);
X fprintf (stderr, "\timag (arg)\t\t%25.20le\n", argument.imag);
X fprintf (stderr, "\treal (result)\t\t%25.20le\n",result.real);
X fprintf (stderr, "\timag (result)\t\t%25.20le\n",result.imag);
X fprintf (stderr, "\treal (expected)\t\t%25.20le\n",expected.real);
X fprintf (stderr, "\timag (expected)\t\t%25.20le\n",expected.imag);
X }
X if (abs_err > testp -> max_err) {
X testp -> max_err = abs_err;
X }
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
X
Xoptions (argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int flag;
X extern int getopt ();
X extern char *optarg;
X
X ENTER ("options");
X eflag = sflag = vflag = FALSE;
X while ((flag = getopt (argc, argv, "#:el:sv")) != EOF) {
X switch (flag) {
X case '#':
X DEBUGPUSH (optarg);
X break;
X case 'e':
X eflag = TRUE;
X break;
X case 'l':
X sscanf (optarg, "%le", &max_abs_err);
X DEBUG3 ("args", "max_abs_err = %le", max_abs_err);
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup (funcname)
Xchar *funcname;
X{
X struct test *testp;
X struct test *rtnval;
X
X ENTER ("lookup");
X rtnval = (struct test *) NULL;
X for (testp = tests; testp -> name != NULL && rtnval == NULL; testp++) {
X if (!strcmp (testp -> name, funcname)) {
X rtnval = testp;
X }
X }
X LEAVE ();
X return (rtnval);
X}
X
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics ()
X{
X struct test *tp;
X
X ENTER ("statistics");
X if (sflag) {
X for (tp = tests; tp -> name != NULL; tp++) {
X printf ("%s:\tmaximum relative error %le\n",
X tp -> name, tp -> max_err);
X }
X }
X LEAVE ();
X}
END_OF_tests/c2c.c
if test 11638 -ne `wc -c <tests/c2c.c`; then
echo shar: \"tests/c2c.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/cc2c.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/cc2c.c\"
else
echo shar: Extracting \"tests/cc2c.c\" \(11225 characters\)
sed "s/^X//" >tests/cc2c.c <<'END_OF_tests/cc2c.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FILE
X *
X * cc2c.c test complex to complex math functions
X *
X * KEY WORDS
X *
X * portable math library
X * test functions
X *
X * DESCRIPTION
X *
X * Tests double precision functions for the Portable Math
X * Library. Tests those functions which expect two
X * double precision complex arguments and return a double
X * precision complex result.
X *
X * Most of the test data in the current data file (cc2c.dat)
X * was generated using double precision FORTRAN arithmetic
X * on a Decsystem-20.
X *
X * Note that the ordering of functions is important for
X * optimum error information. Since some functions call
X * others in the library, the functions being called should
X * be tested first. Naturally, an error in a lower level
X * function will cause propagation of errors up to higher
X * level functions.
X *
X * USAGE
X *
X * cc2c [-esv] [-l limit]
X *
X * -e => force error for each test
X * to verify error handling
X *
X * -l => report errors greater than
X * specified limit (default 10**-6)
X *
X * -s => print summary after tests
X *
X * -v => print each function, argument,
X * and result
X *
X * Test directives are read from the standard input, which
X * may be redirected to the provided test file (cc2c.dat),
X * and any relative errors are automatically written to standard
X * output if they exceed a maximum allowable limit.
X * Each test directive has the form:
X *
X * <name> <arg1> <arg2> <expected result>
X *
X * Each field is separated by a one or more space character(s).
X * Both the argument and the expected result are given
X * in the form:
X *
X * <real> <imag>
X *
X * where <real> is the real part and <imag> is the
X * imaginary part, separated by one or more spaces.
X *
X *
X * NOTE
X *
X * This program uses both the csubt and the cdiv routines,
X * which are pml functions!. Thus if either of these screw
X * up, the results will be unpredictable. BEWARE!
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X
X#ifndef NODBG
X#include <dbg.h>
X#else
X#include "dbg.h"
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define MAX_ABS_ERR 1.0e-6 /* Set to catch only gross errors */
X
Xstatic int vflag; /* Flag for verbose option */
Xstatic int eflag; /* Simulate an error to error printout */
Xstatic int sflag; /* Flag to show final statistics */
X
Xstatic double max_abs_err = MAX_ABS_ERR;
X
X/*
X * External functions which are used internally.
X */
X
Xextern char *strtok ();
Xextern double atof ();
Xextern double cabs ();
Xextern COMPLEX csubt ();
Xextern COMPLEX cdiv ();
X
X/*
X * External functions to be tested.
X */
X
Xextern COMPLEX cadd();
Xextern COMPLEX csubt();
Xextern COMPLEX cdiv();
Xextern COMPLEX cmult();
X
X
X/*
X * Define all recognized test functions. Each function
X * must have an entry in this table, where each
X * entry contains the information specified in the
X * structure "test".
X *
X */
X
Xstruct test { /* Structure of each function to be tested */
X char *name; /* Name of the function to test */
X COMPLEX (*func)(); /* Pointer to the function's entry point */
X double max_err; /* Error accumulator for this function */
X};
X
Xstatic struct test tests[] = { /* Table of all recognized functions */
X "cadd", cadd, 0.0, /* Complex addition */
X "csubt", csubt, 0.0, /* Complex subtraction */
X "cdiv", cdiv, 0.0, /* Complex division */
X "cmult", cmult, 0.0, /* Complex multiplication */
X NULL, NULL, 0.0 /* Function list end marker */
X};
X
X
X/*
X * FUNCTION
X *
X * main entry point for cc2c test utility
X *
X * PSEUDO CODE
X *
X * Begin main
X * Process any options in command line.
X * Do all tests requested by stdin directives.
X * Report final statistics (if enabled).
X * End main
X *
X */
X
Xmain (argc, argv)
Xint argc;
Xchar *argv[];
X{
X ENTER ("main");
X DEBUGWHO (argv[0]);
X options (argc, argv);
X dotests (argv);
X statistics ();
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * dotests process each test from stdin directives
X *
X * ERROR REPORTING
X *
X * Note that in most cases, the error criterion is based
X * on relative error, defined as:
X *
X * error = (result - expected) / expected
X *
X * Naturally, if the expected result is zero, some
X * other criterion must be used. In this case, the
X * absolute error is used. That is:
X *
X * error = result
X *
X * PSEUDO CODE
X *
X * Begin dotests
X * While a test directive is successfully read from stdin
X * Default function name to "{null}"
X * Default real part of argument to 0.0
X * Default imaginary part of argument to 0.0
X * Default expected result to 0.0
X * Extract function name, argument and expected result
X * Lookup test in available test list
X * If no test was found then
X * Tell user that unknown function was specified
X * Else
X * Call function with argument and save result
X * If the verify flag is set then
X * Print function name, argument, and result
X * End if
X * If the expected result is not zero then
X * Compute the relative error
X * Else
X * Use the absolute error
X * End if
X * Get absolute value of error
X * If error exceeds limit or error force flag set
X * Print error notification on stderr
X * End if
X * If this error is max for given function then
X * Remember this error for summary
X * End if
X * End if
X * End while
X * End dotests
X *
X */
X
X
Xdotests (argv)
Xchar *argv[];
X{
X char buffer[256]; /* Directive buffer */
X char function[64]; /* Specified function name */
X COMPLEX arg1; /* Specified function argument 1 */
X COMPLEX arg2; /* Specified function argument 2 */
X COMPLEX expected; /* Specified expected result */
X COMPLEX result; /* Actual result */
X COMPLEX error; /* Relative or absolute error */
X double abs_err; /* Absolute value of error */
X struct test *testp; /* Pointer to function test */
X struct test *lookup (); /* Returns function test pointer */
X register char *strp; /* Pointer to next token in string */
X
X ENTER ("dotests");
X while (fgets (buffer, sizeof(buffer), stdin) != NULL) {
X strcpy (function, "{null}");
X arg1.real = arg1.imag = 0.0;
X arg2.real = arg2.imag = 0.0;
X expected.real = expected.imag = 0.0;
X sscanf (buffer, "%s %le %le %le %le %le %le", function,
X &arg1.real, &arg1.imag, &arg2.real, &arg2.imag,
X &expected.real, &expected.imag);
X testp = lookup (function);
X if (testp == NULL) {
X fprintf (stderr, "%s: unknown function \"%s\".\n",
X argv[0], function);
X } else {
X result = (*testp -> func)(arg1, arg2);
X if (vflag) {
X printf ("%s(%le + j %le, %le + j %le)\n",
X function, arg1.real, arg1.imag, arg2.real, arg2.imag);
X printf (" = %30.23le + j %30.23le.\n", result.real,
X result.imag);
X }
X if (expected.real != 0.0 || expected.imag != 0.0) {
X error = csubt (result, expected);
X error = cdiv (error, expected);
X } else {
X error = result;
X }
X abs_err = cabs (error);
X if ((abs_err > max_abs_err) || eflag) {
X fprintf (stderr, "%s: error in \"%s\"\n", argv[0], function);
X fprintf (stderr, "\treal (arg1)\t\t%25.20le\n", arg1.real);
X fprintf (stderr, "\timag (arg1)\t\t%25.20le\n", arg1.imag);
X fprintf (stderr, "\treal (arg2)\t\t%25.20le\n", arg2.real);
X fprintf (stderr, "\timag (arg2)\t\t%25.20le\n", arg2.imag);
X fprintf (stderr, "\treal (result)\t\t%25.20le\n",result.real);
X fprintf (stderr, "\timag (result)\t\t%25.20le\n",result.imag);
X fprintf (stderr, "\treal (expected)\t\t%25.20le\n",expected.real);
X fprintf (stderr, "\timag (expected)\t\t%25.20le\n",expected.imag);
X }
X if (abs_err > testp -> max_err) {
X testp -> max_err = abs_err;
X }
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * options process command line options
X *
X * PSEUDO CODE
X *
X * Begin options
X * Reset all flags to FALSE by default
X * Initialize flag argument scan pointer
X * If there is a second command line argument then
X * If the argument specifies flags then
X * While there is an unprocessed flag
X * Switch on flag
X * Case "force error flag":
X * Set the "force error" flag
X * Break switch
X * Case "print summary":
X * Set "print summary" flag
X * Break switch
X * Case "verbose":
X * Set "verbose" flag
X * Break switch
X * Default:
X * Tell user unknown flag
X * Break switch
X * End switch
X * End while
X * End if
X * End if
X * End options
X *
X */
X
X
Xoptions (argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int flag;
X extern int getopt ();
X extern char *optarg;
X
X ENTER ("options");
X eflag = sflag = vflag = FALSE;
X while ((flag = getopt (argc, argv, "#:el:sv")) != EOF) {
X switch (flag) {
X case '#':
X DEBUGPUSH (optarg);
X break;
X case 'e':
X eflag = TRUE;
X break;
X case 'l':
X sscanf (optarg, "%le", &max_abs_err);
X DEBUG3 ("args", "max_abs_err = %le", max_abs_err);
X break;
X case 's':
X sflag = TRUE;
X break;
X case 'v':
X vflag = TRUE;
X break;
X }
X }
X LEAVE ();
X}
X
X
X/*
X * FUNCTION
X *
X * loopup lookup test in known test list
X *
X * DESCRIPTION
X *
X * Given the name of a desired test, looks up the test
X * in the known test list and returns a pointer to the
X * test structure.
X *
X * Since the table is so small we simply use a linear
X * search.
X *
X * PSEUDO CODE
X *
X * Begin lookup
X * For each known test
X * If the test's name matches the desired test name
X * Return pointer to the test structure
X * End if
X * End for
X * End lookup
X *
X */
X
Xstruct test *lookup (funcname)
Xchar *funcname;
X{
X struct test *testp;
X struct test *rtnval;
X
X ENTER ("lookup");
X rtnval = (struct test *) NULL;
X for (testp = tests; testp -> name != NULL && rtnval == NULL; testp++) {
X if (!strcmp (testp -> name, funcname)) {
X rtnval = testp;
X }
X }
X LEAVE ();
X return (rtnval);
X}
X
X
X/*
X * FUNCTION
X *
X * statistics print final statistics if desired
X *
X * PSEUDO CODE
X *
X * Begin statistics
X * If a final statistics (summary) is desired then
X * For each test in the known test list
X * Print the maximum error encountered
X * End for
X * End if
X * End statistics
X *
X */
X
Xstatistics ()
X{
X struct test *tp;
X
X ENTER ("statistics");
X if (sflag) {
X for (tp = tests; tp -> name != NULL; tp++) {
X printf ("%s:\tmaximum relative error %le\n",
X tp -> name, tp -> max_err);
X }
X }
X LEAVE ();
X}
END_OF_tests/cc2c.c
if test 11225 -ne `wc -c <tests/cc2c.c`; then
echo shar: \"tests/cc2c.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 6 \(of 6\).
cp /dev/null ark6isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =