bammi@dsrgsun.ces.cwru.edu (Jwahar R. Bammi) (12/14/88)
#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# README.GCC
# README
# pmlsrc/*
# This archive created: Wed Dec 14 04:41:08 1988
# By: Jwahar R. Bammi(Case Western Reserve University)
# Uucp: {decvax,sun,att}!cwjcc!dsrgsun!bammi
# Csnet: bammi@dsrgsun.ces.CWRU.edu
# Arpa: bammi@dsrgsun.ces.CWRU.edu
#
export PATH; PATH=/bin:$PATH
echo shar: extracting "'README.GCC'" '(1628 characters)'
if test -f 'README.GCC'
then
echo shar: over-writing existing file "'README.GCC'"
fi
sed 's/^X//' << \SHAR_EOF > 'README.GCC'
XThis is a upgrade package for Gcc V1.31 for both MinixST gcc
Xand Tos gcc V1.31. This upgrade does not make any changes to the compiler
Xitself, but
X
X - fixes bugs in the libraries (both tos and minix gcc libraries)
X - adds a few functions to the libraries
X - adds a math library based on Fred Fish's Portable Math Library
X package that appeared on comp.sources.misc.
X
XBefore you begin, please check that you have gcc V1.31. None of this
Xwill work with any prior version of gcc (this warning applies to
Xboth gcc-tos and gcc-minixst).
X
XTo install:
X -Gcc-tos users:
X Upgrade your gcc library by
X - applying the patches in tos.cdiff.
X - copying the files ldexp.c, modf.c, frexp.c and dflonum.c
X into the library source directory.
X - compile and install the gcc library again.
X - cd to pmlsrc, and compile and install libpml.a
X using Makefile.tos.
X - cd to pmltests and make tests using Makefile.tos
X - run tests
X
X -Gcc-minix users:
X Upgrade your gcc library by
X - applying the patches in minix.cdiff.
X - copying the files ldexp.c, modf.c, frexp.c and dflonum.c
X into the library source directory.
X - compile and install the gcc library again.
X - cd to pmlsrc, and compile and install libpml.a and libpml32.a
X using Makefile.minix
X - cd to pmltests and make tests using Makefile.minix
X - run tests
X
X to run a test enter:
X TEST -v -s < TEST.dat
X where TEST is one of c2c cc2c c2d d2d or dd2d
X
XPlease direct you comments and suggestions to:
X--
Xusenet: {decvax,sun}!cwjcc!dsrgsun!bammi jwahar r. bammi
Xcsnet: bammi@dsrgsun.ces.CWRU.edu
Xarpa: bammi@dsrgsun.ces.CWRU.edu
XcompuServe: 71515,155
SHAR_EOF
if test 1628 -ne "`wc -c 'README.GCC'`"
then
echo shar: error transmitting "'README.GCC'" '(should have been 1628 characters)'
fi
echo shar: extracting "'README'" '(3902 characters)'
if test -f 'README'
then
echo shar: over-writing existing file "'README'"
fi
sed 's/^X//' << \SHAR_EOF > '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
SHAR_EOF
if test 3902 -ne "`wc -c 'README'`"
then
echo shar: error transmitting "'README'" '(should have been 3902 characters)'
fi
if test ! -d 'pmlsrc'
then
echo shar: creating directory "'pmlsrc'"
mkdir 'pmlsrc'
fi
echo shar: entering directory "'pmlsrc'"
cd 'pmlsrc'
echo shar: extracting "'acos.c'" '(2962 characters)'
if test -f 'acos.c'
then
echo shar: over-writing existing file "'acos.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
SHAR_EOF
if test 2962 -ne "`wc -c 'acos.c'`"
then
echo shar: error transmitting "'acos.c'" '(should have been 2962 characters)'
fi
echo shar: extracting "'acosh.c'" '(2642 characters)'
if test -f 'acosh.c'
then
echo shar: over-writing existing file "'acosh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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
SHAR_EOF
if test 2642 -ne "`wc -c 'acosh.c'`"
then
echo shar: error transmitting "'acosh.c'" '(should have been 2642 characters)'
fi
echo shar: extracting "'asin.c'" '(2675 characters)'
if test -f 'asin.c'
then
echo shar: over-writing existing file "'asin.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2675 -ne "`wc -c 'asin.c'`"
then
echo shar: error transmitting "'asin.c'" '(should have been 2675 characters)'
fi
echo shar: extracting "'asinh.c'" '(2315 characters)'
if test -f 'asinh.c'
then
echo shar: over-writing existing file "'asinh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2315 -ne "`wc -c 'asinh.c'`"
then
echo shar: error transmitting "'asinh.c'" '(should have been 2315 characters)'
fi
echo shar: extracting "'atan.c'" '(5187 characters)'
if test -f 'atan.c'
then
echo shar: over-writing existing file "'atan.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 5187 -ne "`wc -c 'atan.c'`"
then
echo shar: error transmitting "'atan.c'" '(should have been 5187 characters)'
fi
echo shar: extracting "'atan2.c'" '(2363 characters)'
if test -f 'atan2.c'
then
echo shar: over-writing existing file "'atan2.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2363 -ne "`wc -c 'atan2.c'`"
then
echo shar: error transmitting "'atan2.c'" '(should have been 2363 characters)'
fi
echo shar: extracting "'atanh.c'" '(2301 characters)'
if test -f 'atanh.c'
then
echo shar: over-writing existing file "'atanh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2301 -ne "`wc -c 'atanh.c'`"
then
echo shar: error transmitting "'atanh.c'" '(should have been 2301 characters)'
fi
echo shar: extracting "'cabs.c'" '(1722 characters)'
if test -f 'cabs.c'
then
echo shar: over-writing existing file "'cabs.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'cabs.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 * cabs double precision complex absolute value
X *
X * KEY WORDS
X *
X * cabs
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision absolute value of a double
X * precision complex argument, where "absolute value"
X * is taken to mean magnitude. The result replaces the
X * argument.
X *
X * USAGE
X *
X * double cabs (z)
X * COMPLEX z;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az
X *
X * INTERNALS
X *
X * Computes cabs (z) where z = (x) + j(y) from:
X *
X * cabs (z) = sqrt (x*x + y*y)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble cabs (z)
XCOMPLEX z;
X{
X double result;
X extern double sqrt ();
X
X ENTER ("cabs");
X DEBUG4 ("cabsin", "arg %le +j %le", z.real, z.imag);
X result = sqrt ((z.real * z.real) + (z.imag * z.imag));
X DEBUG3 ("cabsout", "result %le", result);
X LEAVE ();
X return (result);
X}
SHAR_EOF
if test 1722 -ne "`wc -c 'cabs.c'`"
then
echo shar: error transmitting "'cabs.c'" '(should have been 1722 characters)'
fi
echo shar: extracting "'cacos.c'" '(1983 characters)'
if test -f 'cacos.c'
then
echo shar: over-writing existing file "'cacos.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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
SHAR_EOF
if test 1983 -ne "`wc -c 'cacos.c'`"
then
echo shar: error transmitting "'cacos.c'" '(should have been 1983 characters)'
fi
echo shar: extracting "'cadd.c'" '(1820 characters)'
if test -f 'cadd.c'
then
echo shar: over-writing existing file "'cadd.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'cadd.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 * cadd double precision complex addition
X *
X * KEY WORDS
X *
X * cadd
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex result of addition of
X * first double precision complex argument with second double
X * precision complex argument.
X *
X * Note that the complex addition function is
X * so simple that it would not normally be called as a function
X * but simply done "inline". It is supplied mostly for
X * completeness.
X *
X * USAGE
X *
X * COMPLEX cadd (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 cadd(z1,z2) from:
X *
X * 1. Let z1 = a + j b
X * Let z2 = c + j d
X *
X * 2. Then cadd(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 cadd (z1, z2)
XCOMPLEX z1;
XCOMPLEX z2;
X{
X ENTER ("cadd");
X z1.real += z2.real;
X z1.imag += z2.imag;
X LEAVE ();
X return (z1);
X}
SHAR_EOF
if test 1820 -ne "`wc -c 'cadd.c'`"
then
echo shar: error transmitting "'cadd.c'" '(should have been 1820 characters)'
fi
echo shar: extracting "'casin.c'" '(1876 characters)'
if test -f 'casin.c'
then
echo shar: over-writing existing file "'casin.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'casin.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 * casin complex double precision arc sine
X *
X * KEY WORDS
X *
X * casin
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex arc sine of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX casin (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 sine of z = x + j y from:
X *
X * casin(z) = -j * clog(csqrt(1-z*z) + j*z)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX casin (z)
XCOMPLEX z;
X{
X COMPLEX temp;
X extern COMPLEX csqrt (), clog (), cmult ();
X
X ENTER ("casin");
X DEBUG4 ("casinin", "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 temp.real -= z.imag;
X temp.imag += z.real;
X temp = clog (temp);
X z.real = temp.imag;
X z.imag = -temp.real;
X DEBUG4 ("casinout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1876 -ne "`wc -c 'casin.c'`"
then
echo shar: error transmitting "'casin.c'" '(should have been 1876 characters)'
fi
echo shar: extracting "'catan.c'" '(1857 characters)'
if test -f 'catan.c'
then
echo shar: over-writing existing file "'catan.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'catan.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 * catan complex double precision arc tangent
X *
X * KEY WORDS
X *
X * catan
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex arc tangent of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX catan (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 tangent of z = x + j y from:
X *
X * catan(z) = -j/2 * clog( (j+z) / (j-z) )
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX catan (z)
XCOMPLEX z;
X{
X COMPLEX temp;
X double swaptemp;
X extern COMPLEX cdiv (), clog ();
X
X ENTER ("catan");
X DEBUG4 ("catanin", "arg %le %le", z.real, z.imag);
X temp.real = -z.real;
X temp.imag = 1.0 - z.imag;
X z.imag += 1.0;
X z = cdiv (z, temp);
X z = clog (z);
X swaptemp = z.real;
X z.real = -0.5 * z.imag;
X z.imag = 0.5 * swaptemp;
X DEBUG4 ("catanout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1857 -ne "`wc -c 'catan.c'`"
then
echo shar: error transmitting "'catan.c'" '(should have been 1857 characters)'
fi
echo shar: extracting "'ccos.c'" '(1881 characters)'
if test -f 'ccos.c'
then
echo shar: over-writing existing file "'ccos.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'ccos.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 * ccos complex double precision cosine
X *
X * KEY WORDS
X *
X * ccos
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex cosine of a double
X * precision complex argument.
X *
X * USAGE
X *
X * COMPLEX ccos (z)
X * COMPLEX z;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-12
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex cosine of z = x + j y from:
X *
X * 1. r_ccos = cos(x) * cosh(y)
X *
X * 2. i_ccos = -sin(x) * sinh(y)
X *
X * 3. ccos(z) = r_ccos + j i_ccos
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX ccos (z)
XCOMPLEX z;
X{
X COMPLEX result;
X extern double sin(), cos(), sinh(), cosh();
X
X ENTER ("ccos");
X DEBUG4 ("ccosin", "arg %le %le", z.real, z.imag);
X result.real = cos(z.real) * cosh(z.imag);
X result.imag = -sin(z.real) * sinh(z.imag);
X DEBUG4 ("ccosout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
SHAR_EOF
if test 1881 -ne "`wc -c 'ccos.c'`"
then
echo shar: error transmitting "'ccos.c'" '(should have been 1881 characters)'
fi
echo shar: extracting "'ccosh.c'" '(1842 characters)'
if test -f 'ccosh.c'
then
echo shar: over-writing existing file "'ccosh.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'ccosh.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 * ccosh complex double precision hyperbolic cosine
X *
X * KEY WORDS
X *
X * ccosh
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex hyperbolic cosine of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX ccosh (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 hyperbolic cosine of Z = x + j y from:
X *
X * ccosh(z) = 0.5 * ( cexp(z) + cexp(-z) )
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX ccosh (z)
XCOMPLEX z;
X{
X COMPLEX cexpmz;
X extern COMPLEX cexp ();
X
X ENTER ("ccosh");
X DEBUG4 ("ccoshin", "arg %le %le", z.real, z.imag);
X cexpmz.real = -z.real;
X cexpmz.imag = -z.imag;
X cexpmz = cexp (cexpmz);
X z = cexp (z);
X z.real += cexpmz.real;
X z.imag += cexpmz.imag;
X z.real *= 0.5;
X z.imag *= 0.5;
X DEBUG4 ("ccoshout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1842 -ne "`wc -c 'ccosh.c'`"
then
echo shar: error transmitting "'ccosh.c'" '(should have been 1842 characters)'
fi
echo shar: extracting "'cdiv.c'" '(2460 characters)'
if test -f 'cdiv.c'
then
echo shar: over-writing existing file "'cdiv.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2460 -ne "`wc -c 'cdiv.c'`"
then
echo shar: error transmitting "'cdiv.c'" '(should have been 2460 characters)'
fi
echo shar: extracting "'cexp.c'" '(1910 characters)'
if test -f 'cexp.c'
then
echo shar: over-writing existing file "'cexp.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'cexp.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 * cexp complex double precision exponential
X *
X * KEY WORDS
X *
X * cexp
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex exponential of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX cexp (z)
X * COMPLEX z;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-13
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex exponential of z = x + j y from:
X *
X * 1. r_cexp = exp(x) * cos(y)
X *
X * 2. i_cexp = exp(x) * sin(y)
X *
X * 3. cexp(z) = r_cexp + j i_cexp
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX cexp (z)
XCOMPLEX z;
X{
X COMPLEX result;
X double expx;
X extern double exp(), sin(), cos();
X
X ENTER ("cexp");
X DEBUG4 ("cexpin", "arg %le %le", z.real, z.imag);
X expx = exp (z.real);
X result.real = expx * cos (z.imag);
X result.imag = expx * sin (z.imag);
X DEBUG4 ("cexpout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
SHAR_EOF
if test 1910 -ne "`wc -c 'cexp.c'`"
then
echo shar: error transmitting "'cexp.c'" '(should have been 1910 characters)'
fi
echo shar: extracting "'clog.c'" '(1870 characters)'
if test -f 'clog.c'
then
echo shar: over-writing existing file "'clog.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'clog.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 * clog complex double precision natural logarithm
X *
X * KEY WORDS
X *
X * clog
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex natural logarithm of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX clog (z)
X * COMPLEX z;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-13
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex natural logarithm of z = x + j y from:
X *
X * 1. r_clog = log(cabs(z))
X *
X * 2. i_clog = atan2(x,y)
X *
X * 3. clog(z) = r_clog + j i_clog
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX clog (z)
XCOMPLEX z;
X{
X double temp;
X extern double cabs (), atan2(), log ();
X
X ENTER ("clog");
X DEBUG4 ("clogin", "arg %le %le", z.real, z.imag);
X temp = log (cabs (z));
X z.imag = atan2 (z.real, z.imag);
X z.real = temp;
X DEBUG4 ("clogout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1870 -ne "`wc -c 'clog.c'`"
then
echo shar: error transmitting "'clog.c'" '(should have been 1870 characters)'
fi
echo shar: extracting "'cmult.c'" '(1994 characters)'
if test -f 'cmult.c'
then
echo shar: over-writing existing file "'cmult.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 1994 -ne "`wc -c 'cmult.c'`"
then
echo shar: error transmitting "'cmult.c'" '(should have been 1994 characters)'
fi
echo shar: extracting "'cos.c'" '(5011 characters)'
if test -f 'cos.c'
then
echo shar: over-writing existing file "'cos.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 5011 -ne "`wc -c 'cos.c'`"
then
echo shar: error transmitting "'cos.c'" '(should have been 5011 characters)'
fi
echo shar: extracting "'cosh.c'" '(2424 characters)'
if test -f 'cosh.c'
then
echo shar: over-writing existing file "'cosh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2424 -ne "`wc -c 'cosh.c'`"
then
echo shar: error transmitting "'cosh.c'" '(should have been 2424 characters)'
fi
echo shar: extracting "'crcp.c'" '(1919 characters)'
if test -f 'crcp.c'
then
echo shar: over-writing existing file "'crcp.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'crcp.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 * crcp complex double precision reciprocal
X *
X * KEY WORDS
X *
X * crcp
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex reciprocal of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX crcp (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 reciprocal of z = x + j y from:
X *
X * 1. Compute denom = x*x + y*y
X *
X * 2. If denom = 0.0 then flag error
X * and return MAX_POS_DBLF + j 0.0
X *
X * 3. Else crcp(z) = (x/denom) + j (-y/denom)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX crcp (z)
XCOMPLEX z;
X{
X double denom;
X
X ENTER ("crcp");
X DEBUG4 ("crcpin", "arg %le %le", z.real, z.imag);
X denom = (z.real * z.real) + (z.imag * z.imag);
X if (denom == 0.0) {
X/*****
X pmlerr(CRCP_OF_ZERO);
X z.real = MAX_POS_DBLF;
X******/
X z.real = 0.0; /* TERRIBLY WRONG */
X z.imag = 0.0;
X } else {
X z.real /= denom;
X z.imag /= -denom;
X }
X DEBUG4 ("crcpout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1919 -ne "`wc -c 'crcp.c'`"
then
echo shar: error transmitting "'crcp.c'" '(should have been 1919 characters)'
fi
echo shar: extracting "'csin.c'" '(1876 characters)'
if test -f 'csin.c'
then
echo shar: over-writing existing file "'csin.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'csin.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 * csin complex double precision sine
X *
X * KEY WORDS
X *
X * csin
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex sine of a double
X * precision complex argument.
X *
X * USAGE
X *
X * COMPLEX csin (z)
X * COMPLEX z;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-12
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X * INTERNALS
X *
X * Computes complex sine of z = x + j y from:
X *
X * 1. r_csin = sin(x) * cosh(y)
X *
X * 2. i_csin = cos(x) * sinh(y)
X *
X * 3. csin(z) = r_csin + j i_csin
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX csin (z)
XCOMPLEX z;
X{
X COMPLEX result;
X extern double sin(), cos(), sinh(), cosh();
X
X ENTER ("csin");
X DEBUG4 ("csinin", "arg %le %le", z.real, z.imag);
X result.real = sin (z.real) * cosh (z.imag);
X result.imag = cos (z.real) * sinh (z.imag);
X DEBUG4 ("csinout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
SHAR_EOF
if test 1876 -ne "`wc -c 'csin.c'`"
then
echo shar: error transmitting "'csin.c'" '(should have been 1876 characters)'
fi
echo shar: extracting "'csinh.c'" '(1836 characters)'
if test -f 'csinh.c'
then
echo shar: over-writing existing file "'csinh.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'csinh.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 * csinh complex double precision hyperbolic sine
X *
X * KEY WORDS
X *
X * csinh
X * machine independent routines
X * complex functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex hyperbolic sine of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX csinh (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 hyperbolic sine of z = x + j y from:
X *
X * csinh(z) = 0.5 * ( cexp(z) - cexp(-z) )
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX csinh (z)
XCOMPLEX z;
X{
X COMPLEX cexpmz;
X extern COMPLEX cexp ();
X
X ENTER ("csinh");
X DEBUG4 ("csinhin", "arg %le %le", z.real, z.imag);
X cexpmz.real = -z.real;
X cexpmz.imag = -z.imag;
X cexpmz = cexp (cexpmz);
X z = cexp (z);
X z.real -= cexpmz.real;
X z.imag -= cexpmz.imag;
X z.real *= 0.5;
X z.imag *= 0.5;
X DEBUG4 ("csinhout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1836 -ne "`wc -c 'csinh.c'`"
then
echo shar: error transmitting "'csinh.c'" '(should have been 1836 characters)'
fi
echo shar: extracting "'csqrt.c'" '(2679 characters)'
if test -f 'csqrt.c'
then
echo shar: over-writing existing file "'csqrt.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2679 -ne "`wc -c 'csqrt.c'`"
then
echo shar: error transmitting "'csqrt.c'" '(should have been 2679 characters)'
fi
echo shar: extracting "'csubt.c'" '(2005 characters)'
if test -f 'csubt.c'
then
echo shar: over-writing existing file "'csubt.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2005 -ne "`wc -c 'csubt.c'`"
then
echo shar: error transmitting "'csubt.c'" '(should have been 2005 characters)'
fi
echo shar: extracting "'ctan.c'" '(1942 characters)'
if test -f 'ctan.c'
then
echo shar: over-writing existing file "'ctan.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'ctan.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 * ctan complex double precision tangent
X *
X * KEY WORDS
X *
X * ctan
X * complex functions
X * machine independent functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex tangent of a double
X * precision complex argument.
X *
X * USAGE
X *
X * COMPLEX ctan (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 tangent of z = x + j y from:
X *
X * 1. Compute ccos(z)
X *
X * 2. If ccos(z) = 0 + j0 then the
X * result is MAX_POS_DBLF + j0
X *
X * 3. Else ctan(z) = csin(z) / ccos(z)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX ctan (z)
XCOMPLEX z;
X{
X COMPLEX ccosz;
X extern COMPLEX ccos (), csin (), cdiv ();
X
X ENTER ("ctan");
X DEBUG4 ("ctanin", "arg %le %le", z.real, z.imag);
X ccosz = ccos (z);
X if (ccosz.real == 0.0 && ccosz.imag == 0.0) {
X/*****
X z.real = MAX_POS_DBLF;
X******/
X z.real = 0.0; /* TERRIBLY WRONG! */
X z.imag = 0.0;
X } else {
X z = csin (z);
X z = cdiv (z, ccosz);
X }
X DEBUG4 ("ctanout", "result %le %le", z.real, z.imag);
X LEAVE ();
X return (z);
X}
SHAR_EOF
if test 1942 -ne "`wc -c 'ctan.c'`"
then
echo shar: error transmitting "'ctan.c'" '(should have been 1942 characters)'
fi
echo shar: extracting "'ctanh.c'" '(1898 characters)'
if test -f 'ctanh.c'
then
echo shar: over-writing existing file "'ctanh.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'ctanh.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 * ctanh complex double precision hyperbolic tangent
X *
X * KEY WORDS
X *
X * ctanh
X * complex functions
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Computes double precision complex hyperbolic tangent of
X * a double precision complex argument.
X *
X * USAGE
X *
X * COMPLEX ctanh (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 hyperbolic tangent of Z = x + j y from:
X *
X * ctanh(z) = (1 - cexp(-2z)) / (1 + cexp(-2z))
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
XCOMPLEX ctanh (z)
XCOMPLEX z;
X{
X COMPLEX result;
X extern COMPLEX cexp (), cdiv ();
X
X ENTER ("ctanh");
X DEBUG4 ("ctanhin", "arg %le %le", z.real, z.imag);
X result.real = -2.0 * z.real;
X result.imag = -2.0 * z.imag;
X result = cexp (result);
X z.real = 1.0 - result.real;
X z.imag = -result.imag;
X result.real += 1.0;
X result = cdiv (z, result);
X DEBUG4 ("ctanhout", "result %le %le", result.real, result.imag);
X LEAVE ();
X return (result);
X}
SHAR_EOF
if test 1898 -ne "`wc -c 'ctanh.c'`"
then
echo shar: error transmitting "'ctanh.c'" '(should have been 1898 characters)'
fi
echo shar: done with directory "'pmlsrc'"
cd ..
# End of shell archive
exit 0
usenet: {decvax,sun}!cwjcc!dsrgsun!bammi jwahar r. bammi
csnet: bammi@dsrgsun.ces.CWRU.edu
arpa: bammi@dsrgsun.ces.CWRU.edu
compuServe: 71515,155