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:
# pmlsrc/* (continued)
# This archive created: Wed Dec 14 04:42:19 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
if test ! -d 'pmlsrc'
then
echo shar: creating directory "'pmlsrc'"
mkdir 'pmlsrc'
fi
echo shar: entering directory "'pmlsrc'"
cd 'pmlsrc'
echo shar: extracting "'Makefile.16'" '(1580 characters)'
if test -f 'Makefile.16'
then
echo shar: over-writing existing file "'Makefile.16'"
fi
sed 's/^X//' << \SHAR_EOF > 'Makefile.16'
X#
X# FILE
X#
X# Makefile build and install the pml library (16 bit ints)
X#
X# SYNOPSIS
X#
X# make funcs make version of library in local directory
X# make install install the library (must be root)
X#
X# WARNING
X#
X# The order of the modules listed in the "LEVEL<n>" macros
X# is significant since these are the orders in which
X# they will be loaded into the library archive. Although
X# some machines support randomly ordered libraries, ordering
X# them correctly doesn't hurt...
XCROSSDIR = /dsrg/bammi/cross-minix
XCROSSLIB = $(CROSSDIR)/lib
XCROSSBIN = $(CROSSDIR)/bin
XCROSSINC = $(CROSSDIR)/include
X
XAR = /dsrg/bammi/cross-gcc/bin/car
XCC = $(CROSSBIN)/mgcc
XCFLAGS = -mshort -O -DIEEE -DNO_DBUG -I.
X
XLIB = $(CROSSLIB)
X
XINC = $(CROSSINC)
X
X
XLEVEL0 = matherr.o
X
XLEVEL1 = sign.o mod.o poly.o dabs.o sqrt.o
X
XLEVEL2 = acos.o acosh.o asin.o asinh.o atan2.o atan.o \
X log10.o tan.o \
X tanh.o sinh.o cosh.o atanh.o \
X log.o sin.o cos.o exp.o max.o min.o
X
XLEVEL3 = casin.o cacos.o cmult.o catan.o ccosh.o clog.o \
X crcp.o csinh.o csqrt.o ctan.o ctanh.o cexp.o \
X ccos.o csin.o cdiv.o csubt.o cabs.o
X
XLEVEL4 = cadd.o
X
XOBJ = $(LEVEL4) $(LEVEL3) $(LEVEL2) $(LEVEL1) $(LEVEL0)
X
X#
X# The default thing to make.
X#
X
Xdefault: libpml.a
X
Xlibpml.a: $(OBJ)
X rm -f libpml.a
X $(AR) srv libpml.a $(OBJ)
X
X$(OBJ) : pml.h
X
X#
X# Stuff to do installation
X#
X
Xinstall : $(LIB)/libpml.a $(INC)/pmluser.h
X
X$(LIB)/libpml.a: libpml.a
X cp libpml.a $(LIB)/libpml.a
X
X$(INC)/pmluser.h : pmluser.h
X cp pmluser.h $(INC)/pmluser.h
X
X#
X# Clean up the directory.
X#
X
Xclean:
X rm -f *.o *.BAK libpml.a *.tmp *.bak nohup.out
SHAR_EOF
if test 1580 -ne "`wc -c 'Makefile.16'`"
then
echo shar: error transmitting "'Makefile.16'" '(should have been 1580 characters)'
fi
echo shar: extracting "'Makefile.32'" '(1592 characters)'
if test -f 'Makefile.32'
then
echo shar: over-writing existing file "'Makefile.32'"
fi
sed 's/^X//' << \SHAR_EOF > 'Makefile.32'
X#
X# FILE
X#
X# Makefile build and install the pml library (32 bit ints)
X#
X# SYNOPSIS
X#
X# make funcs make version of library in local directory
X# make install install the library (must be root)
X#
X# WARNING
X#
X# The order of the modules listed in the "LEVEL<n>" macros
X# is significant since these are the orders in which
X# they will be loaded into the library archive. Although
X# some machines support randomly ordered libraries, ordering
X# them correctly doesn't hurt...
XCROSSDIR = /dsrg/bammi/cross-minix
XCROSSLIB = $(CROSSDIR)/lib
XCROSSBIN = $(CROSSDIR)/bin
XCROSSINC = $(CROSSDIR)/include
X
XAR = /dsrg/bammi/cross-gcc/bin/car
XCC = $(CROSSBIN)/mgcc
XCFLAGS = -O -DIEEE -DNO_DBUG -I.
X
XLIB = $(CROSSLIB)
X
XINC = $(CROSSINC)
X
X
XLEVEL0 = matherr.o
X
XLEVEL1 = sign.o mod.o poly.o dabs.o sqrt.o
X
XLEVEL2 = acos.o acosh.o asin.o asinh.o atan2.o atan.o \
X log10.o tan.o \
X tanh.o sinh.o cosh.o atanh.o \
X log.o sin.o cos.o exp.o max.o min.o
X
XLEVEL3 = casin.o cacos.o cmult.o catan.o ccosh.o clog.o \
X crcp.o csinh.o csqrt.o ctan.o ctanh.o cexp.o \
X ccos.o csin.o cdiv.o csubt.o cabs.o
X
XLEVEL4 = cadd.o
X
XOBJ = $(LEVEL4) $(LEVEL3) $(LEVEL2) $(LEVEL1) $(LEVEL0)
X
X#
X# The default thing to make.
X#
X
Xdefault: libpml32.a
X
Xlibpml32.a: $(OBJ)
X rm -f libpml32.a
X $(AR) srv libpml32.a $(OBJ)
X
X$(OBJ) : pml.h
X
X#
X# Stuff to do installation
X#
X
Xinstall : $(LIB)/libpml32.a $(INC)/pmluser.h
X
X$(LIB)/libpml32.a: libpml32.a
X cp libpml32.a $(LIB)/libpml32.a
X
X$(INC)/pmluser.h : pmluser.h
X cp pmluser.h $(INC)/pmluser.h
X
X#
X# Clean up the directory.
X#
X
Xclean:
X rm -f *.o *.BAK libpml32.a *.tmp *.bak nohup.out
SHAR_EOF
if test 1592 -ne "`wc -c 'Makefile.32'`"
then
echo shar: error transmitting "'Makefile.32'" '(should have been 1592 characters)'
fi
echo shar: extracting "'Makefile.minix'" '(445 characters)'
if test -f 'Makefile.minix'
then
echo shar: over-writing existing file "'Makefile.minix'"
fi
sed 's/^X//' << \SHAR_EOF > 'Makefile.minix'
XCROSSDIR = /dsrg/bammi/cross-minix
XCROSSLIB = $(CROSSDIR)/lib
XCROSSBIN = $(CROSSDIR)/bin
XCROSSINC = $(CROSSDIR)/include
X
XALL = libpml.a libpml32.a
X
Xall : $(ALL)
X
Xlibpml.a :
X make -f Makefile.16 clean
X make -f Makefile.16 libpml.a
X
Xlibpml32.a :
X make -f Makefile.32 clean
X make -f Makefile.32 libpml32.a
X
Xinstall : $(ALL)
X cp libpml.a libpml32.a $(CROSSLIB)
X cp pmluser.h $(CROSSINC)
X
Xclean:
X make -f Makefile.16 clean
X make -f Makefile.32 clean
SHAR_EOF
if test 445 -ne "`wc -c 'Makefile.minix'`"
then
echo shar: error transmitting "'Makefile.minix'" '(should have been 445 characters)'
fi
echo shar: extracting "'Makefile.sun'" '(1428 characters)'
if test -f 'Makefile.sun'
then
echo shar: over-writing existing file "'Makefile.sun'"
fi
sed 's/^X//' << \SHAR_EOF > 'Makefile.sun'
X#
X# FILE
X#
X# Makefile build and install the pml library (32 bit ints)
X#
X# SYNOPSIS
X#
X# make funcs make version of library in local directory
X# make install install the library (must be root)
X#
X# WARNING
X#
X# The order of the modules listed in the "LEVEL<n>" macros
X# is significant since these are the orders in which
X# they will be loaded into the library archive. Although
X# some machines support randomly ordered libraries, ordering
X# them correctly doesn't hurt...
X
XAR = ar
XCC = cc
XCFLAGS = -O -DIEEE -DNO_DBUG -I.
X
XLIB = $(CROSSLIB)
X
XINC = $(CROSSINC)
X
X
XLEVEL0 = matherr.o
X
XLEVEL1 = sign.o mod.o poly.o dabs.o sqrt.o
X
XLEVEL2 = acos.o acosh.o asin.o asinh.o atan2.o atan.o \
X log10.o tan.o \
X tanh.o sinh.o cosh.o atanh.o \
X log.o sin.o cos.o exp.o max.o min.o
X
XLEVEL3 = casin.o cacos.o cmult.o catan.o ccosh.o clog.o \
X crcp.o csinh.o csqrt.o ctan.o ctanh.o cexp.o \
X ccos.o csin.o cdiv.o csubt.o cabs.o
X
XLEVEL4 = cadd.o
X
XOBJ = $(LEVEL4) $(LEVEL3) $(LEVEL2) $(LEVEL1) $(LEVEL0)
X
X#
X# The default thing to make.
X#
X
Xdefault: libpml.a
X
Xlibpml.a: $(OBJ)
X rm -f libpml.a
X $(AR) rv libpml.a $(OBJ)
X ranlib libpml.a
X
X$(OBJ) : pml.h
X
X#
X# Stuff to do installation
X#
X
Xinstall : $(LIB)/libpml.a $(INC)/pmluser.h
X
X$(LIB)/libpml.a: libpml.a
X cp libpml.a $(LIB)/libpml.a
X
X$(INC)/pmluser.h : pmluser.h
X cp pmluser.h $(INC)/pmluser.h
X
X#
X# Clean up the directory.
X#
X
Xclean:
X rm -f *.o *.BAK libpml.a *.tmp *.bak nohup.out
SHAR_EOF
if test 1428 -ne "`wc -c 'Makefile.sun'`"
then
echo shar: error transmitting "'Makefile.sun'" '(should have been 1428 characters)'
fi
echo shar: extracting "'Makefile.tos'" '(1552 characters)'
if test -f 'Makefile.tos'
then
echo shar: over-writing existing file "'Makefile.tos'"
fi
sed 's/^X//' << \SHAR_EOF > 'Makefile.tos'
X#
X# FILE
X#
X# Makefile build and install the pml library (tos-gcc)
X#
X# SYNOPSIS
X#
X# make funcs make version of library in local directory
X# make install install the library (must be root)
X#
X# WARNING
X#
X# The order of the modules listed in the "LEVEL<n>" macros
X# is significant since these are the orders in which
X# they will be loaded into the library archive. Although
X# some machines support randomly ordered libraries, ordering
X# them correctly doesn't hurt...
XCROSSDIR = /dsrg/bammi/cross-gcc
XCROSSLIB = $(CROSSDIR)/lib
XCROSSBIN = $(CROSSDIR)/bin
XCROSSINC = $(CROSSDIR)/include
X
XAR = $(CROSSBIN)/car
XCC = $(CROSSBIN)/cgcc
XCFLAGS = -O -DIEEE -DNO_DBUG -I.
X
XLIB = $(CROSSLIB)
X
XINC = $(CROSSINC)
X
X
XLEVEL0 = matherr.o
X
XLEVEL1 = sign.o mod.o poly.o dabs.o sqrt.o
X
XLEVEL2 = acos.o acosh.o asin.o asinh.o atan2.o atan.o \
X log10.o tan.o \
X tanh.o sinh.o cosh.o atanh.o \
X log.o sin.o cos.o exp.o max.o min.o
X
XLEVEL3 = casin.o cacos.o cmult.o catan.o ccosh.o clog.o \
X crcp.o csinh.o csqrt.o ctan.o ctanh.o cexp.o \
X ccos.o csin.o cdiv.o csubt.o cabs.o
X
XLEVEL4 = cadd.o
X
XOBJ = $(LEVEL4) $(LEVEL3) $(LEVEL2) $(LEVEL1) $(LEVEL0)
X
X#
X# The default thing to make.
X#
X
Xdefault: libpml.a
X
Xlibpml.a: $(OBJ)
X rm -f libpml.a
X $(AR) srv libpml.a $(OBJ)
X
X$(OBJ) : pml.h
X
X#
X# Stuff to do installation
X#
X
Xinstall : $(LIB)/libpml.a $(INC)/pmluser.h
X
X$(LIB)/libpml.a: libpml.a
X cp libpml.a $(LIB)/libpml.a
X
X$(INC)/pmluser.h : pmluser.h
X cp pmluser.h $(INC)/pmluser.h
X
X#
X# Clean up the directory.
X#
X
Xclean:
X rm -f *.o *.BAK libpml.a *.tmp *.bak nohup.out
SHAR_EOF
if test 1552 -ne "`wc -c 'Makefile.tos'`"
then
echo shar: error transmitting "'Makefile.tos'" '(should have been 1552 characters)'
fi
echo shar: extracting "'README'" '(653 characters)'
if test -f 'README'
then
echo shar: over-writing existing file "'README'"
fi
sed 's/^X//' << \SHAR_EOF > 'README'
XThis the src directory for PML (Fred Fish's Portable Math Library)
X
XBefore you make the library
X - You must have GCC V1.31
X - You should have applied all patches to the gcc library
X (the ones that came with this distribution, and the one
X that came with the gcc V1.31 upgrade).
X - You should have compiled the library again.
X
XTo make and install the library
X GCC-TOS : use Makefile.tos
X GCC-MINIX : use Makefile.minix
XIn both Makefiles adjust CC and CROSSDIR to suit your environment.
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 653 -ne "`wc -c 'README'`"
then
echo shar: error transmitting "'README'" '(should have been 653 characters)'
fi
echo shar: extracting "'dabs.c'" '(1653 characters)'
if test -f 'dabs.c'
then
echo shar: over-writing existing file "'dabs.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'dabs.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 * dabs,fabs double precision absolute value
X *
X * KEY WORDS
X *
X * dabs
X * fabs
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns absolute value of double precision number.
X *
X * The fabs routine is supplied for compatibility with unix
X * libraries where for some bizarre reason, the double precision
X * absolute value routine is called fabs.
X *
X * USAGE
X *
X * double dabs (x)
X * double x;
X *
X * double fabs (x)
X * double x;
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "dabs";
X
X
Xdouble dabs (x)
Xdouble x;
X{
X DBUG_ENTER (funcname);
X DBUG_3 ("dabsin", "arg %le", x);
X if (x < 0.0) {
X x = -x;
X }
X DBUG_3 ("dabsout", "result %le", x);
X DBUG_RETURN (x);
X}
X
X
Xdouble fabs (x)
Xdouble x;
X{
X return (dabs(x));
X}
SHAR_EOF
if test 1653 -ne "`wc -c 'dabs.c'`"
then
echo shar: error transmitting "'dabs.c'" '(should have been 1653 characters)'
fi
echo shar: extracting "'exp.c'" '(5456 characters)'
if test -f 'exp.c'
then
echo shar: over-writing existing file "'exp.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'exp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * exp double precision exponential
X *
X * KEY WORDS
X *
X * exp
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision exponential of double precision
X * floating point number.
X *
X * USAGE
X *
X * double exp (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 96-104.
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.4). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes exponential from:
X *
X * exp(x) = 2**y * 2**z * 2**w
X *
X * Where:
X *
X * y = int ( x * log2(e) )
X *
X * v = 16 * frac ( x * log2(e))
X *
X * z = (1/16) * int (v)
X *
X * w = (1/16) * frac (v)
X *
X * Note that:
X *
X * 0 =< v < 16
X *
X * z = {0, 1/16, 2/16, ...15/16}
X *
X * 0 =< w < 1/16
X *
X * Then:
X *
X * 2**z is looked up in a table of 2**0, 2**1/16, ...
X *
X * 2**w is computed from an approximation:
X *
X * 2**w = (A + B) / (A - B)
X *
X * A = C + (D * w * w)
X *
X * B = w * (E + (F * w * w))
X *
X * C = 20.8137711965230361973
X *
X * D = 1.0
X *
X * E = 7.2135034108448192083
X *
X * F = 0.057761135831801928
X *
X * Coefficients are from HART, table #1121, pg 206.
X *
X * Effective multiplication by 2**y is done by a
X * floating point scale with y as scale argument.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X# define C 20.8137711965230361973 /* Polynomial approx coeff. */
X# define D 1.0 /* Polynomial approx coeff. */
X# define E 7.2135034108448192083 /* Polynomial approx coeff. */
X# define F 0.057761135831801928 /* Polynomial approx coeff. */
X
X
X/************************************************************************
X * *
X * This table is fixed in size and reasonably hardware *
X * independent. The given constants were generated on a *
X * DECSYSTEM 20 using double precision FORTRAN. *
X * *
X ************************************************************************
X */
X
Xstatic double fpof2tbl[] = {
X 1.00000000000000000000, /* 2 ** 0/16 */
X 1.04427378242741384020, /* 2 ** 1/16 */
X 1.09050773266525765930, /* 2 ** 2/16 */
X 1.13878863475669165390, /* 2 ** 3/16 */
X 1.18920711500272106640, /* 2 ** 4/16 */
X 1.24185781207348404890, /* 2 ** 5/16 */
X 1.29683955465100966610, /* 2 ** 6/16 */
X 1.35425554693689272850, /* 2 ** 7/16 */
X 1.41421356237309504880, /* 2 ** 8/16 */
X 1.47682614593949931110, /* 2 ** 9/16 */
X 1.54221082540794082350, /* 2 ** 10/16 */
X 1.61049033194925430820, /* 2 ** 11/16 */
X 1.68179283050742908600, /* 2 ** 12/16 */
X 1.75625216037329948340, /* 2 ** 13/16 */
X 1.83400808640934246360, /* 2 ** 14/16 */
X 1.91520656139714729380 /* 2 ** 15/16 */
X};
X
Xstatic char funcname[] = "exp";
X
X
Xdouble exp (x)
Xdouble x;
X{
X register int y;
X register int index;
X auto double w;
X auto double v;
X auto double a;
X auto double b;
X auto double t;
X auto double temp;
X auto double wpof2;
X auto double zpof2;
X auto double rtnval;
X extern double dabs ();
X extern double ldexp ();
X extern double modf ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("expin", "arg %le", x);
X if (x > LOGE_MAXDOUBLE) {
X xcpt.type = OVERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = MAXDOUBLE;
X }
X } else if (x <= LOGE_MINDOUBLE) {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: OVERFLOW error\n", funcname);
X errno = ERANGE;
X xcpt.retval = 0.0;
X }
X } else {
X t = x * LOG2E;
X v = 16.0 * modf (t, &temp);
X y = temp;
X w = modf (v, &temp) / 16.0;
X index = temp;
X
X if (x < 0.0) {
X zpof2 = 1.0 / fpof2tbl[-index];
X } else {
X zpof2 = fpof2tbl[index];
X }
X a = C + (D * w * w);
X b = w * (E + (F * w * w));
X wpof2 = (a + b) / (a - b);
X xcpt.retval = ldexp ((wpof2 * zpof2), y);
X }
X DBUG_3 ("expout", "result %le", rtnval);
X DBUG_RETURN (xcpt.retval);
X}
SHAR_EOF
if test 5456 -ne "`wc -c 'exp.c'`"
then
echo shar: error transmitting "'exp.c'" '(should have been 5456 characters)'
fi
echo shar: extracting "'log.c'" '(4649 characters)'
if test -f 'log.c'
then
echo shar: over-writing existing file "'log.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'log.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * log double precision natural log
X *
X * KEY WORDS
X *
X * log
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision natural log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 105-111.
X *
X * RESTRICTIONS
X *
X * The absolute error in the approximating polynomial is
X * 10**(-19.38). Note that contrary to DEC's assertion
X * in the F4P user's guide, the error is absolute and not
X * relative.
X *
X * This error bound assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes log(X) from:
X *
X * (1) If argument is zero then flag an error
X * and return minus infinity (or rather its
X * machine representation).
X *
X * (2) If argument is negative then flag an
X * error and return minus infinity.
X *
X * (3) Given that x = m * 2**k then extract
X * mantissa m and exponent k.
X *
X * (4) Transform m which is in range [0.5,1.0]
X * to range [1/sqrt(2), sqrt(2)] by:
X *
X * s = m * sqrt(2)
X *
X * (5) Compute z = (s - 1) / (s + 1)
X *
X * (6) Now use the approximation from HART
X * page 111 to find log(s):
X *
X * log(s) = z * ( P(z**2) / Q(z**2) )
X *
X * Where:
X *
X * P(z**2) = SUM [ Pj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * Q(z**2) = SUM [ Qj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = -0.240139179559210509868484e2
X * P1 = 0.30957292821537650062264e2
X * P2 = -0.96376909336868659324e1
X * P3 = 0.4210873712179797145
X * Q0 = -0.120069589779605254717525e2
X * Q1 = 0.19480966070088973051623e2
X * Q2 = -0.89111090279378312337e1
X * Q3 = 1.0000
X *
X * (coefficients from HART table #2705 pg 229)
X *
X * (7) Finally, compute log(x) from:
X *
X * log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double log_pcoeffs[] = {
X -0.24013917955921050986e2,
X 0.30957292821537650062e2,
X -0.96376909336868659324e1,
X 0.4210873712179797145
X};
X
Xstatic double log_qcoeffs[] = {
X -0.12006958977960525471e2,
X 0.19480966070088973051e2,
X -0.89111090279378312337e1,
X 1.0000
X};
X
Xstatic char funcname[] = "log";
X
X
Xdouble log (x)
Xdouble x;
X{
X auto int k;
X auto double s;
X auto double z;
X auto double zt2;
X auto double pqofz;
X auto struct exception xcpt;
X extern double frexp ();
X extern double poly ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("login", "arg %le", x);
X if (x == 0.0) {
X xcpt.type = SING;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: SINGULARITY error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else {
X s = SQRT2 * frexp (x, &k);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "s = %le", s);
X z = (s - 1.0) / (s + 1.0);
X DBUG_3 ("log", "z = %le", z);
X zt2 = z * z;
X DBUG_3 ("log", "zt2 = %le", zt2);
X pqofz = z * (poly (3, log_pcoeffs, zt2) / poly (3, log_qcoeffs, zt2));
X DBUG_3 ("pqofz", "pqofz = %le", pqofz);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "LN2 = %le", LN2);
X x = k * LN2;
X DBUG_3 ("log", "x = %le", x);
X x -= LNSQRT2;
X DBUG_3 ("log", "x = %le", x);
X x += pqofz;
X DBUG_3 ("log", "x = %le", x);
X xcpt.retval = x;
X }
X DBUG_3 ("logout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
SHAR_EOF
if test 4649 -ne "`wc -c 'log.c'`"
then
echo shar: error transmitting "'log.c'" '(should have been 4649 characters)'
fi
echo shar: extracting "'log10.c'" '(1744 characters)'
if test -f 'log10.c'
then
echo shar: over-writing existing file "'log10.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'log10.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 * log10 double precision common log
X *
X * KEY WORDS
X *
X * log10
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision common log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log10 (x)
X * double x;
X *
X * REFERENCES
X *
X * PDP-11 Fortran IV-plus users guide, Digital Equip. Corp.,
X * 1975, pp. B-3.
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 log10(x) from:
X *
X * log10(x) = log10(e) * log(x)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "log10";
X
X
Xdouble log10 (x)
Xdouble x;
X{
X extern double log ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("log10in", "arg %le", x);
X x = LOG10E * log (x);
X DBUG_3 ("log10out", "result %le", x);
X DBUG_RETURN (x);
X}
SHAR_EOF
if test 1744 -ne "`wc -c 'log10.c'`"
then
echo shar: error transmitting "'log10.c'" '(should have been 1744 characters)'
fi
echo shar: extracting "'matherr.c'" '(1609 characters)'
if test -f 'matherr.c'
then
echo shar: over-writing existing file "'matherr.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'matherr.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 * matherr default math error handler function
X *
X * KEY WORDS
X *
X * error handler
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Default math error handler for the math library. This routine
X * may be replaced by one of the user's own handlers. The default
X * is to do nothing and returns zero. If the user wishes to supply
X * the return value for the function, and to suppress default
X * error handling by the function, his matherr routine must return
X * non-zero.
X *
X * USAGE
X *
X * int matherr (xcpt)
X * struct exception *xcpt;
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic char funcname[] = "matherr";
X
Xint matherr (xcpt)
Xstruct exception *xcpt;
X{
X DBUG_ENTER (funcname);
X DBUG_RETURN (0);
X}
SHAR_EOF
if test 1609 -ne "`wc -c 'matherr.c'`"
then
echo shar: error transmitting "'matherr.c'" '(should have been 1609 characters)'
fi
echo shar: extracting "'max.c'" '(1411 characters)'
if test -f 'max.c'
then
echo shar: over-writing existing file "'max.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'max.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 * max double precision maximum of two arguments
X *
X * KEY WORDS
X *
X * max
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns maximum value of two double precision numbers.
X *
X * USAGE
X *
X * double max (x,y)
X * double x;
X * double y;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble max (x, y)
Xdouble x;
Xdouble y;
X{
X ENTER ("max");
X DEBUG4 ("maxin", "x = %le y = %le", x, y);
X if (x < y) {
X x = y;
X }
X DEBUG3 ("maxout", "result %le", x);
X LEAVE ();
X return (x);
X}
SHAR_EOF
if test 1411 -ne "`wc -c 'max.c'`"
then
echo shar: error transmitting "'max.c'" '(should have been 1411 characters)'
fi
echo shar: extracting "'min.c'" '(1412 characters)'
if test -f 'min.c'
then
echo shar: over-writing existing file "'min.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'min.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 * min double precision minimum of two arguments
X *
X * KEY WORDS
X *
X * min
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns minimum value of two double precision numbers.
X *
X * USAGE
X *
X * double min (x, y)
X * double x;
X * double y;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble min (x, y)
Xdouble x;
Xdouble y;
X{
X ENTER ("min");
X DEBUG4 ("minin", "x = %le y = %le", x, y);
X if (x > y) {
X x = y;
X }
X DEBUG3 ("minout", "result %le", x);
X LEAVE ();
X return (x);
X}
SHAR_EOF
if test 1412 -ne "`wc -c 'min.c'`"
then
echo shar: error transmitting "'min.c'" '(should have been 1412 characters)'
fi
echo shar: extracting "'mod.c'" '(1520 characters)'
if test -f 'mod.c'
then
echo shar: over-writing existing file "'mod.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'mod.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 * mod double precision modulo
X *
X * KEY WORDS
X *
X * mod
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision modulo of two double
X * precision arguments.
X *
X * USAGE
X *
X * double mod (value, base)
X * double value;
X * double base;
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xdouble mod (value, base)
Xdouble value;
Xdouble base;
X{
X auto double intpart;
X extern double modf ();
X
X DBUG_ENTER ("mod");
X DBUG_4 ("modin", "args %le %le", value, base);
X value /= base;
X value = modf (value, &intpart);
X value *= base;
X DBUG_3 ("modout", "result %le", value);
X DBUG_RETURN (value);
X}
SHAR_EOF
if test 1520 -ne "`wc -c 'mod.c'`"
then
echo shar: error transmitting "'mod.c'" '(should have been 1520 characters)'
fi
echo shar: extracting "'pml.h'" '(4239 characters)'
if test -f 'pml.h'
then
echo shar: over-writing existing file "'pml.h'"
fi
sed 's/^X//' << \SHAR_EOF > 'pml.h'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * This file gets included with all of the floating point math
X * library routines when they are compiled. Note that
X * this is the proper place to put machine dependencies
X * whenever possible.
X *
X * It should be pointed out that for simplicity's sake, the
X * environment parameters are defined as floating point constants,
X * rather than octal or hexadecimal initializations of allocated
X * storage areas. This means that the range of allowed numbers
X * may not exactly match the hardware's capabilities. For example,
X * if the maximum positive double precision floating point number
X * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
X * defined to be 1.11E100 then the numbers between 1.11E100 and
X * 1.11...E100 are considered to be undefined. For most
X * applications, this will cause no problems.
X *
X * An alternate method is to allocate a global static "double" variable,
X * say "maxdouble", and use a union declaration and initialization
X * to initialize it with the proper bits for the EXACT maximum value.
X * This was not done because the only compilers available to the
X * author did not fully support union initialization features.
X *
X */
X
X#ifndef NO_DBUG
X# include <dbug.h>
X#else
X# define DBUG_ENTER(a1)
X# define DBUG_RETURN(a1) return(a1)
X# define DBUG_VOID_RETURN return
X# define DBUG_EXECUTE(keyword,a1)
X# define DBUG_2(keyword,format)
X# define DBUG_3(keyword,format,a1)
X# define DBUG_4(keyword,format,a1,a2)
X# define DBUG_5(keyword,format,a1,a2,a3)
X# define DBUG_PUSH(a1)
X# define DBUG_POP()
X# define DBUG_PROCESS(a1)
X# define DBUG_FILE (stderr)
X# define ENTER(a1)
X# define LEAVE()
X# define DEBUG3(keyword,format,a1)
X# define DEBUG4(keyword,format,a1,a2)
X# define DEBUGPUSH(a)
X# define DEBUGWHO(w)
X#endif
X
X#include <errno.h>
Xextern int errno;
X#ifdef atarist /* #defied automatically by gcc-tos only */
X#define EDOM EERROR
X#define ERANGE EERROR
X#endif
X
X/*
X * MAXDOUBLE => Maximum double precision number
X * MINDOUBLE => Minimum double precision number
X * DMAXEXP => Maximum exponent of a double
X * DMINEXP => Minimum exponent of a double
X */
X
X#define MAXDOUBLE 1.7e+308
X#define MINDOUBLE 3.0e-308
X#define DMAXEXP 1023
X#define DMINEXP (-1022)
X
X#define LOG2_MAXDOUBLE 1024
X#define LOG2_MINDOUBLE (-1023)
X#define LOGE_MAXDOUBLE 7.09782712893383970e+02
X#define LOGE_MINDOUBLE -7.09089565712824080e+02
X
X/*
X * The following are hacks which should be fixed when I understand all
X * the issues a little better. |tanh(TANH_MAXARG)| = 1.0
X */
X#define TANH_MAXARG 16
X#define SQRT_MAXDOUBLE 1.304380e19
X
X#define PI 3.14159265358979323846
X#define TWOPI 6.28318530717958620
X#define HALFPI 1.57079632679489660
X#define FOURTHPI 0.785398163397448280
X#define SIXTHPI 0.523598775598298820
X
X#define LOG2E 1.4426950408889634074 /* Log to base 2 of e */
X#define LOG10E 0.4342944819032518276
X#define SQRT2 1.41421356237309504880
X#define SQRT3 1.7320508075688772935
X#define LN2 0.69314718055994530942
X#define LNSQRT2 0.3465735902799726547
X
X
X/*
X * MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X */
X
X#ifdef IEEE
X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */
X#endif
X
X#ifdef TRUE
X#undef TRUE
X#endif
X
X#ifdef FALSE
X#undef FALSE
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define VOID void
SHAR_EOF
if test 4239 -ne "`wc -c 'pml.h'`"
then
echo shar: error transmitting "'pml.h'" '(should have been 4239 characters)'
fi
echo shar: extracting "'pmlerr.c'" '(7682 characters)'
if test -f 'pmlerr.c'
then
echo shar: over-writing existing file "'pmlerr.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'pmlerr.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/************************************************************************
X * *
X * PORTABLE MATH LIBRARY -- ERROR HANDLER *
X * *
X * This is a sample PML library error handler. *
X * It may be used as is, or another of the user's choice *
X * substituted. *
X * *
X * In any case, the global "pmlerr" must be defined somewhere *
X * in the user task, since many Portable Math Library routines *
X * reference it. The other routines in this file are not called *
X * by any library routines and may be omitted. *
X * *
X ************************************************************************
X */
X
X# include <stdio.h>
X# include <pmluser.h>
X# include "pml.h"
X
Xstatic struct pml_err {
X int flags; /* Flag word; bits defined in pml.h */
X char *msg; /* Error message */
X char *func; /* Function in which error occured */
X};
X
Xstatic struct pml_err pml_errs[] = {
X CONTINUE | COUNT | LOG, "overflow", "exp",
X CONTINUE | COUNT | LOG, "underflow", "exp",
X CONTINUE | COUNT | LOG, "exponent overflow", "scale",
X CONTINUE | COUNT | LOG, "negative argument", "sqrt",
X CONTINUE | COUNT | LOG, "zero argument", "log",
X CONTINUE | COUNT | LOG, "negative argument", "log",
X CONTINUE | COUNT | LOG, "argument magnitude greater than 1.0", "acos",
X CONTINUE | COUNT | LOG, "argument magnitude greater than 1.0", "asin",
X CONTINUE | COUNT | LOG, "overflow", "tan",
X CONTINUE | COUNT | LOG, "overflow", "cosh",
X CONTINUE | COUNT | LOG, "underflow", "cosh",
X CONTINUE | COUNT | LOG, "overflow", "sinh",
X CONTINUE | COUNT | LOG, "underflow", "sinh",
X CONTINUE | COUNT | LOG, "overflow", "asinh",
X CONTINUE | COUNT | LOG, "argument less than 1.0", "acosh",
X CONTINUE | COUNT | LOG, "overflow", "acosh",
X CONTINUE | COUNT | LOG, "argument magnitude not < 1.0", "atanh",
X CONTINUE | COUNT | LOG, "underflow", "atan",
X CONTINUE | COUNT | LOG, "complex division by zero", "cdiv",
X CONTINUE | COUNT | LOG, "complex reciprocal of zero", "crcp",
X CONTINUE | COUNT | LOG, "exponent underflow", "scale",
X CONTINUE | COUNT | LOG, "argument has no fractional part", "dint",
X};
X
Xstatic int err_count = 0; /* Counter for PML errors */
Xstatic int err_limit = MAX_ERRORS; /* PML error limit */
X
X/*
X * FUNCTION
X *
X * pmlcfs Clear specified PML error handler flags
X *
X * KEY WORDS
X *
X * pmlcfs
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Clear the specified PML error handler flags for the
X * specified error. Two or more flags may be cleared simultaneously
X * by "or-ing" them in the call, for example "LOG | CONTINUE".
X * The manifest constants for the flags and error codes are
X * defined in <pmluser.h>.
X *
X * USAGE
X *
X * pmlcfs(err_code,flags)
X * int err_code;
X * int flags;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xpmlcfs(err_code,flag_word)
Xregister int err_code;
Xregister int flag_word;
X{
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"pmlcfs: invalid error code %d\n",err_code);
X } else {
X pml_errs[err_code].flags &= ~flag_word;
X }
X}
X
X
X/*
X * FUNCTION
X *
X * pmlcnt get PML error count and reset it to zero
X *
X * KEY WORDS
X *
X * pmlcnt
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns the total number of PML errors seen
X * prior to the call, and resets the error count to zero.
X *
X * USAGE
X *
X * int pmlcnt()
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xint pmlcnt()
X{
X register int rtn_val;
X
X rtn_val = err_count;
X err_count = 0;
X return(rtn_val);
X}
X
X
X/*
X * FUNCTION
X *
X * pmlerr Portable Math Library error handler
X *
X * KEY WORDS
X *
X * pmlerr
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Provides a sample PML error handler. Does
X * not use any available hardware "traps" so is machine
X * independent. Generally only called internally by the
X * other PML routines.
X *
X * There are currently three flags which control the
X * response for specific errors:
X *
X * (1) LOG When set an error message is sent
X * to the user terminal.
X *
X * (2) COUNT When set the error is counted
X * against the PML error limit.
X *
X * (3) CONTINUE When set the task continues
X * providing the error count has not
X * exceeded the PML error limit.
X *
X * Each of these flags can be set or reset independently
X * by "pmlsfs" or "pmlcfs" respectively.
X *
X * USAGE
X *
X * pmlerr(err_code)
X * int err_code;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xpmlerr(err_code)
Xregister int err_code;
X{
X register struct pml_err *err;
X
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"pmlerr: invalid error code %d\n",err_code);
X } else {
X err = &pml_errs[err_code];
X if (err->flags & LOG) {
X fprintf(stderr,"pml: %s in function \"%s\"\n",err->msg,err->func);
X }
X if (err->flags & COUNT) {
X err_count++;
X }
X if ((err->flags & CONTINUE) && (err_count <= err_limit)) {
X return;
X } else {
X fprintf(stderr,"pml: error limit exceeded\n");
X fprintf(stderr,"pml: task aborted with %d error(s)\n",
X err_count);
X exit(-1);
X }
X }
X}
X
X
X/*
X * FUNCTION
X *
X * pmllim Set Portable Math Library error limit
X *
X * KEY WORDS
X *
X * pmllim
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Sets the PML error limit to the specified value
X * and returns it previous value.
X * Does not affect the current error count (which may be reset
X * to zero by a call to "pmlcnt"). Note that the default error
X * limit is set at compile time by the value in "pml.h".
X *
X * USAGE
X *
X * int pmllim(limit)
X * int limit;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
Xint pmllim(limit)
Xregister int limit;
X{
X register int rtn_val;
X
X rtn_val = err_limit;
X err_limit = limit;
X return(rtn_val);
X}
X
X
X/*
X * FUNCTION
X *
X * pmlsfs Set specified PML error handler flags
X *
X * KEY WORDS
X *
X * pmlsfs
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Set the specified PML error handler flags for the
X * specified error. Two or more flags may be set simultaneously
X * by "or-ing" them in the call, for example "LOG | CONTINUE".
X * The manifest constants for the flags and error codes are
X * defined in <pmluser.h>.
X *
X * USAGE
X *
X * pmlsfs(err_code,flags)
X * int err_code;
X * int flags;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X *
X */
X
Xpmlsfs(err_code,flag_word)
Xregister int err_code;
Xregister int flag_word;
X{
X if (err_code < 0 || err_code > (sizeof(pml_errs)/sizeof(struct pml_err))) {
X fprintf(stderr,"? pmlsfs --- invalid error code %d.\n",err_code);
X } else {
X pml_errs[err_code].flags |= flag_word;
X }
X}
SHAR_EOF
if test 7682 -ne "`wc -c 'pmlerr.c'`"
then
echo shar: error transmitting "'pmlerr.c'" '(should have been 7682 characters)'
fi
echo shar: extracting "'pmluser.h'" '(3251 characters)'
if test -f 'pmluser.h'
then
echo shar: over-writing existing file "'pmluser.h'"
fi
sed 's/^X//' << \SHAR_EOF > 'pmluser.h'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FILE
X *
X * pmluser.h include file for users of portable math library
X *
X * SYNOPSIS
X *
X * #include <pmluser.h>
X *
X * DESCRIPTION
X *
X * This file should be included in any user compilation module
X * which accesses routines from the Portable Math Library (PML).
X *
X */
X
X
X/*
X * Create the type "COMPLEX". This is an obvious extension that I
X * hope becomes a part of standard C someday.
X *
X */
X
Xtypedef struct cmplx { /* Complex structure */
X double real; /* Real part */
X double imag; /* Imaginary part */
X} COMPLEX;
X
X/* exceptions ++jrb */
Xtypedef enum { DOMAIN, OVERFLOW, UNDERFLOW, PLOSS, SING } exception_type;
X
Xstruct exception {
X exception_type type; /* exception type */
X char *name; /* function in which it occured */
X double arg1; /* an arg */
X double retval; /* val to return */
X};
X
X#ifdef __GNUC__
Xdouble acos(double);
Xdouble acosh(double);
Xdouble asin(double);
Xdouble asinh(double);
Xdouble atan(double);
Xdouble atan2(double, double);
Xdouble atanh(double);
Xdouble cabs(COMPLEX);
Xdouble cos(double);
Xdouble cosh(double);
Xdouble dabs(double);
Xdouble fabs(double);
Xdouble exp(double);
Xdouble log(double);
Xdouble log10(double);
Xint matherr(struct exception *);
Xdouble max(double, double);
Xdouble min(double, double);
Xdouble mod(double, double);
Xint pmlcfs(int, int);
Xint pmlcnt();
Xint pmlerr(int);
Xint pmllim(int);
Xint pmlsfs(int, int);
Xdouble poly(int, double *, double);
Xdouble sign(double, double);
Xdouble sin(double);
Xdouble sinh(double);
Xdouble sqrt(double);
Xdouble tan(double);
Xdouble tanh(double);
X
X#else
X
Xdouble acos(/* double */);
Xdouble acosh(/* double */);
Xdouble asin(/* double */);
Xdouble asinh(/* double */);
Xdouble atan(/* double */);
Xdouble atan2(/* double, double */);
Xdouble atanh(/* double */);
Xdouble cabs(/* COMPLEX */);
Xdouble cos(/* double */);
Xdouble cosh(/* double */);
Xdouble dabs(/* double */);
Xdouble fabs(/* double */);
Xdouble exp(/* double */);
Xdouble log(/* double */);
Xdouble log10(/* double */);
Xint matherr(/* struct exception * */);
Xdouble max(/* double, double */);
Xdouble min(/* double, double */);
Xdouble mod(/* double, double */);
Xint pmlcfs(/* int, int */);
Xint pmlcnt(/* */);
Xint pmlerr(/* int */);
Xint pmllim(/* int */);
Xint pmlsfs(/* int, int */);
Xdouble poly(/* int, double *, double */);
Xdouble sign(/* double, double */);
Xdouble sin(/* double */);
Xdouble sinh(/* double */);
Xdouble sqrt(/* double */);
Xdouble tan(/* double */);
Xdouble tanh(/* double */);
X
X#endif /* __GNUC__ */
SHAR_EOF
if test 3251 -ne "`wc -c 'pmluser.h'`"
then
echo shar: error transmitting "'pmluser.h'" '(should have been 3251 characters)'
fi
echo shar: extracting "'poly.c'" '(2049 characters)'
if test -f 'poly.c'
then
echo shar: over-writing existing file "'poly.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2049 -ne "`wc -c 'poly.c'`"
then
echo shar: error transmitting "'poly.c'" '(should have been 2049 characters)'
fi
echo shar: extracting "'sign.c'" '(1580 characters)'
if test -f 'sign.c'
then
echo shar: over-writing existing file "'sign.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'sign.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 * sign transfer of sign
X *
X * KEY WORDS
X *
X * sign
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns first argument with same sign as second argument.
X *
X * USAGE
X *
X * double sign (x, y)
X * double x;
X * double y;
X *
X * PROGRAMMER
X *
X * Fred Fish
X * Tempe, Az 85281
X * (602) 966-8871
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
Xdouble sign (x, y)
Xdouble x;
Xdouble y;
X{
X double rtnval;
X
X ENTER ("sign");
X DEBUG4 ("signin", "args %le %le", x, y);
X if (x > 0.0) {
X if (y > 0.0) {
X rtnval = x;
X } else {
X rtnval = -x;
X }
X } else {
X if (y < 0.0) {
X rtnval = x;
X } else {
X rtnval = -x;
X }
X }
X DEBUG3 ("signout", "result %le", rtnval);
X LEAVE ();
X return (rtnval);
X}
SHAR_EOF
if test 1580 -ne "`wc -c 'sign.c'`"
then
echo shar: error transmitting "'sign.c'" '(should have been 1580 characters)'
fi
echo shar: extracting "'sin.c'" '(4983 characters)'
if test -f 'sin.c'
then
echo shar: over-writing existing file "'sin.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'sin.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sin double precision sine
X *
X * KEY WORDS
X *
X * sin
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision sine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sin (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes sin(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call sin recursively
X * using relation sin(x) = -sin(x - PI).
X *
X * (3) If x < -PI/2 then call sin recursively
X * using relation sin(x) = -sin(x + PI).
X *
X * (4) If x > PI/4 then call cos using
X * relation sin(x) = cos(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation sin(x) = -cos(PI/2 + x).
X *
X * (6) If x is small enough that polynomial
X * evaluation would cause underflow
X * then return x, since sin(x)
X * approaches x as x approaches zero.
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 118 can be used:
X *
X * sin(x) = y * ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.206643433369958582409167054e+7
X * P1 = -0.18160398797407332550219213e+6
X * P2 = 0.359993069496361883172836e+4
X * P3 = -0.2010748329458861571949e+2
X * Q0 = 0.263106591026476989637710307e+7
X * Q1 = 0.3927024277464900030883986e+5
X * Q2 = 0.27811919481083844087953e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3063 pg 234)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double sin_pcoeffs[] = {
X 0.20664343336995858240e7,
X -0.18160398797407332550e6,
X 0.35999306949636188317e4,
X -0.20107483294588615719e2
X};
X
Xstatic double sin_qcoeffs[] = {
X 0.26310659102647698963e7,
X 0.39270242774649000308e5,
X 0.27811919481083844087e3,
X 1.0
X};
X
Xstatic char funcname[] = "sin";
X
X
Xdouble sin (x)
Xdouble x;
X{
X double y;
X double yt2;
X double rtnval;
X extern double mod ();
X extern double cos ();
X extern double poly ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("sinin", "arg %le", x);
X if (x < -PI || x > PI) {
X x = mod (x, TWOPI);
X if (x > PI) {
X x = x - TWOPI;
X } else if (x < -PI) {
X x = x + TWOPI;
X }
X }
X if (x > HALFPI) {
X xcpt.retval = -(sin (x - PI));
X } else if (x < -HALFPI) {
X xcpt.retval = -(sin (x + PI));
X } else if (x > FOURTHPI) {
X xcpt.retval = cos (HALFPI - x);
X } else if (x < -FOURTHPI) {
X xcpt.retval = -(cos (HALFPI + x));
X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
X xcpt.retval = x;
X } else {
X y = x / FOURTHPI;
X yt2 = y * y;
X xcpt.retval = y * (poly (3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2));
X }
X DBUG_3 ("sinout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
SHAR_EOF
if test 4983 -ne "`wc -c 'sin.c'`"
then
echo shar: error transmitting "'sin.c'" '(should have been 4983 characters)'
fi
echo shar: extracting "'sinh.c'" '(2423 characters)'
if test -f 'sinh.c'
then
echo shar: over-writing existing file "'sinh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2423 -ne "`wc -c 'sinh.c'`"
then
echo shar: error transmitting "'sinh.c'" '(should have been 2423 characters)'
fi
echo shar: extracting "'sqrt.c'" '(4591 characters)'
if test -f 'sqrt.c'
then
echo shar: over-writing existing file "'sqrt.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'sqrt.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sqrt double precision square root
X *
X * KEY WORDS
X *
X * sqrt
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision square root of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sqrt (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 89-96.
X *
X * RESTRICTIONS
X *
X * The relative error is 10**(-30.1) after three applications
X * of Heron's iteration for the square root.
X *
X * However, this assumes exact arithmetic in the iterations
X * and initial approximation. Additional errors may occur
X * due to truncation, rounding, or machine precision limits.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes square root by:
X *
X * (1) Range reduction of argument to [0.5,1.0]
X * by application of identity:
X *
X * sqrt(x) = 2**(k/2) * sqrt(x * 2**(-k))
X *
X * k is the exponent when x is written as
X * a mantissa times a power of 2 (m * 2**k).
X * It is assumed that the mantissa is
X * already normalized (0.5 =< m < 1.0).
X *
X * (2) An approximation to sqrt(m) is obtained
X * from:
X *
X * u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
X *
X * P0 = 0.594604482
X * P1 = 2.54164041
X * Q0 = 2.13725758
X * Q1 = 1.0
X *
X * (coefficients from HART table #350 pg 193)
X *
X * (3) Three applications of Heron's iteration are
X * performed using:
X *
X * y[n+1] = 0.5 * (y[n] + (m/y[n]))
X *
X * where y[0] = u = approx sqrt(m)
X *
X * (4) If the value of k was odd then y is either
X * multiplied by the square root of two or
X * divided by the square root of two for k positive
X * or negative respectively. This rescales y
X * by multiplying by 2**frac(k/2).
X *
X * (5) Finally, y is rescaled by int(k/2) which
X * is equivalent to multiplication by 2**int(k/2).
X *
X * The result of steps 4 and 5 is that the value
X * of y between 0.5 and 1.0 has been rescaled by
X * 2**(k/2) which removes the original rescaling
X * done prior to finding the mantissa square root.
X *
X * NOTES
X *
X * The Convergent Technologies compiler optimizes division
X * by powers of two to "arithmetic shift right" instructions.
X * This is ok for positive integers but gives different
X * results than other C compilers for negative integers.
X * For example, (-1)/2 is -1 on the CT box, 0 on every other
X * machine I tried.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#define P0 0.594604482 /* Approximation coeff */
X#define P1 2.54164041 /* Approximation coeff */
X#define Q0 2.13725758 /* Approximation coeff */
X#define Q1 1.0 /* Approximation coeff */
X
X#define ITERATIONS 3 /* Number of iterations */
X
Xstatic char funcname[] = "sqrt";
X
X
Xdouble sqrt (x)
Xdouble x;
X{
X auto int k;
X register int bugfix;
X register int kmod2;
X register int count;
X auto int exponent;
X auto double m;
X auto double u;
X auto double y;
X auto double rtnval;
X auto struct exception xcpt;
X extern double frexp ();
X extern double ldexp ();
X
X DBUG_ENTER ("sqrt");
X DBUG_3 ("sqrtin", "arg %le", x);
X if (x == 0.0) {
X rtnval = 0.0;
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X } else {
X m = frexp (x, &k);
X u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
X for (count = 0, y = u; count < ITERATIONS; count++) {
X y = 0.5 * (y + (m / y));
X }
X if ((kmod2 = (k % 2)) < 0) {
X y /= SQRT2;
X } else if (kmod2 > 0) {
X y *= SQRT2;
X }
X bugfix = 2;
X xcpt.retval = ldexp (y, k/bugfix);
X }
X DBUG_3 ("sqrtout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
SHAR_EOF
if test 4591 -ne "`wc -c 'sqrt.c'`"
then
echo shar: error transmitting "'sqrt.c'" '(should have been 4591 characters)'
fi
echo shar: extracting "'tan.c'" '(2208 characters)'
if test -f 'tan.c'
then
echo shar: over-writing existing file "'tan.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2208 -ne "`wc -c 'tan.c'`"
then
echo shar: error transmitting "'tan.c'" '(should have been 2208 characters)'
fi
echo shar: extracting "'tanh.c'" '(2312 characters)'
if test -f 'tanh.c'
then
echo shar: over-writing existing file "'tanh.c'"
fi
sed 's/^X//' << \SHAR_EOF > '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}
SHAR_EOF
if test 2312 -ne "`wc -c 'tanh.c'`"
then
echo shar: error transmitting "'tanh.c'" '(should have been 2312 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