Path: utzoo!utgpu!jarvis.csri.toronto.edu!clyde.concordia.ca!uunet!munnari.oz.au!cluster!metro!sunaus.oz!richb
From: richb@sunaus.oz (Rich Burridge)
Newsgroups: comp.sources.bugs
Subject: Official patch #6 for calctool v2.4 (Part 2 of 3).
Message-ID: <1990Feb2.022139.8453@sunaus.oz>
Date: 2 Feb 90 02:21:39 GMT
Organization: Sun Microsystems - Australia
Lines: 1347

---- 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