Portable Math Library in C (Part 2 of 6)

[f...@mcdsun.uucp]

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 =