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 3 of 3).
Message-ID: <1990Feb2.022225.8521@sunaus.oz>
Date: 2 Feb 90 02:22:25 GMT
Organization: Sun Microsystems - Australia
Lines: 1225

---- Cut Here and unpack ----
#!/bin/sh
# this is part 3 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file mathlib.c continued
#
CurArch=3
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.c"
sed 's/^X//' << 'SHAR_EOF' >> mathlib.c
X
X
X/*  FUNCTION
X *
X *  sin  double precision sine
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 *  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#ifdef   NEED_SIN
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
Xdouble
Xsin(x)
Xdouble x ;
X{
X  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 = -(sin(x - PI)) ;
X  else if (x < -HALFPI)   retval = -(sin(x + PI)) ;
X  else if (x > FOURTHPI)  retval = cos(HALFPI - x) ;
X  else if (x < -FOURTHPI) retval = -(cos(HALFPI + x)) ;
X  else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) retval = x ;
X  else
X    {
X      y = x / FOURTHPI ;
X      yt2 = y * y ;
X      retval = y * (poly(3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2)) ;
X    }
X  return(retval) ;
X}
X#endif /*NEED_SIN*/
X
X
X/*  FUNCTION
X *
X *      sinh   double precision hyperbolic sine
X *
X *  DESCRIPTION
X *
X *      Returns double precision hyperbolic sine of double precision
X *      floating point number.
X *
X *  USAGE
X *
X *      double sinh(x)
X *      double x ;
X *
X *  REFERENCES
X *
X *      Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X *  RESTRICTIONS
X *
X *      Inputs greater than ln(MAXDOUBLE) result in overflow.
X *      Inputs less than ln(MINDOUBLE) result in underflow.
X *
X *      For precision information refer to documentation of the
X *      floating point library routines called.
X *
X *  INTERNALS
X *
X *      Computes hyperbolic sine from:
X *
X *              sinh(x) = 0.5 * (exp(x) - 1.0/exp(x))
X */
X
X#ifdef   NEED_SINH
Xdouble
Xsinh(x)
Xdouble x ;
X{
X  auto double retval ;
X
X  if (x > LOGE_MAXDOUBLE)
X    {
X      doerr("sinh", "OVERFLOW", ERANGE) ;
X      retval = MAXDOUBLE ;
X    }
X  else if (x < LOGE_MINDOUBLE)
X    {
X      doerr("sinh", "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_SINH*/
X
X
X/*  FUNCTION
X *
X *      sqrt   double precision square root
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 *  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#ifdef   NEED_SQRT
X
X#define  ITERATIONS  3                        /* Number of iterations */
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
Xdouble
Xsqrt(x)
Xdouble x ;
X{
X  register int bugfix, count, kmod2 ;
X  auto int exponent, k ;
X  auto double m, u, y ;
X  auto double retval ;
X  
X  if (x == 0.0) retval = 0.0 ;
X  else if (x < 0.0)
X    {
X      doerr("sqrt", "DOMAIN", EDOM) ;
X      retval = 0.0 ;
X    }
X  else
X    {
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           if ((kmod2 = (k % 2)) < 0) y /= SQRT2 ;
X      else if (kmod2 > 0)             y *= SQRT2 ;
X      bugfix = 2 ;
X      retval = ldexp(y, k / bugfix) ;
X    }
X  return(retval) ;
X}
X#endif /*NEED_SQRT*/
X
X
X/*  FUNCTION
X *
X *      tan   Double precision tangent
X *
X *  DESCRIPTION
X *
X *      Returns tangent of double precision floating point number.
X *
X *  USAGE
X *
X *      double tan(x)
X *      double x ;
X *
X *  INTERNALS
X *
X *      Computes the tangent from tan(x) = sin(x) / cos(x).
X *
X *      If cos(x) = 0 and sin(x) >= 0, then returns largest
X *      floating point number on host machine.
X *
X *      If cos(x) = 0 and sin(x) < 0, then returns smallest
X *      floating point number on host machine.
X *
X *  REFERENCES
X *
X *      Fortran IV plus user's guide, Digital Equipment Corp. pp. B-8
X */
X
X#ifdef   NEED_TAN
Xdouble
Xtan(x)
Xdouble x ;
X{
X  double cosx, sinx ;
X  auto double retval ;
X
X  sinx = sin(x) ;
X  cosx = cos(x) ;
X  if (cosx == 0.0)
X    {
X      doerr("tan", "OVERFLOW", ERANGE) ;
X      if (sinx >= 0.0) retval = MAXDOUBLE ;
X      else             retval = -MAXDOUBLE ;
X    }
X  else retval = sinx / cosx ;
X  return(retval) ;
X}
X#endif /*NEED_TAN*/
X
X
X/*  FUNCTION
X *
X *  tanh   double precision hyperbolic tangent
X *
X *  DESCRIPTION
X *
X *  Returns double precision hyperbolic tangent of double precision
X *  floating point number.
X *
X *  USAGE
X *
X *  double tanh(x)
X *  double x ;
X *
X *  REFERENCES
X *
X *  Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X *  RESTRICTIONS
X *
X *  For precision information refer to documentation of the
X *  floating point library routines called.
X *  
X *  INTERNALS
X *
X *  Computes hyperbolic tangent from:
X *
X *    tanh(x) = 1.0 for x > TANH_MAXARG
X *      = -1.0 for x < -TANH_MAXARG
X *      =  sinh(x) / cosh(x) otherwise
X */
X
X#ifdef   NEED_TANH
Xdouble
Xtanh(x)
Xdouble x ;
X{
X  auto double retval ;
X  register int positive ;
X
X  if (x > TANH_MAXARG || x < -TANH_MAXARG)
X    {
X      if (x > 0.0) positive = 1 ;
X      else positive = 0 ;
X      doerr("tanh", "PLOSS", ERANGE) ;
X      if (positive) retval = 1.0 ;
X      else          retval = -1.0 ;
X    }
X  else retval = sinh(x) / cosh(x) ;
X  return(retval) ;
X}
X#endif   /*NEED_TANH*/
X
X
X/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved.
X *
X * Redistribution and use in source and binary forms are permitted
X * provided that the above copyright notice and this paragraph are
X * duplicated in all such forms and that any documentation,
X * advertising materials, and other materials related to such
X * distribution and use acknowledge that the software was developed
X * by the University of California, Berkeley.  The name of the
X * University may not be used to endorse or promote products derived
X * from this software without specific prior written permission.
X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
X * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
X *
X * All recipients should regard themselves as participants in an ongoing
X * research project and hence should feel obligated to report their
X * experiences (good or bad) with these elementary function codes, using
X * the sendbug(8) program, to the authors.
X */
X
X#ifdef   NEED_POW
X/* POW(X,Y)
X * RETURN X**Y
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 1/8/85;
X * REVISED BY K.C. NG on 7/10/85.
X *
X * Required system supported functions:
X *      scalb(x,n)
X *      logb(x)
X *      copysign(x,y)
X *      finite(x)
X *      drem(x,y)
X *
X * Required kernel functions:
X *      exp__E(a, c)    ...return  exp(a+c) - 1 - a*a/2
X *      log__L(x)       ...return  (log(1+x) - 2s)/s, s=x/(2+x)
X *      pow_p(x, y)     ...return  +(anything)**(finite non zero)
X *
X * Method
X *      1. Compute and return log(x) in three pieces:
X *              log(x) = n*ln2 + hi + lo,
X *         where n is an integer.
X *      2. Perform y * log(x) by simulating muti-precision arithmetic and
X *         return the answer in three pieces:
X *              y * log(x) = m * ln2 + hi + lo,
X *         where m is an integer.
X *      3. Return x ** y = exp(y * log(x))
X *              = 2^m * (exp(hi + lo)).
X *
X * Special cases:
X *      (anything) ** 0  is 1 ;
X *      (anything) ** 1  is itself;
X *      (anything) ** NaN is NaN;
X *      NaN ** (anything except 0) is NaN;
X *      +-(anything > 1) ** +INF is +INF;
X *      +-(anything > 1) ** -INF is +0;
X *      +-(anything < 1) ** +INF is +0;
X *      +-(anything < 1) ** -INF is +INF;
X *      +-1 ** +-INF is NaN and signal INVALID;
X *      +0 ** +(anything except 0, NaN)  is +0;
X *      -0 ** +(anything except 0, NaN, odd integer)  is +0;
X *      +0 ** -(anything except 0, NaN)  is +INF and signal DIV-BY-ZERO;
X *      -0 ** -(anything except 0, NaN, odd integer)  is +INF with signal;
X *      -0 ** (odd integer) = -( +0 ** (odd integer) );
X *      +INF ** +(anything except 0,NaN) is +INF;
X *      +INF ** -(anything except 0,NaN) is +0;
X *      -INF ** (odd integer) = -( +INF ** (odd integer) );
X *      -INF ** (even integer) = ( +INF ** (even integer) );
X *      -INF ** -(anything except integer,NaN) is NaN with signal;
X *      -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
X *      -(anything except 0) ** (non-integer) is NaN with signal;
X *
X * Accuracy:
X *      pow(x, y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
X *      and a Zilog Z8000,
X *                      pow(integer, integer)
X *      always returns the correct integer provided it is representable.
X *      In a test run with 100,000 random arguments with 0 < x, y < 20.0
X *      on a VAX, the maximum observed error was 1.79 ulps (units in the
X *      last place).
X *
X * Constants :
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe)        /* VAX D format */
X#include <errno.h>
Xextern double infnan() ;
X#ifdef vax
X#define _0x(A,B)        0x/**/A/**/B
X#else   /* vax */
X#define _0x(A,B)        0x/**/B/**/A
X#endif  /* vax */
X/* static double */
X/* ln2hi  =  6.9314718055829871446E-1    , Hex  2^  0   *  .B17217F7D00000 */
X/* ln2lo  =  1.6465949582897081279E-12   , Hex  2^-39   *  .E7BCD5E4F1D9CC */
X/* invln2 =  1.4426950408889634148E0     , Hex  2^  1   *  .B8AA3B295C17F1 */
X/* sqrt2  =  1.4142135623730950622E0     ; Hex  2^  1   *  .B504F333F9DE65 */
Xstatic long     ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)} ;
Xstatic long     ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)} ;
Xstatic long    invln2x[] = { _0x(aa3b,40b8), _0x(17f1,295c)} ;
Xstatic long     sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)} ;
X#define  ln2hi    (*(double*) ln2hix)
X#define  ln2lo    (*(double*) ln2lox)
X#define  invln2   (*(double*) invln2x)
X#define  sqrt2    (*(double*) sqrt2x)
X#else   /* defined(vax)||defined(tahoe) */
Xstatic double
Xln2hi  =  6.9314718036912381649E-1    , /* Hex  2^ -1   *  1.62E42FEE00000 */
Xln2lo  =  1.9082149292705877000E-10   , /* Hex  2^-33   *  1.A39EF35793C76 */
Xinvln2 =  1.4426950408889633870E0     , /* Hex  2^  0   *  1.71547652B82FE */
Xsqrt2  =  1.4142135623730951455E0     ; /* Hex  2^  0   *  1.6A09E667F3BCD */
X#endif  /* defined(vax) || defined(tahoe) */
X
Xstatic double zero = 0.0 ;
X
Xdouble
Xpow(x, y)
Xdouble x, y ;
X{
X  double drem(), pow_p(), copysign(), t ;
X  int finite() ;
X
X       if (y == 0.0) return(1.0) ;
X#if !defined(vax) && !defined(tahoe)
X  else if (y == 1.0 || x != x) return(x) ;   /* if x is NaN or y = 1 */
X#else
X  else if (y == 1.0) return(x) ;             /* if y = 1 */
X#endif  /* !defined(vax) && !defined(tahoe) */
X
X#if !defined(vax) && !defined(tahoe)
X  else if (y != y) return(y) ;               /* if y is NaN */
X#endif  /* !defined(vax) && !defined(tahoe) */
X  else if (!finite(y))                       /* if y is INF */
X         if ((t = copysign(x, 1.0)) == 1.0)
X           return(0.0 / zero) ;
X    else if (t > 1.0) return((y > 0.0) ? y : 0.0) ;
X    else return((y < 0.0) ? -y : 0.0) ;
X  else if (y == 2.0) return(x * x) ;
X  else if (y == -1.0) return(1.0 / x) ;
X
X/* sign(x) = 1 */
X
X  else if (copysign(1.0, x) == 1.0) return(pow_p(x, y)) ;
X
X/* sign(x)= -1 */
X/* if y is an even integer */
X
X  else if ((t = drem(y, 2.0)) == 0.0) return(pow_p(-x, y)) ;
X
X/* if y is an odd integer */
X
X  else if (copysign(t, 1.0) == 1.0) return(-pow_p(-x, y)) ;
X
X/* Henceforth y is not an integer */
X
X  else if (x == 0.0) return((y > 0.0) ? -x : 1.0 / (-x)) ;   /* x is -0 */
X  else                                  /* return NaN */
X    {
X#if defined(vax) || defined(tahoe)
X      return(infnan(EDOM)) ;            /* NaN */
X#else   /* defined(vax) || defined(tahoe) */
X      return(0.0 / zero) ;
X#endif  /* defined(vax) || defined(tahoe) */
X     }
X}
X
X
X/* pow_p(x,y) return x**y for x with sign = 1 and finite y */
X
Xstatic double
Xpow_p(x, y)
Xdouble x, y ;
X{
X  double logb(), scalb(), copysign(), log__L(), exp__E() ;
X  double c, s, t, z, tx, ty ;
X
X#ifdef tahoe
X  double tahoe_tmp ;
X#endif  /* tahoe */
X  float sx, sy ;
X  long k = 0 ;
X  int n, m ;
X
X  if (x == 0.0 || !finite(x))             /* if x is +INF or +0 */
X    {
X
X#if defined(vax) || defined(tahoe)
X      return((y > 0.0) ? x : infnan(ERANGE)) ;  /* if y < 0.0, return +INF */
X#else   /* defined(vax) || defined(tahoe) */
X      return((y > 0.0) ? x : 1.0 / x) ;
X#endif  /* defined(vax) || defined(tahoe) */
X     }
X  if (x == 1.0) return(x) ;     /* if x = 1.0, return 1 since y is finite */
X
X/* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
X
X    z = scalb(x, -(n = logb(x))) ;
X
X#if !defined(vax) && !defined(tahoe)      /* IEEE double; subnormal number */
X  if (n <= -1022)
X    {
X      n += (m = logb(z)) ;
X      z = scalb(z, -m) ;
X    }
X#endif  /* !defined(vax) && !defined(tahoe) */
X
X  if (z >= sqrt2)
X    {
X      n += 1 ;
X      z *= 0.5 ;
X    }
X  z -= 1.0 ;
X
X/* log(x) = nlog2 + log(1 + z) ~ nlog2 + t + tx */
X
X  s = z / (2.0 + z) ;
X  c = z * z * 0.5 ;
X  tx = s * (c + log__L(s * s)) ;
X  t = z - (c - tx) ;
X  tx += (z - t) - c ;
X
X/* if y * log(x) is neither too big nor too small */
X
X  if ((s = logb(y) + logb(n + t)) < 12.0)
X    if (s > -60.0)
X      {
X
X/* compute y * log(x) ~ mlog2 + t + c */
X
X        s = y * (n + invln2 * t) ;
X        m = s + copysign(0.5, s) ;      /* m := nint(y * log(x)) */
X        k = y ;
X        if ((double) k == y)            /* if y is an integer */
X          {
X            k = m - k * n ;
X            sx = t ;
X            tx += (t - sx) ;
X          }
X        else                            /* if y is not an integer */
X          {
X            k = m ;
X            tx += n * ln2lo ;
X            sx = (c = n * ln2hi) + t ;
X            tx += (c - sx) + t ;
X          }
X
X/* end of checking whether k == y */
X
X        sy = y ;
X        ty = y - sy ;                   /* y ~ sy + ty */
X
X#ifdef tahoe
X        s = (tahoe_tmp = sx) * sy - k * ln2hi ;
X#else   /* tahoe */
X        s = (double) sx * sy - k * ln2hi ;  /* (sy + ty) * (sx + tx) - kln2 */
X#endif  /* tahoe */
X
X        z = (tx * ty - k * ln2lo) ;
X        tx = tx * sy ; ty = sx * ty ;
X        t = ty + z ; t += tx ; t += s ;
X        c = -((((t-s)-tx)-ty)-z) ;
X
X/* return exp(y * log(x)) */
X
X        t += exp__E(t,c) ;
X        return(scalb(1.0 + t, m)) ;
X      }
X
X/* end of if log(y * log(x)) > -60.0 */
X
X  else                  /* exp(+- tiny) = 1 with inexact flag */
X    {
X      ln2hi + ln2lo ;
X      return(1.0) ;
X    }
X  else if (copysign(1.0, y) * (n + invln2 * t) < 0.0)
X    return(scalb(1.0, -5000)) ;     /* exp(-(big#)) underflows to zero */
X  else
X    return(scalb(1.0, 5000)) ;      /* exp(+(big#)) overflows to INF */
X}
X
X
X/* Some IEEE standard 754 recommended functions and remainder and sqrt for
X * supporting the C elementary functions.
X ******************************************************************************
X * WARNING:
X *      These codes are developed (in double) to support the C elementary
X * functions temporarily. They are not universal, and some of them are very
X * slow (in particular, drem and sqrt is extremely inefficient). Each
X * computer system should have its implementation of these functions using
X * its own assembler.
X ******************************************************************************
X *
X * IEEE 754 required operations:
X *     drem(x, p)
X *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
X *              nearest x/y; in half way case, choose the even one.
X *     sqrt(x)
X *              returns the square root of x correctly rounded according to
X *              the rounding mod.
X *
X * IEEE 754 recommended functions:
X * (a) copysign(x, y)
X *              returns x with the sign of y.
X * (b) scalb(x, N)
X *              returns  x * (2 ** N), for integer values N.
X * (c) logb(x)
X *              returns the unbiased exponent of x, a signed integer in
X *              double precision, except that logb(0) is -INF, logb(INF)
X *              is +INF, and logb(NAN) is that NAN.
X * (d) finite(x)
X *              returns the value TRUE if -INF < x < +INF and returns
X *              FALSE otherwise.
X *
X * CODED IN C BY K.C. NG, 11/25/84;
X * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
X */
X
X#if      defined(vax) || defined(tahoe)      /* VAX D format */
X  static unsigned short msign = 0x7fff, mexp = 0x7f80 ;
X  static short  prep1 = 57, gap = 7, bias = 129 ;
X  static double novf = 1.7E38, nunf = 3.0E-39 ;
X#else /* defined(vax) || defined(tahoe) */
X  static unsigned short msign = 0x7fff, mexp = 0x7ff0 ;
X  static short prep1 = 54, gap = 4, bias = 1023 ;
X  static double novf = 1.7E308, nunf = 3.0E-308 ;
X#endif  /* defined(vax) || defined(tahoe) */
X
Xdouble
Xscalb(x, N)
Xdouble x ;
Xint N ;
X{
X  int k ;
X  double scalb() ;
X
X#ifdef national
X  unsigned short *px = (unsigned short *) &x + 3 ;
X#else   /* national */
X  unsigned short *px = (unsigned short *) &x ;
X#endif  /* national */
X
X  if (x == 0.0) return(x) ;
X
X#if defined(vax) || defined(tahoe)
X
X  if ((k = *px & mexp) != ~msign)
X    {
X      if (N < -260) return(nunf * nunf) ;
X      else if (N > 260)
X        {
X          extern double infnan(), copysign() ;
X          return(copysign(infnan(ERANGE), x)) ;
X        }
X#else  /* defined(vax) || defined(tahoe) */
X
X  if ((k = *px & mexp) != mexp)
X    {
X           if (N < -2100) return(nunf * nunf) ;
X      else if (N > 2100) return(novf + novf) ;
X      if (k == 0)
X        {
X          x *= scalb(1.0, (int) prep1) ;
X          N -= prep1 ;
X          return(scalb(x, N)) ;
X        }
X#endif  /* defined(vax) || defined(tahoe) */
X
X      if ((k = (k >> gap) + N) > 0)
X        if (k < (mexp >> gap)) *px = (*px & ~mexp) | (k << gap) ;
X        else x = novf + novf ;               /* overflow */
X      else
X        if (k > -prep1)
X          {                                  /* gradual underflow */
X            *px = (*px & ~mexp) | (short) (1 << gap) ;
X            x *= scalb(1.0, k-1) ;
X          }
X        else return(nunf * nunf) ;
X    }
X  return(x) ;
X}
X
X
Xdouble
Xcopysign(x, y)
Xdouble x, y ;
X{
X#ifdef   national
X  unsigned short *px = (unsigned short *) &x + 3,
X                 *py = (unsigned short *) &y + 3 ;
X#else /* national */
X  unsigned short *px = (unsigned short *) &x,
X                 *py = (unsigned short *) &y ;
X#endif  /* national */
X
X#if defined(vax) || defined(tahoe)
X  if ((*px & mexp) == 0) return(x) ;
X#endif  /* defined(vax) || defined(tahoe) */
X
X  *px = (*px & msign) | (*py & ~msign) ;
X  return(x) ;
X}
X
X
Xdouble
Xlogb(x)
Xdouble x ;
X{
X#ifdef national
X  short *px = (short *) &x + 3, k ;
X#else   /* national */
X  short *px = (short *) &x, k ;
X#endif  /* national */
X
X#if defined(vax) || defined(tahoe)
X  return (int) (((*px & mexp) >> gap) - bias) ;
X#else   /* defined(vax) || defined(tahoe) */
X
X  if ((k = *px & mexp) != mexp)
X         if (k != 0)   return((k >> gap) - bias) ;
X    else if (x != 0.0) return(-1022.0) ;
X    else               return(-(1.0 / zero)) ;
X  else if (x != x) return(x) ;
X  else
X    {
X      *px &= msign ;
X      return(x) ;
X    }
X#endif  /* defined(vax) || defined(tahoe) */
X}
X
X
Xfinite(x)
Xdouble x ;
X{
X#if defined(vax) || defined(tahoe)
X  return(1) ;
X#else   /* defined(vax) || defined(tahoe) */
X#ifdef national
X  return((*((short *) &x + 3) & mexp) != mexp) ;
X#else   /* national */
X  return((*((short *) &x) & mexp) != mexp) ;
X#endif  /* national */
X#endif  /* defined(vax) || defined(tahoe) */
X}
X
X
Xdouble
Xdrem(x, p)
Xdouble x, p ;
X{
X  short sign ;
X  double hp, dp, tmp, drem(), scalb() ;
X  unsigned short k ;
X
X#ifdef national
X  unsigned short *px = (unsigned short *) &x   + 3,
X                 *pp = (unsigned short *) &p   + 3,
X                 *pd = (unsigned short *) &dp  + 3,
X                 *pt = (unsigned short *) &tmp + 3 ;
X#else   /* national */
X  unsigned short *px = (unsigned short *) &x,
X                 *pp = (unsigned short *) &p  ,
X                 *pd = (unsigned short *) &dp ,
X                 *pt = (unsigned short *) &tmp ;
X#endif  /* national */
X
X  *pp &= msign ;
X
X#if defined(vax) || defined(tahoe)
X  if ((*px & mexp) == ~msign)  /* is x a reserved operand? */
X#else   /* defined(vax) || defined(tahoe) */
X  if ((*px & mexp) == mexp)
X#endif  /* defined(vax) || defined(tahoe) */
X    return (x-p) - (x-p) ;     /* create nan if x is inf */
X
X  if (p == 0.0)
X    {
X      doerr("drem", "SINGULARITY", EDOM) ;
X#if defined(vax) || defined(tahoe)
X      extern double infnan() ;
X      return(infnan(EDOM)) ;
X#else   /* defined(vax) || defined(tahoe) */
X      return(0.0 / zero) ;
X#endif  /* defined(vax) || defined(tahoe) */
X    }
X
X#if defined(vax) || defined(tahoe)
X  if ((*pp & mexp) == ~msign)  /* is p a reserved operand? */
X#else   /* defined(vax) || defined(tahoe) */
X  if ((*pp & mexp) == mexp)
X#endif  /* defined(vax)||defined(tahoe) */
X    {
X      if (p != p) return p ;
X      else return x ;
X    }
X  else if (((*pp & mexp) >> gap) <= 1)
X
X/* subnormal p, or almost subnormal p */
X
X    {
X      double b ;
X      b = scalb(1.0, (int) prep1) ;
X      p *= b ;
X      x = drem(x, p) ;
X      x *= b ;
X      return(drem(x, p) / b) ;
X    }
X  else if (p >= novf / 2)
X    {
X      p /= 2 ;
X      x /= 2 ;
X      return(drem(x, p) * 2) ;
X    }
X  else
X    {
X      dp = p + p ;
X      hp = p / 2 ;
X      sign = *px & ~msign ;
X      *px &= msign ;
X      while (x > dp)
X        {
X          k = (*px & mexp) - (*pd & mexp) ;
X          tmp = dp ;
X          *pt += k ;
X
X#if defined(vax) || defined(tahoe)
X          if (x < tmp) *pt -= 128 ;
X#else /* defined(vax) || defined(tahoe) */
X          if (x < tmp) *pt -= 16 ;
X#endif  /* defined(vax) || defined(tahoe) */
X
X            x -= tmp ;
X        }
X      if (x > hp)
X        {
X          x -= p ;
X          if (x >= hp) x -= p ;
X        }
X
X#if defined(vax) || defined(tahoe)
X      if (x)
X#endif  /* defined(vax) || defined(tahoe) */
X        *px ^= sign ;
X      return(x) ;
X    }
X}
X
X
X/* exp__E(x,c)
X * ASSUMPTION: c << x  SO THAT  fl(x+c)=x.
X * (c is the correction term for x)
X * exp__E RETURNS
X *
X *                       /  exp(x+c) - 1 - x ,  1E-19 < |x| < .3465736
X *       exp__E(x,c) =  |
X *                       \  0 ,  |x| < 1E-19.
X *
X * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
X * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/31/85;
X * REVISED BY K.C. NG on 3/16/85, 4/16/85.
X *
X * Required system supported function:
X *      copysign(x,y)
X *
X * Method:
X *      1. Rational approximation. Let r=x+c.
X *         Based on
X *                                   2 * sinh(r/2)
X *                exp(r) - 1 =   ----------------------   ,
X *                               cosh(r/2) - sinh(r/2)
X *         exp__E(r) is computed using
X *                   x*x            (x/2)*W - ( Q - ( 2*P  + x*P ) )
X *                   --- + (c + x*[---------------------------------- + c ])
X *                    2                          1 - W
X *         where  P := p1*x^2 + p2*x^4,
X *                Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6)
X *                W := x/2-(Q-x*P),
X *
X *         (See the listing below for the values of p1,p2,q1,q2,q3. The poly-
X *          nomials P and Q may be regarded as the approximations to sinh
X *          and cosh :
X *              sinh(r/2) =  r/2 + r * P  ,  cosh(r/2) =  1 + Q . )
X *
X *         The coefficients were obtained by a special Remez algorithm.
X *
X * Approximation error:
X *
X *   |  exp(x) - 1                         |        2**(-57),  (IEEE double)
X *   | ------------  -  (exp__E(x,0)+x)/x  |  <=
X *   |       x                             |        2**(-69).  (VAX D)
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe)        /* VAX D format */
X#ifdef vax
X#define _0x(A,B)        0x/**/A/**/B
X#else   /* vax */
X#define _0x(A,B)        0x/**/B/**/A
X#endif  /* vax */
X
X/* static double */
X/* p1     =  1.5150724356786683059E-2    , Hex  2^ -6   *  .F83ABE67E1066A */
X/* p2     =  6.3112487873718332688E-5    , Hex  2^-13   *  .845B4248CD0173 */
X/* q1     =  1.1363478204690669916E-1    , Hex  2^ -3   *  .E8B95A44A2EC45 */
X/* q2     =  1.2624568129896839182E-3    , Hex  2^ -9   *  .A5790572E4F5E7 */
X/* q3     =  1.5021856115869022674E-6    ; Hex  2^-19   *  .C99EB4604AC395 */
X
Xstatic long        p1x[] = { _0x(3abe,3d78), _0x(066a,67e1) } ;
Xstatic long        p2x[] = { _0x(5b42,3984), _0x(0173,48cd) } ;
Xstatic long        q1x[] = { _0x(b95a,3ee8), _0x(ec45,44a2) } ;
Xstatic long        q2x[] = { _0x(7905,3ba5), _0x(f5e7,72e4) } ;
Xstatic long        q3x[] = { _0x(9eb4,36c9), _0x(c395,604a) } ;
X#define  p1  (*(double*) p1x)
X#define  p2  (*(double*) p2x)
X#define  q1  (*(double*) q1x)
X#define  q2  (*(double*) q2x)
X#define  q3  (*(double*) q3x)
X#else   /* defined(vax) || defined(tahoe) */
X
Xstatic double
Xp1 = 1.3887401997267371720E-2,     /* Hex  2^ -7   *  1.C70FF8B3CC2CF */
Xp2 = 3.3044019718331897649E-5,     /* Hex  2^-15   *  1.15317DF4526C4 */
Xq1 = 1.1110813732786649355E-1,     /* Hex  2^ -4   *  1.C719538248597 */
Xq2 = 9.9176615021572857300E-4 ;    /* Hex  2^-10   *  1.03FC4CB8C98E8 */
X#endif  /* defined(vax) || defined(tahoe) */
X
Xdouble
Xexp__E(x, c)
Xdouble x, c ;
X{
X  static double small = 1.0E-19 ;
X  double copysign(), z, p, q, xp, xh, w ;
X
X  if (copysign(x, 1.0) > small)
X     {
X       z = x * x ;
X       p = z * (p1 + z * p2) ;
X
X#if defined(vax) || defined(tahoe)
X       q = z * (q1 + z * (q2 + z * q3)) ;
X#else   /* defined(vax) || defined(tahoe) */
X       q = z * (q1 + z * q2) ;
X#endif  /* defined(vax) || defined(tahoe) */
X       xp = x * p ;
X       xh = x * 0.5 ;
X       w = xh - (q - xp) ;
X       p = p + p ;
X       c += x * ((xh * w - (q - (p + xp))) / (1.0 - w) + c) ;
X       return(z * 0.5 + c) ;
X     }
X
X/* end of |x| > small */
X
X   else
X     {
X       if (x != 0.0) 1.0 + small ;      /* raise the inexact flag */
X       return(copysign(0.0, x)) ;
X     }
X}
X
X
X/* log__L(Z)
X *              LOG(1+X) - 2S                          X
X * RETURN      ---------------  WHERE Z = S*S,  S = ------- , 0 <= Z <= .0294... *                    S                              2 + X
X *
X * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
X * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
X *
X * Method :
X *      1. Polynomial approximation: let s = x/(2+x).
X *         Based on log(1+x) = log(1+s) - log(1-s)
X *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
X *
X *         (log(1+x) - 2s)/s is computed by
X *
X *             z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
X *
X *         where z=s*s. (See the listing below for Lk's values.) The
X *         coefficients are obtained by a special Remez algorithm.
X *
X * Accuracy:
X *      Assuming no rounding error, the maximum magnitude of the approximation
X *      error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
X *      for VAX D format.
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe)        /* VAX D format (56 bits) */
X#ifdef vax
X#define _0x(A,B)        0x/**/A/**/B
X#else   /* vax */
X#define _0x(A,B)        0x/**/B/**/A
X#endif  /* vax */
X
X/* static double */
X/* L1     =  6.6666666666666703212E-1    , Hex  2^  0   *  .AAAAAAAAAAAAC5 */
X/* L2     =  3.9999999999970461961E-1    , Hex  2^ -1   *  .CCCCCCCCCC2684 */
X/* L3     =  2.8571428579395698188E-1    , Hex  2^ -1   *  .92492492F85782 */
X/* L4     =  2.2222221233634724402E-1    , Hex  2^ -2   *  .E38E3839B7AF2C */
X/* L5     =  1.8181879517064680057E-1    , Hex  2^ -2   *  .BA2EB4CC39655E */
X/* L6     =  1.5382888777946145467E-1    , Hex  2^ -2   *  .9D8551E8C5781D */
X/* L7     =  1.3338356561139403517E-1    , Hex  2^ -2   *  .8895B3907FCD92 */
X/* L8     =  1.2500000000000000000E-1    , Hex  2^ -2   *  .80000000000000 */
X
Xstatic long L1x[] = { _0x(aaaa,402a), _0x(aac5,aaaa) } ;
Xstatic long L2x[] = { _0x(cccc,3fcc), _0x(2684,cccc) } ;
Xstatic long L3x[] = { _0x(4924,3f92), _0x(5782,92f8) } ;
Xstatic long L4x[] = { _0x(8e38,3f63), _0x(af2c,39b7) } ;
Xstatic long L5x[] = { _0x(2eb4,3f3a), _0x(655e,cc39) } ;
Xstatic long L6x[] = { _0x(8551,3f1d), _0x(781d,e8c5) } ;
Xstatic long L7x[] = { _0x(95b3,3f08), _0x(cd92,907f) } ;
Xstatic long L8x[] = { _0x(0000,3f00), _0x(0000,0000) } ;
X
X#define  L1  (*(double*) L1x)
X#define  L2  (*(double*) L2x)
X#define  L3  (*(double*) L3x)
X#define  L4  (*(double*) L4x)
X#define  L5  (*(double*) L5x)
X#define  L6  (*(double*) L6x)
X#define  L7  (*(double*) L7x)
X#define  L8  (*(double*) L8x)
X#else   /* defined(vax) || defined(tahoe) */
X
Xstatic double
XL1 =  6.6666666666667340202E-1,       /* Hex  2^ -1   *  1.5555555555592 */
XL2 =  3.9999999999416702146E-1,       /* Hex  2^ -2   *  1.999999997FF24 */
XL3 =  2.8571428742008753154E-1,       /* Hex  2^ -2   *  1.24924941E07B4 */
XL4 =  2.2222198607186277597E-1,       /* Hex  2^ -3   *  1.C71C52150BEA6 */
XL5 =  1.8183562745289935658E-1,       /* Hex  2^ -3   *  1.74663CC94342F */
XL6 =  1.5314087275331442206E-1,       /* Hex  2^ -3   *  1.39A1EC014045B */
XL7 =  1.4795612545334174692E-1 ;      /* Hex  2^ -3   *  1.2F039F0085122 */
X#endif  /* defined(vax) || defined(tahoe) */
X
Xdouble
Xlog__L(z)
Xdouble z ;
X{
X#if defined(vax) || defined(tahoe)
X  return(z * (L1 + z * (L2 + z * (L3 + z * (L4 + z *
X             (L5 + z * (L6 + z * (L7 + z * L8)))))))) ;
X#else   /* defined(vax) || defined(tahoe) */
X  return(z * (L1 + z * (L2 + z * (L3 + z * (L4 + z *
X             (L5 + z * (L6 + z * L7))))))) ;
X#endif  /* defined(vax) || defined(tahoe) */
X}
X#endif /*NEED_POW*/
SHAR_EOF
echo "File mathlib.c is complete"
chmod 0444 mathlib.c || echo "restore of mathlib.c fails"
set `wc -c mathlib.c`;Sum=$1
if test "$Sum" != "66590"
then echo original size 66590, current size $Sum;fi
rm -f s2_seq_.tmp
echo "You have unpacked the last part"
exit 0