Official patch #6 for calctool v2.4 (Part 2 of 3).
Rich Burridge
richb at sunaus.oz
Fri Feb 2 13:21:39 AEST 1990
---- Cut Here and unpack ----
#!/bin/sh
# this is part 2 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file mathlib.h continued
#
CurArch=2
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file mathlib.h"
sed 's/^X//' << 'SHAR_EOF' >> mathlib.h
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
Xextern double acos(), acosh(), asin(), asinh(), atan(), atanh() ;
Xextern double cos(), cosh(), exp(), fabs(), log(), log10(), pow() ;
Xextern double sin(), sinh(), sqrt(), tan(), tanh() ;
X
X/*START============start of definitions from <values.h>============START
X *
X * If your system has a /usr/include/values.h, or has another include
X * file which defines:
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 * you can comment out these definitions down to the END line below.
X */
X
X#ifndef BITSPERBYTE
X/* These values work with any binary representation of integers
X * where the high-order bit contains the sign. */
X
X/* a number used normally for size of a shift */
X#if gcos
X#define BITSPERBYTE 9
X#else
X#define BITSPERBYTE 8
X#endif
X#define BITS(type) (BITSPERBYTE * (int)sizeof(type))
X
X/* Various values that describe the binary floating-point representation
X * MAXDOUBLE - the largest double
X * ((_EXPBASE ** DMAXEXP) * (1 - (_EXPBASE ** -DSIGNIF)))
X * MINDOUBLE - the smallest double (_EXPBASE ** (DMINEXP - 1))
X * DMAXEXP - the maximum exponent of a double (as returned by frexp())
X * DMINEXP - the minimum exponent of a double (as returned by frexp())
X * DSIGNIF - the number of significant bits in a double
X * _IEEE - 1 if IEEE standard representation is used
X * _DEXPLEN - the number of bits for the exponent of a double
X * _HIDDENBIT - 1 if high-significance bit of mantissa is implicit
X */
X
X#if u3b || u3b5 || sun
X#define MAXDOUBLE 1.79769313486231470e+308
X#define MINDOUBLE 4.94065645841246544e-324
X#define _IEEE 1
X#define _DEXPLEN 11
X#define _HIDDENBIT 1
X#define DMINEXP (-(DMAXEXP + DSIGNIF - _HIDDENBIT - 3))
X#endif
X
X#if pdp11 || vax
X#define MAXDOUBLE 1.701411834604692293e+38
X#define MINDOUBLE (0.01 * 2.938735877055718770e-37)
X#define _IEEE 0
X#define _DEXPLEN 8
X#define _HIDDENBIT 1
X#define DMINEXP (-DMAXEXP)
X#endif
X
X#if gcos
X#define MAXDOUBLE 1.7014118346046923171e+38
X#define MINDOUBLE 2.9387358770557187699e-39
X#define _IEEE 0
X#define _DEXPLEN 8
X#define _HIDDENBIT 0
X#define DMINEXP (-(DMAXEXP + 1))
X#endif
X
X#define DSIGNIF (BITS(double) - _DEXPLEN + _HIDDENBIT - 1)
X#define DMAXEXP ((1 << _DEXPLEN - 1) - 1 + _IEEE)
X
X#endif /*BITSPERBYTE*/
X
X/*END==============end of definitions from <values.h>==================END*/
X
X
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
X#define TANH_MAXARG 16
X#define SQRT_MAXDOUBLE 1.304380e19
X
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/* MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X * Apologies for the double negative. I want as few definitions
X * needed in the default case as possible.
X */
X
X#ifndef NOIEEE
X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */
X#endif /*NOIEEE*/
X
X/* It is hoped that your system supplies all the mathematical functions
X * required by calctool. If not then, it is possible to use the needed
X * ones from the portable math library routines that comes with these
X * sources.
X *
X * There is one definition for each routine used by calctool. These are
X * commented out by default to signify that this system has that routine.
X * If you are missing any, then uncomment the appropriate definitions.
X */
X
X/*#define NEED_ACOS */
X/*#define NEED_ACOSH */
X/*#define NEED_ASIN */
X/*#define NEED_ASINH */
X/*#define NEED_ATAN */
X/*#define NEED_ATANH */
X/*#define NEED_COS */
X/*#define NEED_COSH */
X/*#define NEED_EXP */
X/*#define NEED_LOG */
X/*#define NEED_LOG10 */
X/*#define NEED_POW */
X/*#define NEED_SIN */
X/*#define NEED_SINH */
X/*#define NEED_SQRT */
X/*#define NEED_TAN */
X/*#define NEED_TANH */
SHAR_EOF
echo "File mathlib.h is complete"
chmod 0444 mathlib.h || echo "restore of mathlib.h fails"
set `wc -c mathlib.h`;Sum=$1
if test "$Sum" != "8385"
then echo original size 8385, current size $Sum;fi
echo "x - extracting mathlib.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > mathlib.c &&
X
X/* @(#)mathlib.c 1.2 90/02/02
X *
X * These are the mathematical routines used by calctool.
X *
X * This is being done because libm.a doesn't appear to be as portable
X * as originally assumed.
X *
X * It is hoped that your system supplies all the mathematical functions
X * required by calctool. If not then, it is possible to use the needed
X * ones from this library. Look in mathlib.h for a set of definitions,
X * and uncomment and set appropriately.
X *
X * These routines are taken from two sources:
X *
X * 1/ Fred Fishs' portable maths library which was posted to the
X * comp.sources.unix newsgroup on April 1987.
X *
X * acos, acosh, asin, asinh, atan, atanh, cos, cosh, dabs,
X * exp, log, log10, mod, poly, sin, sinh, sqrt, tan, tanh.
X *
X * 2/ The BSD4.3 maths library found on uunet in
X * bsd_sources/src/usr.lib/libm.
X *
X * pow, pow_p, scalb, logb, copysign, finite, drem, exp__E,
X * log__L.
X *
X * Customised and condensed by Rich Burridge,
X * Sun Microsystems, Australia.
X *
X * Permission is given to distribute these sources, as long as the
X * copyright messages are not removed, and no monies are exchanged.
X *
X * No responsibility is taken for any errors or inaccuracies inherent
X * either to the comments or the code of this program, but if
X * reported to me then an attempt will be made to fix them.
X */
X
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#include <stdio.h>
X#include <errno.h>
X#include "mathlib.h"
X#include "calctool.h"
X
Xdouble atan(), cos(), cosh(), dabs(), exp(), frexp() ;
Xdouble ldexp(), log(), mod(), modf(), poly(), sin() ;
Xdouble sinh(), sqrt() ;
X
Xextern double frexp(), ldexp(), modf() ;
X
X
X/* FUNCTION
X *
X * acos double precision arc cosine
X *
X * DESCRIPTION
X *
X * Returns double precision arc cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double acos(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-1.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * in atan is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes arccosine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then call
X * doerr and return 0.0 by default.
X *
X * (2) If x = 0.0 then acos(x) = PI/2.
X *
X * (3) If x = 1.0 then acos(x) = 0.0
X *
X * (4) If x = -1.0 then acos(x) = PI.
X *
X * (5) If 0.0 < x < 1.0 then
X * acos(x) = atan(Y)
X * Y = sqrt[1-(x**2)] / x
X *
X * (6) If -1.0 < x < 0.0 then
X * acos(x) = atan(Y) + PI
X * Y = sqrt[1-(x**2)] / x
X */
X
X#ifdef NEED_ACOS
Xdouble
Xacos(x)
Xdouble x ;
X{
X double y ;
X auto double retval ;
X
X if (x > 1.0 || x < -1.0)
X {
X doerr("acos", "DOMAIN", EDOM) ;
X retval = 0.0 ;
X }
X else if (x == 0.0) retval = HALFPI ;
X else if (x == 1.0) retval = 0.0 ;
X else if (x == -1.0) retval = PI ;
X else
X {
X y = atan(sqrt(1.0 - (x * x)) / x) ;
X if (x > 0.0) retval = y ;
X else retval = y + PI ;
X }
X return(retval) ;
X}
X#endif /*NEED_ACOS*/
X
X
X/* FUNCTION
X *
X * acosh double precision hyperbolic arc cosine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double acosh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The range of the ACOSH function is all real numbers greater
X * than or equal to 1.0 however large arguments may cause
X * overflow in the x squared portion of the function evaluation.
X *
X * For precision information refer to documentation of the
X * floating point library primatives called.
X *
X * INTERNALS
X *
X * Computes acosh(x) from:
X *
X * 1. If x < 1.0 then report illegal
X * argument and return zero.
X *
X * 2. If x > sqrt(MAXDOUBLE) then
X * set x = sqrt(MAXDOUBLE and
X * continue after reporting overflow.
X *
X * 3. acosh(x) = log [x+sqrt(x**2 - 1)]
X */
X
X#ifdef NEED_ACOSH
Xdouble
Xacosh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x < 1.0)
X {
X doerr("acosh", "DOMAIN", ERANGE) ;
X retval = 0.0 ;
X }
X else if (x > SQRT_MAXDOUBLE)
X {
X doerr("acosh", "OVERFLOW", ERANGE) ;
X x = SQRT_MAXDOUBLE ;
X retval = log(2* SQRT_MAXDOUBLE) ;
X }
X else retval = log(x + sqrt(x*x - 1.0)) ;
X return(retval) ;
X}
X#endif /*NEED_ACOSH*/
X
X
X/*
X * FUNCTION
X *
X * asin double precision arc sine
X *
X * DESCRIPTION
X *
X * Returns double precision arc sine of double precision
X * floating point argument.
X *
X * If argument is less than -1.0 or greater than +1.0, calls
X * doerr with a DOMAIN error. If doerr does not handle
X * the error then prints error message and returns 0.
X *
X * USAGE
X *
X * double asin(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-2.
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the floating
X * point library primatives called.
X *
X * INTERNALS
X *
X * Computes arcsine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then
X * call doerr and return 0.0 by default.
X *
X * (2) If x = 0.0 then asin(x) = 0.0
X *
X * (3) If x = 1.0 then asin(x) = PI/2.
X *
X * (4) If x = -1.0 then asin(x) = -PI/2
X *
X * (5) If -1.0 < x < 1.0 then
X * asin(x) = atan(y)
X * y = x / sqrt[1-(x**2)]
X */
X
X#ifdef NEED_ASIN
Xdouble
Xasin(x)
Xdouble x ;
X{
X auto double retval ;
X
X if ( x > 1.0 || x < -1.0)
X {
X doerr("asin", "DOMAIN", EDOM) ;
X retval = 0.0 ;
X }
X else if (x == 0.0) retval = 0.0 ;
X else if (x == 1.0) retval = HALFPI ;
X else if (x == -1.0) retval = -HALFPI ;
X else retval = atan(x / sqrt (1.0 - (x * x))) ;
X return(retval) ;
X}
X#endif /*NEED_ASIN*/
X
X
X/* FUNCTION
X *
X * asinh double precision hyperbolic arc sine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc sine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double asinh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The domain of the ASINH function is the entire real axis
X * however the evaluation of x squared may cause overflow
X * for large magnitudes.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes asinh(x) from:
X *
X * 1. Let xmax = sqrt(MAXDOUBLE - 1)
X *
X * 2. If x < -xmax or xmax < x then
X * let x = xmax and flag overflow.
X *
X * 3. asinh(x) = log [x+sqrt(x**2 + 1)]
X */
X
X#ifdef NEED_ASINH
Xdouble
Xasinh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x < -SQRT_MAXDOUBLE || x > SQRT_MAXDOUBLE)
X {
X doerr("asinh", "OVERFLOW", ERANGE) ;
X retval = log(2 * SQRT_MAXDOUBLE) ;
X }
X else retval = log(x + sqrt(x*x + 1.0)) ;
X return(retval) ;
X}
X#endif /*NEED_ASINH*/
X
X
X/* FUNCTION
X *
X * atan double precision arc tangent
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 * 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#ifdef NEED_ATAN
X
X#define LAST_BOUND 0.2679491924311227074725 /* tan (PI/12) */
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
Xdouble
Xatan(x)
Xdouble x ;
X{
X register int order ;
X double t1, t2, xt2 ;
X auto double retval ;
X
X if (x < 0.0) retval = -(atan(-x)) ;
X else if (x > 1.0)
X {
X if (x < MAXDOUBLE && x > -MAXDOUBLE)
X retval = HALFPI - atan(1.0 / x) ;
X else
X {
X doerr("atan", "UNDERFLOW", EDOM) ;
X retval = 0.0 ;
X }
X }
X else if (x > LAST_BOUND)
X {
X t1 = x * SQRT3 - 1.0 ;
X t2 = SQRT3 + x ;
X retval = SIXTHPI + atan(t1 / t2) ;
X }
X else if (x < X16_UNDERFLOWS)
X {
X doerr("atan", "PLOSS", EDOM) ;
X retval = x ;
X }
X else
X {
X xt2 = x * x ;
X order = sizeof(atan_coeffs) / sizeof(double) ;
X order -= 1 ;
X retval = x * poly(order, atan_coeffs, xt2) ;
X }
X return(retval) ;
X}
X#endif /*NEED_ATAN*/
X
X
X/* FUNCTION
X *
X * atanh double precision hyperbolic arc tangent
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc tangent of double precision
X * floating point number.
X *
X * USAGE
X *
X * double atanh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The range of the atanh function is -1.0 to +1.0 exclusive.
X * Certain pathological cases near 1.0 may cause overflow
X * in evaluation of 1+x / 1-x, depending upon machine exponent
X * range and mantissa precision.
X *
X * For precision information refer to documentation of the
X * other floating point library routines called.
X *
X * INTERNALS
X *
X * Computes atanh(x) from:
X *
X * 1. If x <= -1.0 or x >= 1.0
X * then report argument out of range and return 0.0
X *
X * 2. atanh(x) = 0.5 * log((1+x)/(1-x))
X */
X
X#ifdef NEED_ATANH
Xdouble
Xatanh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x <= -1.0 || x >= 1.0)
X {
X doerr("atan", "DOMAIN", ERANGE) ;
X retval = 0.0 ;
X }
X else retval = 0.5 * log((1+x) / (1-x)) ;
X return(retval) ;
X}
X#endif /*NEED_ATANH*/
X
X
X/* FUNCTION
X *
X * cos - double precision cosine
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 * 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#ifdef NEED_COS
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
Xdouble
Xcos(x)
Xdouble x ;
X{
X auto double y, yt2 ;
X auto double retval ;
X
X if (x < -PI || x > PI)
X {
X x = mod(x, TWOPI) ;
X if (x > PI) x = x - TWOPI ;
X else if (x < -PI) x = x + TWOPI ;
X }
X if (x > HALFPI) retval = -(cos(x - PI)) ;
X else if (x < -HALFPI) retval = -(cos(x + PI)) ;
X else if (x > FOURTHPI) retval = sin(HALFPI - x) ;
X else if (x < -FOURTHPI) retval = sin(HALFPI + x) ;
X else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS)
X retval = sqrt(1.0 - (x * x)) ;
X else
X {
X y = x / FOURTHPI ;
X yt2 = y * y ;
X retval = poly(3, cos_pcoeffs, yt2) / poly(3, cos_qcoeffs, yt2) ;
X }
X return(retval) ;
X}
X#endif /*NEED_COS*/
X
X
X/* FUNCTION
X *
X * cosh double precision hyperbolic cosine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double cosh(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-4
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes hyperbolic cosine from:
X *
X * cosh(X) = 0.5 * (exp(X) + exp(-X))
X */
X
X
X#ifdef NEED_COSH
Xdouble
Xcosh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x > LOGE_MAXDOUBLE)
X {
X doerr("cosh", "OVERFLOW", ERANGE) ;
X retval = MAXDOUBLE ;
X }
X else if (x < LOGE_MINDOUBLE)
X {
X doerr("cosh", "UNDERFLOW", ERANGE) ;
X retval = MINDOUBLE ;
X }
X else
X {
X x = exp(x) ;
X retval = 0.5 * (x + 1.0 / x) ;
X }
X return(retval) ;
X}
X#endif /*NEED_COSH*/
X
X
X/* FUNCTION
X *
X * dabs double precision absolute value
X *
X * DESCRIPTION
X *
X * Returns absolute value of double precision number.
X *
X * USAGE
X *
X * double dabs(x)
X * double x ;
X */
X
X#ifdef NEED_EXP
Xdouble
Xdabs(x)
Xdouble x ;
X{
X if (x < 0.0) x = -x ;
X return(x) ;
X}
X#endif /*NEED_EXP*/
X
X
X/* FUNCTION
X *
X * exp double precision exponential
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 *
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#ifdef NEED_EXP
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 * 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
Xdouble
Xexp(x)
Xdouble x ;
X{
X register int ival, y ;
X auto double a, b, rtnval, t, temp, v, w, wpof2, zpof2 ;
X auto double retval ;
X
X if (x > LOGE_MAXDOUBLE)
X {
X doerr("exp", "OVERFLOW", ERANGE) ;
X retval = MAXDOUBLE ;
X }
X else if (x <= LOGE_MINDOUBLE)
X {
X doerr("exp", "UNDERFLOW", ERANGE) ;
X retval = 0.0 ;
X }
X else
X {
X t = x * LOG2E ;
X (void) modf(t, &temp) ;
X y = temp ;
X v = 16.0 * modf(t, &temp) ;
X (void) modf(dabs(v), &temp) ;
X ival = temp ;
X if (x < 0.0) zpof2 = 1.0 / fpof2tbl[ival] ;
X else zpof2 = fpof2tbl[ival] ;
X w = modf(v, &temp) / 16.0 ;
X a = C + (D * w * w) ;
X b = w * (E + (F * w * w)) ;
X wpof2 = (a + b) / (a - b) ;
X retval = ldexp((wpof2 * zpof2), y) ;
X }
X return(retval) ;
X}
X#endif /*NEED_EXP*/
X
X
X/* FUNCTION
X *
X * log double precision natural log
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 * 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#ifdef NEED_LOG
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
Xdouble
Xlog(x)
Xdouble x ;
X{
X auto int k ;
X auto double pqofz, s, z, zt2 ;
X auto double retval ;
X
X if (x == 0.0)
X {
X doerr("log", "SINGULARITY", EDOM) ;
X retval = -MAXDOUBLE ;
X }
X else if (x < 0.0)
X {
X doerr("log", "DOMAIN", EDOM) ;
X retval = -MAXDOUBLE ;
X }
X else
X {
X s = SQRT2 * frexp(x, &k) ;
X z = (s - 1.0) / (s + 1.0) ;
X zt2 = z * z ;
X pqofz = z * (poly(3, log_pcoeffs, zt2) / poly(3, log_qcoeffs, zt2)) ;
X x = k * LN2 ;
X x -= LNSQRT2 ;
X x += pqofz ;
X retval = x ;
X }
X return(retval) ;
X}
X#endif /*NEED_LOG*/
X
X
X/* FUNCTION
X *
X * log10 double precision common log
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 * INTERNALS
X *
X * Computes log10(x) from:
X *
X * log10(x) = log10(e) * log(x)
X */
X
X#ifdef NEED_LOG10
Xdouble
Xlog10(x)
Xdouble x ;
X{
X x = LOG10E * log(x) ;
X return(x) ;
X}
X#endif /*NEED_LOG10*/
X
X
X/* FUNCTION
X *
X * mod double precision modulo
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
X#ifdef NEED_COS || NEED_SIN
Xdouble mod(value, base)
Xdouble value ;
Xdouble base ;
X{
X auto double intpart ;
X
X value /= base ;
X value = modf(value, &intpart) ;
X value *= base ;
X return(value) ;
X}
X#endif /*NEED_COS || NEED_SIN*/
X
X
X/* FUNCTION
X *
X * poly double precision polynomial evaluation
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 * 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#ifdef NEED_ATAN || NEED_COS || NEED_LOG || NEED_SIN
Xdouble
Xpoly(order, coeffs, x)
Xregister int order ;
Xdouble *coeffs ;
Xdouble x ;
X{
X auto double curr_coeff, rtn_value ;
X
X if (order <= 0) rtn_value = *coeffs ;
X else
X {
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 return(rtn_value) ;
X}
X#endif /*NEED_ATAN || NEED_COS || NEED_LOG || NEED_SIN*/
SHAR_EOF
echo "End of part 2"
echo "File mathlib.c is continued in part 3"
echo "3" > s2_seq_.tmp
exit 0
More information about the Comp.sources.bugs
mailing list