Official patch #6 for calctool v2.4 (Part 3 of 3).
Rich Burridge
richb at sunaus.oz
Fri Feb 2 13:22:25 AEST 1990
---- 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
More information about the Comp.sources.bugs
mailing list