Portable Math Library in C (Part 3 of 6)
f...@mcdsun.uucp
This is a portable math library written entirely in C. Since it has been
several years since I had any interest in doing any more work on it, and
people may find it useful, I have decided to post it. There should be
a lead-in posting (part 0 of 6?) containing a README file and commands
to make the directories for the regular parts 1-6. Be sure to read the
README file in 'part 0'.
-Fred Fish
=============== Cut here and feed to the shell ========================
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 3 (of 6)."
# Contents: funcs/src/atan.c funcs/src/cos.c funcs/src/log.c
# funcs/src/pml.h funcs/src/sin.c funcs/src/sqrt.c
# funcs/unused/scale.c funcs/unused/xexp.c funcs/unused/xmant.c
# tests/c2c.dat tests/cc2c.dat
# Wrapped by fnf@mcdsun on Fri Apr 10 16:21:43 1987
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f funcs/src/atan.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/atan.c\"
else
echo shar: Extracting \"funcs/src/atan.c\" \(5187 characters\)
sed "s/^X//" >funcs/src/atan.c <<'END_OF_funcs/src/atan.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * atan double precision arc tangent
X *
X * KEY WORDS
X *
X * atan
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision arc tangent of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double atan (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 120-130.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes arctangent(x) from:
X *
X * (1) If x < 0 then negate x, perform steps
X * 2, 3, and 4, and negate the returned
X * result. This makes use of the identity
X * atan(-x) = -atan(x).
X *
X * (2) If argument x > 1 at this point then
X * test to be sure that x can be inverted
X * without underflowing. If not, reduce
X * x to largest possible number that can
X * be inverted, issue warning, and continue.
X * Perform steps 3 and 4 with arg = 1/x
X * and subtract returned result from
X * pi/2. This makes use of the identity
X * atan(x) = pi/2 - atan(1/x) for x>0.
X *
X * (3) At this point 0 <= x <= 1. If
X * x > tan(pi/12) then perform step 4
X * with arg = (x*sqrt(3)-1)/(sqrt(3)+x)
X * and add pi/6 to returned result. This
X * last transformation maps arguments
X * greater than tan(pi/12) to arguments
X * in range 0...tan(pi/12).
X *
X * (4) At this point the argument has been
X * transformed so that it lies in the
X * range -tan(pi/12)...tan(pi/12).
X * Since very small arguments may cause
X * underflow in the polynomial evaluation,
X * a final check is performed. If the
X * argument is less than the underflow
X * bound, the function returns x, since
X * atan(x) approaches asin(x) which
X * approaches x, as x goes to zero.
X *
X * (5) atan(x) is now computed by one of the
X * approximations given in the cited
X * reference (Hart). That is:
X *
X * atan(x) = x * SUM [ C[i] * x**(2*i) ]
X * over i = {0,1,...8}.
X *
X * Where:
X *
X * C[0] = .9999999999999999849899
X * C[1] = -.333333333333299308717
X * C[2] = .1999999999872944792
X * C[3] = -.142857141028255452
X * C[4] = .11111097898051048
X * C[5] = -.0909037114191074
X * C[6] = .0767936869066
X * C[7] = -.06483193510303
X * C[8] = .0443895157187
X * (coefficients from HART table #4945 pg 267)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double atan_coeffs[] = {
X .9999999999999999849899, /* P0 must be first */
X -.333333333333299308717,
X .1999999999872944792,
X -.142857141028255452,
X .11111097898051048,
X -.0909037114191074,
X .0767936869066,
X -.06483193510303,
X .0443895157187 /* Pn must be last */
X};
X
Xstatic char funcname[] = "atan";
X
X#define LAST_BOUND 0.2679491924311227074725 /* tan (PI/12) */
X
X
Xdouble atan (x)
Xdouble x;
X{
X register int order;
X double xt2;
X double t1;
X double t2;
X extern double poly ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("atanin", "arg %le", x);
X if (x < 0.0) {
X xcpt.retval = -(atan (-x));
X } else if (x > 1.0) {
X if (x < MAXDOUBLE && x > -MAXDOUBLE) {
X xcpt.retval = HALFPI - atan (1.0 / x);
X } else {
X xcpt.type = UNDERFLOW;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: UNDERFLOW error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X }
X } else if (x > LAST_BOUND) {
X t1 = x * SQRT3 - 1.0;
X t2 = SQRT3 + x;
X xcpt.retval = SIXTHPI + atan (t1 / t2);
X } else if (x < X16_UNDERFLOWS) {
X xcpt.type = PLOSS;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: PLOSS error\n", funcname);
X errno = EDOM;
X xcpt.retval = x;
X }
X } else {
X xt2 = x * x;
X order = sizeof (atan_coeffs) / sizeof(double);
X order -= 1;
X xcpt.retval = x * poly (order, atan_coeffs, xt2);
X }
X DBUG_3 ("atanout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/atan.c
if test 5187 -ne `wc -c <funcs/src/atan.c`; then
echo shar: \"funcs/src/atan.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/cos.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/cos.c\"
else
echo shar: Extracting \"funcs/src/cos.c\" \(5011 characters\)
sed "s/^X//" >funcs/src/cos.c <<'END_OF_funcs/src/cos.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * cos double precision cosine
X *
X * KEY WORDS
X *
X * cos
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double cos (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes cos(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call cos recursively
X * using relation cos(x) = -cos(x - PI).
X *
X * (3) If x < -PI/2 then call cos recursively
X * using relation cos(x) = -cos(x + PI).
X *
X * (4) If x > PI/4 then call sin using
X * relation cos(x) = sin(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation cos(x) = sin(PI/2 + x).
X *
X * (6) If x would cause underflow in approx
X * evaluation arithmetic then return
X * sqrt(1.0 - x**2).
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 119 can be used:
X *
X * cos(x) = ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.12905394659037374438571854e+7
X * P1 = -0.3745670391572320471032359e+6
X * P2 = 0.134323009865390842853673e+5
X * P3 = -0.112314508233409330923e+3
X * Q0 = 0.12905394659037373590295914e+7
X * Q1 = 0.234677731072458350524124e+5
X * Q2 = 0.2096951819672630628621e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3843 pg 244)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X *
X */
X
X# include <stdio.h>
X# include <pmluser.h>
X# include "pml.h"
X
Xstatic double cos_pcoeffs[] = {
X 0.12905394659037374438e7,
X -0.37456703915723204710e6,
X 0.13432300986539084285e5,
X -0.11231450823340933092e3
X};
X
Xstatic double cos_qcoeffs[] = {
X 0.12905394659037373590e7,
X 0.23467773107245835052e5,
X 0.20969518196726306286e3,
X 1.0
X};
X
Xstatic char funcname[] = "cos";
X
X
Xdouble cos (x)
Xdouble x;
X{
X auto double y;
X auto double yt2;
X auto double rtnval;
X extern double mod ();
X extern double sin ();
X extern double sqrt ();
X extern double poly ();
X struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("cosin", "arg %le", x);
X if (x < -PI || x > PI) {
X x = mod (x, TWOPI);
X if (x > PI) {
X x = x - TWOPI;
X } else if (x < -PI) {
X x = x + TWOPI;
X }
X }
X if (x > HALFPI) {
X xcpt.retval = -(cos (x - PI));
X } else if (x < -HALFPI) {
X xcpt.retval = -(cos (x + PI));
X } else if (x > FOURTHPI) {
X xcpt.retval = sin (HALFPI - x);
X } else if (x < -FOURTHPI) {
X xcpt.retval = sin (HALFPI + x);
X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
X xcpt.retval = sqrt (1.0 - (x * x));
X } else {
X y = x / FOURTHPI;
X yt2 = y * y;
X xcpt.retval = poly (3, cos_pcoeffs, yt2) / poly (3, cos_qcoeffs, yt2);
X }
X DBUG_3 ("cosout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/cos.c
if test 5011 -ne `wc -c <funcs/src/cos.c`; then
echo shar: \"funcs/src/cos.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/log.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/log.c\"
else
echo shar: Extracting \"funcs/src/log.c\" \(4649 characters\)
sed "s/^X//" >funcs/src/log.c <<'END_OF_funcs/src/log.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * log double precision natural log
X *
X * KEY WORDS
X *
X * log
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision natural log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 105-111.
X *
X * RESTRICTIONS
X *
X * The absolute error in the approximating polynomial is
X * 10**(-19.38). Note that contrary to DEC's assertion
X * in the F4P user's guide, the error is absolute and not
X * relative.
X *
X * This error bound assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes log(X) from:
X *
X * (1) If argument is zero then flag an error
X * and return minus infinity (or rather its
X * machine representation).
X *
X * (2) If argument is negative then flag an
X * error and return minus infinity.
X *
X * (3) Given that x = m * 2**k then extract
X * mantissa m and exponent k.
X *
X * (4) Transform m which is in range [0.5,1.0]
X * to range [1/sqrt(2), sqrt(2)] by:
X *
X * s = m * sqrt(2)
X *
X * (5) Compute z = (s - 1) / (s + 1)
X *
X * (6) Now use the approximation from HART
X * page 111 to find log(s):
X *
X * log(s) = z * ( P(z**2) / Q(z**2) )
X *
X * Where:
X *
X * P(z**2) = SUM [ Pj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * Q(z**2) = SUM [ Qj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = -0.240139179559210509868484e2
X * P1 = 0.30957292821537650062264e2
X * P2 = -0.96376909336868659324e1
X * P3 = 0.4210873712179797145
X * Q0 = -0.120069589779605254717525e2
X * Q1 = 0.19480966070088973051623e2
X * Q2 = -0.89111090279378312337e1
X * Q3 = 1.0000
X *
X * (coefficients from HART table #2705 pg 229)
X *
X * (7) Finally, compute log(x) from:
X *
X * log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double log_pcoeffs[] = {
X -0.24013917955921050986e2,
X 0.30957292821537650062e2,
X -0.96376909336868659324e1,
X 0.4210873712179797145
X};
X
Xstatic double log_qcoeffs[] = {
X -0.12006958977960525471e2,
X 0.19480966070088973051e2,
X -0.89111090279378312337e1,
X 1.0000
X};
X
Xstatic char funcname[] = "log";
X
X
Xdouble log (x)
Xdouble x;
X{
X auto int k;
X auto double s;
X auto double z;
X auto double zt2;
X auto double pqofz;
X auto struct exception xcpt;
X extern double frexp ();
X extern double poly ();
X
X DBUG_ENTER (funcname);
X DBUG_3 ("login", "arg %le", x);
X if (x == 0.0) {
X xcpt.type = SING;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: SINGULARITY error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = -MAXDOUBLE;
X }
X } else {
X s = SQRT2 * frexp (x, &k);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "s = %le", s);
X z = (s - 1.0) / (s + 1.0);
X DBUG_3 ("log", "z = %le", z);
X zt2 = z * z;
X DBUG_3 ("log", "zt2 = %le", zt2);
X pqofz = z * (poly (3, log_pcoeffs, zt2) / poly (3, log_qcoeffs, zt2));
X DBUG_3 ("pqofz", "pqofz = %le", pqofz);
X DBUG_3 ("log", "k = %d", k);
X DBUG_3 ("log", "LN2 = %le", LN2);
X x = k * LN2;
X DBUG_3 ("log", "x = %le", x);
X x -= LNSQRT2;
X DBUG_3 ("log", "x = %le", x);
X x += pqofz;
X DBUG_3 ("log", "x = %le", x);
X xcpt.retval = x;
X }
X DBUG_3 ("logout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/log.c
if test 4649 -ne `wc -c <funcs/src/log.c`; then
echo shar: \"funcs/src/log.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/pml.h -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/pml.h\"
else
echo shar: Extracting \"funcs/src/pml.h\" \(3819 characters\)
sed "s/^X//" >funcs/src/pml.h <<'END_OF_funcs/src/pml.h'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * This file gets included with all of the floating point math
X * library routines when they are compiled. Note that
X * this is the proper place to put machine dependencies
X * whenever possible.
X *
X * It should be pointed out that for simplicity's sake, the
X * environment parameters are defined as floating point constants,
X * rather than octal or hexadecimal initializations of allocated
X * storage areas. This means that the range of allowed numbers
X * may not exactly match the hardware's capabilities. For example,
X * if the maximum positive double precision floating point number
X * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
X * defined to be 1.11E100 then the numbers between 1.11E100 and
X * 1.11...E100 are considered to be undefined. For most
X * applications, this will cause no problems.
X *
X * An alternate method is to allocate a global static "double" variable,
X * say "maxdouble", and use a union declaration and initialization
X * to initialize it with the proper bits for the EXACT maximum value.
X * This was not done because the only compilers available to the
X * author did not fully support union initialization features.
X *
X */
X
X#ifndef NO_DBUG
X# include <dbug.h>
X#else
X# define DBUG_ENTER(a1)
X# define DBUG_RETURN(a1) return(a1)
X# define DBUG_VOID_RETURN return
X# define DBUG_EXECUTE(keyword,a1)
X# define DBUG_2(keyword,format)
X# define DBUG_3(keyword,format,a1)
X# define DBUG_4(keyword,format,a1,a2)
X# define DBUG_5(keyword,format,a1,a2,a3)
X# define DBUG_PUSH(a1)
X# define DBUG_POP()
X# define DBUG_PROCESS(a1)
X# define DBUG_FILE (stderr)
X#endif
X
X#include <values.h>
X#include <math.h>
X#include <errno.h>
X
X/*
X * MAXDOUBLE => Maximum double precision number
X * MINDOUBLE => Minimum double precision number
X * DMAXEXP => Maximum exponent of a double
X * DMINEXP => Minimum exponent of a double
X */
X
X#define MINDOUBLE (1.0/MAXDOUBLE)
X#define LOG2_MAXDOUBLE (DMAXEXP + 1)
X#define LOG2_MINDOUBLE (DMINEXP - 1)
X#define LOGE_MAXDOUBLE (LOG2_MAXDOUBLE / LOG2E)
X#define LOGE_MINDOUBLE (LOG2_MINDOUBLE / LOG2E)
X
X/*
X * The following are hacks which should be fixed when I understand all
X * the issues a little better. |tanh(TANH_MAXARG)| = 1.0
X */
X#define TANH_MAXARG 16
X#define SQRT_MAXDOUBLE 1.304380e19
X
X#define PI 3.14159265358979323846
X#define TWOPI (2.0 * PI)
X#define HALFPI (PI / 2.0)
X#define FOURTHPI (PI / 4.0)
X#define SIXTHPI (PI / 6.0)
X#define LOG2E 1.4426950408889634074 /* Log to base 2 of e */
X#define LOG10E 0.4342944819032518276
X#define SQRT2 1.4142135623730950488
X#define SQRT3 1.7320508075688772935
X#define LN2 0.6931471805599453094
X#define LNSQRT2 0.3465735902799726547
X
X
X/*
X * MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X */
X
X#ifdef IEEE
X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */
X#endif
X
X#define TRUE 1 /* This really should be in stdio.h */
X#define FALSE 0 /* This too */
X#define VOID void
END_OF_funcs/src/pml.h
if test 3819 -ne `wc -c <funcs/src/pml.h`; then
echo shar: \"funcs/src/pml.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/sin.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/sin.c\"
else
echo shar: Extracting \"funcs/src/sin.c\" \(4983 characters\)
sed "s/^X//" >funcs/src/sin.c <<'END_OF_funcs/src/sin.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sin double precision sine
X *
X * KEY WORDS
X *
X * sin
X * machine independent routines
X * trigonometric functions
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision sine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sin (x)
X * double x;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes sin(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call sin recursively
X * using relation sin(x) = -sin(x - PI).
X *
X * (3) If x < -PI/2 then call sin recursively
X * using relation sin(x) = -sin(x + PI).
X *
X * (4) If x > PI/4 then call cos using
X * relation sin(x) = cos(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation sin(x) = -cos(PI/2 + x).
X *
X * (6) If x is small enough that polynomial
X * evaluation would cause underflow
X * then return x, since sin(x)
X * approaches x as x approaches zero.
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 118 can be used:
X *
X * sin(x) = y * ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.206643433369958582409167054e+7
X * P1 = -0.18160398797407332550219213e+6
X * P2 = 0.359993069496361883172836e+4
X * P3 = -0.2010748329458861571949e+2
X * Q0 = 0.263106591026476989637710307e+7
X * Q1 = 0.3927024277464900030883986e+5
X * Q2 = 0.27811919481083844087953e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3063 pg 234)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
Xstatic double sin_pcoeffs[] = {
X 0.20664343336995858240e7,
X -0.18160398797407332550e6,
X 0.35999306949636188317e4,
X -0.20107483294588615719e2
X};
X
Xstatic double sin_qcoeffs[] = {
X 0.26310659102647698963e7,
X 0.39270242774649000308e5,
X 0.27811919481083844087e3,
X 1.0
X};
X
Xstatic char funcname[] = "sin";
X
X
Xdouble sin (x)
Xdouble x;
X{
X double y;
X double yt2;
X double rtnval;
X extern double mod ();
X extern double cos ();
X extern double poly ();
X auto struct exception xcpt;
X
X DBUG_ENTER (funcname);
X DBUG_3 ("sinin", "arg %le", x);
X if (x < -PI || x > PI) {
X x = mod (x, TWOPI);
X if (x > PI) {
X x = x - TWOPI;
X } else if (x < -PI) {
X x = x + TWOPI;
X }
X }
X if (x > HALFPI) {
X xcpt.retval = -(sin (x - PI));
X } else if (x < -HALFPI) {
X xcpt.retval = -(sin (x + PI));
X } else if (x > FOURTHPI) {
X xcpt.retval = cos (HALFPI - x);
X } else if (x < -FOURTHPI) {
X xcpt.retval = -(cos (HALFPI + x));
X } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
X xcpt.retval = x;
X } else {
X y = x / FOURTHPI;
X yt2 = y * y;
X xcpt.retval = y * (poly (3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2));
X }
X DBUG_3 ("sinout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/sin.c
if test 4983 -ne `wc -c <funcs/src/sin.c`; then
echo shar: \"funcs/src/sin.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/src/sqrt.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/src/sqrt.c\"
else
echo shar: Extracting \"funcs/src/sqrt.c\" \(4591 characters\)
sed "s/^X//" >funcs/src/sqrt.c <<'END_OF_funcs/src/sqrt.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X
X/*
X * FUNCTION
X *
X * sqrt double precision square root
X *
X * KEY WORDS
X *
X * sqrt
X * machine independent routines
X * math libraries
X *
X * DESCRIPTION
X *
X * Returns double precision square root of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sqrt (x)
X * double x;
X *
X * REFERENCES
X *
X * Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 89-96.
X *
X * RESTRICTIONS
X *
X * The relative error is 10**(-30.1) after three applications
X * of Heron's iteration for the square root.
X *
X * However, this assumes exact arithmetic in the iterations
X * and initial approximation. Additional errors may occur
X * due to truncation, rounding, or machine precision limits.
X *
X * PROGRAMMER
X *
X * Fred Fish
X *
X * INTERNALS
X *
X * Computes square root by:
X *
X * (1) Range reduction of argument to [0.5,1.0]
X * by application of identity:
X *
X * sqrt(x) = 2**(k/2) * sqrt(x * 2**(-k))
X *
X * k is the exponent when x is written as
X * a mantissa times a power of 2 (m * 2**k).
X * It is assumed that the mantissa is
X * already normalized (0.5 =< m < 1.0).
X *
X * (2) An approximation to sqrt(m) is obtained
X * from:
X *
X * u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
X *
X * P0 = 0.594604482
X * P1 = 2.54164041
X * Q0 = 2.13725758
X * Q1 = 1.0
X *
X * (coefficients from HART table #350 pg 193)
X *
X * (3) Three applications of Heron's iteration are
X * performed using:
X *
X * y[n+1] = 0.5 * (y[n] + (m/y[n]))
X *
X * where y[0] = u = approx sqrt(m)
X *
X * (4) If the value of k was odd then y is either
X * multiplied by the square root of two or
X * divided by the square root of two for k positive
X * or negative respectively. This rescales y
X * by multiplying by 2**frac(k/2).
X *
X * (5) Finally, y is rescaled by int(k/2) which
X * is equivalent to multiplication by 2**int(k/2).
X *
X * The result of steps 4 and 5 is that the value
X * of y between 0.5 and 1.0 has been rescaled by
X * 2**(k/2) which removes the original rescaling
X * done prior to finding the mantissa square root.
X *
X * NOTES
X *
X * The Convergent Technologies compiler optimizes division
X * by powers of two to "arithmetic shift right" instructions.
X * This is ok for positive integers but gives different
X * results than other C compilers for negative integers.
X * For example, (-1)/2 is -1 on the CT box, 0 on every other
X * machine I tried.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#define P0 0.594604482 /* Approximation coeff */
X#define P1 2.54164041 /* Approximation coeff */
X#define Q0 2.13725758 /* Approximation coeff */
X#define Q1 1.0 /* Approximation coeff */
X
X#define ITERATIONS 3 /* Number of iterations */
X
Xstatic char funcname[] = "sqrt";
X
X
Xdouble sqrt (x)
Xdouble x;
X{
X auto int k;
X register int bugfix;
X register int kmod2;
X register int count;
X auto int exponent;
X auto double m;
X auto double u;
X auto double y;
X auto double rtnval;
X auto struct exception xcpt;
X extern double frexp ();
X extern double ldexp ();
X
X DBUG_ENTER ("sqrt");
X DBUG_3 ("sqrtin", "arg %le", x);
X if (x == 0.0) {
X rtnval = 0.0;
X } else if (x < 0.0) {
X xcpt.type = DOMAIN;
X xcpt.name = funcname;
X xcpt.arg1 = x;
X if (!matherr (&xcpt)) {
X fprintf (stderr, "%s: DOMAIN error\n", funcname);
X errno = EDOM;
X xcpt.retval = 0.0;
X }
X } else {
X m = frexp (x, &k);
X u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
X for (count = 0, y = u; count < ITERATIONS; count++) {
X y = 0.5 * (y + (m / y));
X }
X if ((kmod2 = (k % 2)) < 0) {
X y /= SQRT2;
X } else if (kmod2 > 0) {
X y *= SQRT2;
X }
X bugfix = 2;
X xcpt.retval = ldexp (y, k/bugfix);
X }
X DBUG_3 ("sqrtout", "result %le", xcpt.retval);
X DBUG_RETURN (xcpt.retval);
X}
END_OF_funcs/src/sqrt.c
if test 4591 -ne `wc -c <funcs/src/sqrt.c`; then
echo shar: \"funcs/src/sqrt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/scale.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/scale.c\"
else
echo shar: Extracting \"funcs/unused/scale.c\" \(4092 characters\)
sed "s/^X//" >funcs/unused/scale.c <<'END_OF_funcs/unused/scale.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * scale scale a double precision number by power of 2
X *
X * KEY WORDS
X *
X * scale
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * double scale (value, scale)
X * double value;
X * integer scale;
X *
X * DESCRIPTION
X *
X * Adds a specified integer to a double precision number's
X * exponent, effectively multiplying by a power of 2 for positive
X * scale values and divided by a power of 2 for negative
X * scale values.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl
X * Tempe, Az. 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and the 68000 the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X * Scale extracts the exponent, shifts it into the lower
X * order bits, adds the scale value to it, checks for
X * exponent overflow or underflow, shifts it back into the
X * high bits, and inserts it back into the value.
X *
X * NOTE: On the DECSYSTEM-20, the exponent for negative
X * numbers is complemented. This is not true for the PDP-11.
X *
X * NOTE: On the 68000, the mantissa is the same for both
X * positive and negative numbers, only the sign is different.
X * On the PDP-11, the mantissa is two's complement. Both
X * use hidden bit normalization.
X *
X */
X
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#ifdef PDP10
X#define EXP_MASK 0377000000000 /* Mask for exponent */
X#define MANT_MASK 0400777777777 /* Mask for mantissa */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#define EXP_SHIFTS 27 /* Shifts for exp in LSBs */
X#endif
X
X#ifdef pdp11
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#endif
X
X#ifdef mc68000
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define MANT_MASK 0x807FFFFF /* Mask for mantissa */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define LEXP_MASK 0377 /* Mask for shifted exponent */
X#endif
X
Xstatic union {
X double dval;
X long lval[2];
X} share;
X
X
X
Xdouble scale(value,scale)
Xdouble value;
Xregister int scale;
X{
X register long temp1, temp2, *lpntr;
X
X lpntr = &share.lval[0];
X share.dval = value;
X temp1 = *lpntr;
X temp2 = (temp1 & EXP_MASK) >> EXP_SHIFTS;
X#ifdef PDP10
X if (value < 0.0) {
X temp2 = ~temp2 & LEXP_MASK;
X }
X#endif
X temp2 += scale;
X if (temp2 > MAX_EXPONENT+128) {
X pmlerr(SCALE_OVERFLOW);
X } else if (temp2 < 0) {
X pmlerr(SCALE_UNDERFLOW);
X } else {
X#ifdef PDP10
X if (value < 0.0) {
X temp2 = ~temp2 & LEXP_MASK;
X }
X#endif
X temp1 &= MANT_MASK;
X temp2 = temp2 << EXP_SHIFTS;
X *lpntr = temp1 | temp2;
X }
X return (share.dval);
X}
END_OF_funcs/unused/scale.c
if test 4092 -ne `wc -c <funcs/unused/scale.c`; then
echo shar: \"funcs/unused/scale.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/xexp.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/xexp.c\"
else
echo shar: Extracting \"funcs/unused/xexp.c\" \(3294 characters\)
sed "s/^X//" >funcs/unused/xexp.c <<'END_OF_funcs/unused/xexp.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * xexp extract double precision number's exponent
X *
X * KEY WORDS
X *
X * xexp
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * int xexp(value)
X * double value;
X *
X * DESCRIPTION
X *
X * Extracts exponent from a double precision number and
X * returns it as a normal length integer.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl.
X * Tempe, Az 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and mc68000 (non IEEE format) the double
X * precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * For the mc68000 using IEEE format the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEEEEMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X
X#ifdef PDP10
X#define EXP_MASK 0377000000000 /* Mask for exponent */
X#define EXP_SHIFTS 27 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X
X#ifdef mc68000
X#ifdef IEEE
X#define EXP_MASK 0x7FF00000 /* Mask for exponent */
X#define EXP_SHIFTS 20 /* Shifts to get into LSB's */
X#define EXP_BIAS 1023 /* Exponent bias */
X#else
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X#endif
X
X#ifdef pdp11
X#define EXP_MASK 0x7F800000 /* Mask for exponent */
X#define EXP_SHIFTS 23 /* Shifts to get into LSB's */
X#define EXP_BIAS 128 /* Exponent bias */
X#endif
X
Xunion dtol {
X double dval;
X int ival[2];
X};
X
X
Xint xexp (value)
Xunion dtol value;
X{
X register int *ipntr;
X
X if (value.dval == 0.0) {
X return (0);
X } else {
X ipntr = &value.ival[0];
X#ifdef PDP10
X if (value.dval < 0.0) { /* Exponent is complemented */
X *ipntr = ~*ipntr;
X }
X#endif
X *ipntr &= EXP_MASK;
X *ipntr >>= EXP_SHIFTS;
X *ipntr -= EXP_BIAS;
X return (*ipntr);
X }
X}
END_OF_funcs/unused/xexp.c
if test 3294 -ne `wc -c <funcs/unused/xexp.c`; then
echo shar: \"funcs/unused/xexp.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f funcs/unused/xmant.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"funcs/unused/xmant.c\"
else
echo shar: Extracting \"funcs/unused/xmant.c\" \(3950 characters\)
sed "s/^X//" >funcs/unused/xmant.c <<'END_OF_funcs/unused/xmant.c'
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/*
X * FUNCTION
X *
X * xmant extract double precision number's mantissa
X *
X * KEY WORDS
X *
X * xmant
X * math libraries
X * machine dependent routines
X *
X * SYNOPSIS
X *
X * double xmant (value)
X * double value;
X *
X * DESCRIPTION
X *
X * Extracts a double precision number's mantissa and returns it
X * as a double precision number with a normalized mantissa and
X * a zero exponent.
X *
X * AUTHOR
X *
X * Fred Fish
X * 1346 W 10th Pl
X * Tempe, Az 85281
X *
X * INTERNALS
X *
X * This routine is highly machine dependent. As such, no
X * attempt was made to make it general, hence it may have
X * to be completely rewritten when transportation of the
X * floating point library is attempted.
X *
X * For the DECSYSTEM-20 the double precision word format is:
X *
X * WORD N => SEEEEEEEEMMMMMMMMMMMMMMMMMMMMMMMMMMM
X * WORD N+1 => XMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
X *
X * For the PDP-11 and the mc68000 the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEMMMMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * For the mc68000 using IEEE format the double precision word
X * format is:
X *
X * WORD N => SEEEEEEEEEEEMMMM
X * WORD N+1 => MMMMMMMMMMMMMMMM
X * WORD N+2 => MMMMMMMMMMMMMMMM
X * WORD N+3 => MMMMMMMMMMMMMMMM
X *
X * Where: S => Sign bit
X * E => Exponent
X * X => Ignored (set to 0)
X * M => Mantissa bit
X *
X * NOTE: Beware of 0.0; on some machines which use excess 128
X * notation for the exponent, if the mantissa is zero the exponent
X * is also.
X *
X */
X
X#include <stdio.h>
X#include <pmluser.h>
X#include "pml.h"
X
X#ifdef PDP10
X#define MANT_MASK 0400777777777 /* Mantissa extraction mask */
X#define ZPOS_MASK 0200000000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0177000000000 /* Negative # mask for exp = 0 */
X#endif
X
X#ifdef pdp11
X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */
X#endif
X
X#ifdef mc68000
X#ifdef IEEE
X#define MANT_MASK 0x800FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x3FF00000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x3FF00000 /* Negative # mask for exp = 0 */
X#else
X#define MANT_MASK 0x807FFFFF /* Mantissa extraction mask */
X#define ZPOS_MASK 0x40000000 /* Positive # mask for exp = 0 */
X#define ZNEG_MASK 0x40000000 /* Negative # mask for exp = 0 */
X#endif
X#endif
X
Xunion dtol {
X double dval;
X int ival[2];
X};
X
X
Xdouble xmant (value)
Xunion dtol value;
X{
X register int *ipntr;
X
X ipntr = &value.ival[0];
X *ipntr &= MANT_MASK;
X *ipntr |= ZPOS_MASK;
X return (value.dval);
X}
X
X
X/****************** OLD WAY ***********
Xdouble xmant (value)
Xunion dtol value;
X{
X register long *ipntr;
X register int neg;
X
X if (value == 0.0) {
X return (value);
X } else {
X if (value < 0.0) {
X neg = TRUE;
X } else {
X neg = FALSE;
X }
X ipntr = &value.ival[0];
X *ipntr &= MANT_MASK;
X if (neg) {
X *ipntr |= ZNEG_MASK;
X } else {
X *ipntr |= ZPOS_MASK;
X }
X return (value.dval);
X }
X}
X**********************/
END_OF_funcs/unused/xmant.c
if test 3950 -ne `wc -c <funcs/unused/xmant.c`; then
echo shar: \"funcs/unused/xmant.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/c2c.dat -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/c2c.dat\"
else
echo shar: Extracting \"tests/c2c.dat\" \(5354 characters\)
sed "s/^X//" >tests/c2c.dat <<'END_OF_tests/c2c.dat'
Xcsin 1.0e00 2.0e00 3.165778547525405883e00 1.959601044654846191e00
Xcsin 0.0e00 2.0e00 0.0e00 3.626860409975051879e00
Xcsin 1.0e00 0.0e00 8.414709866046905517e-01 0.0e00
Xcsin -2.0e00 2.0e00 -3.420954853296279907e00 -1.509306490421295166e00
Xcsin -1.0e00 -2.0e00 -3.165778547525405883e00 -1.959601044654846191e00
Xcsin 1.0e00 -2.0e00 3.165778547525405883e00 -1.959601044654846191e00
Xccos 1.0e00 2.0e00 2.032723009586334228e00 -3.051897794008255004e00
Xccos 0.0e00 2.0e00 3.762195706367492675e00 0.0e00
Xccos 1.0e00 0.0e00 5.403023064136505127e-01 0.0e00
Xccos -2.0e00 2.0e00 -1.565625846385955810e00 3.297894805669784545e00
Xccos -1.0e00 -2.0e00 2.032723009586334228e00 -3.051897794008255004e00
Xccos 1.0e00 -2.0e00 2.032723009586334228e00 3.051897794008255004e00
Xcexp 1.0e00 2.0e00 -1.131204381585121154e00 2.471726655960083007e00
Xcexp 0.0e00 2.0e00 -4.161468371748924255e-01 9.092974215745925903e-01
Xcexp 1.0e00 0.0e00 2.718281835317611694e00 0.0e00
Xcexp -2.0e00 2.0e00 -5.631934991106390953e-02 1.230600243434309959e-01
Xcexp -1.0e00 -2.0e00 -1.530918665230274200e-01 -3.345118276774883270e-01
Xcexp 1.0e00 -2.0e00 -1.131204381585121154e00 -2.471726655960083007e00
Xctanh 1.0e00 2.0e00 1.166736245155334472e00 -2.434581927955150604e-01
Xctanh 0.0e00 2.0e00 0.0e00 -2.185039848089218139e00
Xctanh 1.0e00 0.0e00 7.615941539406776428e-01 0.0e00
Xctanh -2.0e00 2.0e00 -1.023835599422454834e00 -2.839295566082000732e-02
Xctanh -1.0e00 -2.0e00 -1.166736245155334472e00 2.434581927955150604e-01
Xctanh 1.0e00 -2.0e00 1.166736245155334472e00 2.434581927955150604e-01
Xctan 1.0e00 2.0e00 3.381283208727836608e-02 1.014793619513511657e00
Xctan 0.0e00 2.0e00 0.0e00 9.640275761485099792e-01
Xctan 1.0e00 0.0e00 1.557407721877098083e00 0.0e00
Xctan -2.0e00 2.0e00 2.839295566082000732e-02 1.023835599422454834e00
Xctan -1.0e00 -2.0e00 -3.381283208727836608e-02 -1.014793619513511657e00
Xctan 1.0e00 -2.0e00 3.381283208727836608e-02 -1.014793619513511657e00
Xcsqrt 0.0 0.0 0.000000000000000000e00 0.000000000000000000e00
Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01
Xcsqrt 0.0 2.0 1.000000000000000000e00 1.000000000000000000e00
Xcsqrt 1.0 0.0 1.000000000000000000e00 0.000000000000000000e00
Xcsqrt -2.0 2.0 6.435942575335502624e-01 1.553773969411849975e00
Xcsqrt -1.0 2.0 7.861513718962669372e-01 1.272019654512405395e00
Xcsqrt 1.0 2.0 1.272019654512405395e00 7.861513718962669372e-01
Xcsinh 1.0e00 2.0e00 -4.890562593936920166e-01 1.403119236230850219e00
Xcsinh 0.0e00 2.0e00 0.0e00 9.092974215745925903e-01
Xcsinh 1.0e00 0.0e00 1.175201192498207092e00 0.0e00
Xcsinh -2.0e00 2.0e00 1.509306475520133972e00 3.420954853296279907e00
Xcsinh -1.0e00 -2.0e00 4.890562593936920166e-01 -1.403119236230850219e00
Xcsinh 1.0e00 -2.0e00 -4.890562593936920166e-01 -1.403119236230850219e00
Xcrcp 1.0e00 2.0e00 1.999999992549419403e-01 -3.999999985098838806e-01
Xcrcp 0.0e00 2.0e00 0.0e00 -5.0e-01
Xcrcp 1.0e00 0.0e00 1.0e00 0.0e00
Xcrcp -2.0e00 2.0e00 -2.5e-01 -2.5e-01
Xcrcp -1.0e00 -2.0e00 -1.999999992549419403e-01 3.999999985098838806e-01
Xcrcp 1.0e00 -2.0e00 1.999999992549419403e-01 3.999999985098838806e-01
Xclog 1.0e00 2.0e00 8.047189563512802124e-01 1.107148721814155578e00
Xclog 0.0e00 2.0e00 6.931471824645996093e-01 1.570796325802803039e00
Xclog 1.0e00 0.0e00 0.0e00 0.0e00
Xclog -2.0e00 2.0e00 1.039720773696899414e00 2.356194496154785156e00
Xclog -1.0e00 -2.0e00 8.047189563512802124e-01 -2.034443944692611694e00
Xclog 1.0e00 -2.0e00 8.047189563512802124e-01 -1.107148721814155578e00
Xccosh 1.0e00 2.0e0 -6.421481221914291381e-01 1.068607419729232788e00
Xccosh 0.0e00 2.0e00 -4.161468371748924255e-01 0.0e00
Xccosh 1.0e00 0.0e00 1.543080642819404602e00 0.0e00
Xccosh -2.0e00 2.0e00 -1.565625831484794616e00 -3.297894805669784545e00
Xccosh -1.0e00 -2.0e00 -6.421481221914291381e-01 1.068607419729232788e00
Xccosh 1.0e00 -2.0e00 -6.421481221914291381e-01 -1.068607419729232788e00
Xcatan 1.0e00 2.0e00 1.338972523808479309e00 4.023594781756401062e-01
Xcatan 0.0e00 2.0e00 1.570796325802803039e00 5.493061468005180358e-01 /* -real part */
Xcatan 1.0e00 0.0e00 7.853981629014015197e-01 0.0e00
Xcatan -2.0e00 2.0e00 -1.311223268508911132e00 2.388778626918792724e-01
Xcatan -1.0e00 -2.0e00 -1.338972523808479309e00 -4.023594819009304046e-01
Xcatan 1.0e00 -2.0e00 1.338972523808479309e00 -4.023594819009304046e-01
Xcasin 1.0 2.0 4.270785786211490631e-01 1.528570890426635742e00
Xcasin 0.0 2.0 0.000000000000000000e00 1.443635478615760803e00
Xcasin 1.0 0.0 1.570796325802803039e00 0.000000000000000000e00
Xcasin -2.0 2.0 -7.542491331696510314e-01 1.734324455261230468e00
Xcasin -1.0 -2.0 -4.270785823464393615e-01 -1.528570920228958129e00
Xcasin 1.0 -2.0 4.270785823464393615e-01 -1.528570920228958129e00
Xcacos 1.0 2.0 1.143717750906944274e00 -1.528570920228958129e00
Xcacos 0.0 2.0 1.570796325802803039e00 -1.443635478615760803e00
Xcacos 1.0 0.0 0.000000000000000000e00 0.000000000000000000e00
Xcacos -2.0 2.0 2.325045466423034668e00 -1.734324529767036438e00
Xcacos -1.0 -2.0 1.997874900698661804e00 1.528570890426635742e00
Xcacos 1.0 -2.0 1.143717750906944274e00 1.528570890426635742e00
END_OF_tests/c2c.dat
if test 5354 -ne `wc -c <tests/c2c.dat`; then
echo shar: \"tests/c2c.dat\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f tests/cc2c.dat -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"tests/cc2c.dat\"
else
echo shar: Extracting \"tests/cc2c.dat\" \(3825 characters\)
sed "s/^X//" >tests/cc2c.dat <<'END_OF_tests/cc2c.dat'
Xcadd 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.683000028133392334e00 6.440000057220458984e00
Xcadd -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 2.003005027770996093e00 -1.887652993202209472e00
Xcadd 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 3.719695717096328735e00 4.449983000755310058e00
Xcadd 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.954332172870635986e00 0.000000000000000000e00
Xcadd -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -4.234000027179718017e00 -6.334455549716949462e00
Xcadd 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.300034515380859375e03
Xcdiv 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 4.609979800879955291e-01 1.263787318021059036e-01
Xcdiv -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.863302707672119140e-01 6.862613558769226074e-02
Xcdiv 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 5.729839205741882324e-02 -4.476501746103167533e-02
Xcdiv 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 4.663971476256847381e-01 0.000000000000000000e00
Xcdiv -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 5.215928852558135986e-01 8.269466646015644073e-02
Xcdiv 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 3.765932763671875000e04 -3.554461171875000000e05
Xcmult 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -5.596120119094848632e00 1.293470001220703125e01
Xcmult -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 4.705165147781372070e00 1.302998638153076171e01
Xcmult 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 8.510770574212074279e-01 2.011436134576797485e00
Xcmult 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 5.323814511299133300e00 0.000000000000000000e00
Xcmult -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 -5.635822236537933349e00 1.193936669826507568e01
Xcmult 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 -4.487600040435791015e01 4.235604019165039062e02
Xcsubt 1.122999995946884155e00 2.340000003576278686e00 3.560000002384185791e00 4.100000023841857910e00 -2.437000006437301635e00 -1.760000020265579223e00
Xcsubt -1.345669999718666076e00 2.234566986560821533e00 3.348675012588500976e00 -4.122219979763031005e00 -4.694344997406005859e00 6.356786966323852539e00
Xcsubt 3.857562988996505737e-01 1.000030040740966796e-01 3.333939403295516967e00 4.349979996681213378e00 -2.948183119297027587e00 -4.249976992607116699e00
Xcsubt 1.575757563114166259e00 0.000000000000000000e00 3.378574609756469726e00 0.000000000000000000e00 -1.802817046642303466e00 0.000000000000000000e00
Xcsubt -1.233999997377395629e00 -2.334455549716949462e00 -3.000000000000000000e00 -4.000000000000000000e00 1.766000002622604370e00 1.665544450283050537e00
Xcsubt 1.227000000000000000e04 1.300000000000000000e03 0.000000000000000000e00 3.452000021934509277e-02 1.227000000000000000e04 1.299965484619140625e03
END_OF_tests/cc2c.dat
if test 3825 -ne `wc -c <tests/cc2c.dat`; then
echo shar: \"tests/cc2c.dat\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 3 \(of 6\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
= Drug tests; just say *NO*! (Moto just announced new drug testing program) =
= Fred Fish Motorola Computer Division, 3013 S 52nd St, Tempe, Az 85282 USA =
= seismo!noao!mcdsun!fnf (602) 438-5976 =