InfoMagic Source Code 1993 July/THE_SOURCE_CODE_CD_ROM.iso/X/mit/lib/X/XcmsTrig.c

C/C++ Source or Header | 1993-07-21 | 14.7 KB | 592 lines

/* $XConsortium: XcmsTrig.c,v 1.4 91/11/05 11:02:20 rws Exp $" */

/*
 * Code and supporting documentation (c) Copyright 1990 1991 Tektronix, Inc.
 *     All Rights Reserved
 *
 * This file is a component of an X Window System-specific implementation
 * of Xcms based on the TekColor Color Management System.  Permission is
 * hereby granted to use, copy, modify, sell, and otherwise distribute this
 * software and its documentation for any purpose and without fee, provided
 * that this copyright, permission, and disclaimer notice is reproduced in
 * all copies of this software and in supporting documentation.  TekColor
 * is a trademark of Tektronix, Inc.
 *
 * Tektronix makes no representation about the suitability of this software
 * for any purpose.  It is provided "as is" and with all faults.
 *
 * TEKTRONIX DISCLAIMS ALL WARRANTIES APPLICABLE TO THIS SOFTWARE,
 * INCLUDING THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 * PARTICULAR PURPOSE.  IN NO EVENT SHALL TEKTRONIX BE LIABLE FOR ANY
 * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
 * RESULTING FROM LOSS OF USE, DATA, OR PROFITS, WHETHER IN AN ACTION OF
 * CONTRACT, NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR THE PERFORMANCE OF THIS SOFTWARE.
 */

/*
 *    It should be pointed out that for simplicity's sake, the
 *    environment parameters are defined as floating point constants,
 *    rather than octal or hexadecimal initializations of allocated
 *    storage areas.  This means that the range of allowed numbers
 *    may not exactly match the hardware's capabilities.  For example,
 *    if the maximum positive double precision floating point number
 *    is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
 *    defined to be 1.11E100 then the numbers between 1.11E100 and
 *    1.11...E100 are considered to be undefined.  For most
 *    applications, this will cause no problems.
 *
 *    An alternate method is to allocate a global static "double" variable,
 *    say "maxdouble", and use a union declaration and initialization
 *    to initialize it with the proper bits for the EXACT maximum value.
 *    This was not done because the only compilers available to the
 *    author did not fully support union initialization features.
 *
 */

/*
 *    EXTERNS
 */
extern double _XcmsSquareRoot();

/*
 *    FORWARD DECLARATIONS
 */
double _XcmsCosine();
static double _XcmsModulo();
static double _XcmsModuloF();
static double _XcmsPolynomial();
double _XcmsSine();
double _XcmsArcTangent();

/*
 *    DEFINES
 */

#define XCMS_MAXERROR           0.000001
#define XCMS_MAXITER           10000
#define XCMS_PI               3.14159265358979323846264338327950
#define XCMS_TWOPI         6.28318530717958620
#define XCMS_HALFPI        1.57079632679489660
#define XCMS_FOURTHPI        0.785398163397448280
#define XCMS_SIXTHPI        0.523598775598298820
#define XCMS_RADIANS(d)        ((d) * XCMS_PI / 180.0)
#define XCMS_DEGREES(r)        ((r) * 180.0 / XCMS_PI)
#define XCMS_X6_UNDERFLOWS    (4.209340e-52)    /* X**6 almost underflows */
#define XCMS_X16_UNDERFLOWS    (5.421010e-20)    /* X**16 almost underflows*/
#define XCMS_CHAR_BIT        8
#define XCMS_LONG_MAX        0x7FFFFFFF
#define XCMS_DEXPLEN        11
#define XCMS_NBITS(type)    (XCMS_CHAR_BIT * (int)sizeof(type))
#define XCMS_FABS(x)        ((x) < 0.0 ? -(x) : (x))

/* XCMS_DMAXPOWTWO - largest power of two exactly representable as a double */
#define XCMS_DMAXPOWTWO    ((double)(XCMS_LONG_MAX) * \
        (1L << (XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(long) + 1))

/*
 *    LOCAL VARIABLES
 */

static double cos_pcoeffs[] = {
    0.12905394659037374438e7,
   -0.37456703915723204710e6,
    0.13432300986539084285e5,
   -0.11231450823340933092e3
};

static double cos_qcoeffs[] = {
    0.12905394659037373590e7,
    0.23467773107245835052e5,
    0.20969518196726306286e3,
    1.0
};

static double sin_pcoeffs[] = {
    0.20664343336995858240e7,
   -0.18160398797407332550e6,
    0.35999306949636188317e4,
   -0.20107483294588615719e2
};

static double sin_qcoeffs[] = {
    0.26310659102647698963e7,
    0.39270242774649000308e5,
    0.27811919481083844087e3,
    1.0
};

/*
 *
 *  FUNCTION
 *
 *    _XcmsCosine   double precision cosine
 *
 *  KEY WORDS
 *
 *    cos
 *    machine independent routines
 *    trigonometric functions
 *    math libraries
 *
 *  DESCRIPTION
 *
 *    Returns double precision cosine of double precision
 *    floating point argument.
 *
 *  USAGE
 *
 *    double _XcmsCosine (x)
 *    double x;
 *
 *  REFERENCES
 *
 *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *    1968, pp. 112-120.
 *
 *  RESTRICTIONS
 *
 *    The sin and cos routines are interactive in the sense that
 *    in the process of reducing the argument to the range -PI/4
 *    to PI/4, each may call the other.  Ultimately one or the
 *    other uses a polynomial approximation on the reduced
 *    argument.  The sin approximation has a maximum relative error
 *    of 10**(-17.59) and the cos approximation has a maximum
 *    relative error of 10**(-16.18).
 *
 *    These error bounds assume exact arithmetic
 *    in the polynomial evaluation.  Additional rounding and
 *    truncation errors may occur as the argument is reduced
 *    to the range over which the polynomial approximation
 *    is valid, and as the polynomial is evaluated using
 *    finite-precision arithmetic.
 *
 *  PROGRAMMER
 *
 *    Fred Fish
 *
 *  INTERNALS
 *
 *    Computes cos(x) from:
 *
 *        (1)    Reduce argument x to range -PI to PI.
 *
 *        (2)    If x > PI/2 then call cos recursively
 *            using relation cos(x) = -cos(x - PI).
 *
 *        (3)    If x < -PI/2 then call cos recursively
 *            using relation cos(x) = -cos(x + PI).
 *
 *        (4)    If x > PI/4 then call sin using
 *            relation cos(x) = sin(PI/2 - x).
 *
 *        (5)    If x < -PI/4 then call cos using
 *            relation cos(x) = sin(PI/2 + x).
 *
 *        (6)    If x would cause underflow in approx
 *            evaluation arithmetic then return
 *            sqrt(1.0 - x**2).
 *
 *        (7)    By now x has been reduced to range
 *            -PI/4 to PI/4 and the approximation
 *            from HART pg. 119 can be used:
 *
 *            cos(x) = ( p(y) / q(y) )
 *            Where:
 *
 *            y    =    x * (4/PI)
 *
 *            p(y) =  SUM [ Pj * (y**(2*j)) ]
 *            over j = {0,1,2,3}
 *
 *            q(y) =  SUM [ Qj * (y**(2*j)) ]
 *            over j = {0,1,2,3}
 *
 *            P0   =    0.12905394659037374438571854e+7
 *            P1   =    -0.3745670391572320471032359e+6
 *            P2   =    0.134323009865390842853673e+5
 *            P3   =    -0.112314508233409330923e+3
 *            Q0   =    0.12905394659037373590295914e+7
 *            Q1   =    0.234677731072458350524124e+5
 *            Q2   =    0.2096951819672630628621e+3
 *            Q3   =    1.0000...
 *            (coefficients from HART table #3843 pg 244)
 *
 *
 *    **** NOTE ****    The range reduction relations used in
 *    this routine depend on the final approximation being valid
 *    over the negative argument range in addition to the positive
 *    argument range.  The particular approximation chosen from
 *    HART satisfies this requirement, although not explicitly
 *    stated in the text.  This may not be true of other
 *    approximations given in the reference.
 *
 */

double _XcmsCosine (x)
double x;
{
    auto double y;
    auto double yt2;
    double retval;

    if (x < -XCMS_PI || x > XCMS_PI) {
    x = _XcmsModulo (x, XCMS_TWOPI);
        if (x > XCMS_PI) {
        x = x - XCMS_TWOPI;
        } else if (x < -XCMS_PI) {
        x = x + XCMS_TWOPI;
        }
    }
    if (x > XCMS_HALFPI) {
    retval = -(_XcmsCosine (x - XCMS_PI));
    } else if (x < -XCMS_HALFPI) {
    retval = -(_XcmsCosine (x + XCMS_PI));
    } else if (x > XCMS_FOURTHPI) {
    retval = _XcmsSine (XCMS_HALFPI - x);
    } else if (x < -XCMS_FOURTHPI) {
    retval = _XcmsSine (XCMS_HALFPI + x);
    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
    retval = _XcmsSquareRoot (1.0 - (x * x));
    } else {
    y = x / XCMS_FOURTHPI;
    yt2 = y * y;
    retval = _XcmsPolynomial (3, cos_pcoeffs, yt2) / _XcmsPolynomial (3, cos_qcoeffs, yt2);
    }
    return (retval);
}


/*
 *  FUNCTION
 *
 *    _XcmsModulo   double precision modulo
 *
 *  KEY WORDS
 *
 *    _XcmsModulo
 *    machine independent routines
 *    math libraries
 *
 *  DESCRIPTION
 *
 *    Returns double precision modulo of two double
 *    precision arguments.
 *
 *  USAGE
 *
 *    double _XcmsModulo (value, base)
 *    double value;
 *    double base;
 *
 *  PROGRAMMER
 *
 *    Fred Fish
 *
 */
static double _XcmsModulo (value, base)
double value;
double base;
{
    auto double intpart;

    value /= base;
    value = _XcmsModuloF (value, &intpart);
    value *= base;
    return(value);
}


/*
 * frac = (double) _XcmsModuloF(double val, double *dp)
 *    return fractional part of 'val'
 *    set *dp to integer part of 'val'
 *
 * Note  -> only compiled for the CA or KA.  For the KB/MC,
 * "math.c" instantiates a copy of the inline function
 * defined in "math.h".
 */
static double
_XcmsModuloF(val, dp)
double val;
register double *dp;
{
    register double abs;
    register double ip;

    /* should check for illegal values here - nan, inf, etc */
    abs = XCMS_FABS(val);
    if (abs >= XCMS_DMAXPOWTWO) {
        ip = val;
    } else {
        ip = abs + XCMS_DMAXPOWTWO;    /* dump fraction */
        ip -= XCMS_DMAXPOWTWO;    /* restore w/o frac */
        if (ip > abs)        /* if it rounds up */
            ip -= 1.0;    /* fix it */
        ip = XCMS_FABS(ip);
    }
    *dp = ip;
    return (val - ip); /* signed fractional part */
}


/*
 *  FUNCTION
 *
 *    _XcmsPolynomial   double precision polynomial evaluation
 *
 *  KEY WORDS
 *
 *    poly
 *    machine independent routines
 *    math libraries
 *
 *  DESCRIPTION
 *
 *    Evaluates a polynomial and returns double precision
 *    result.  Is passed a the order of the polynomial,
 *    a pointer to an array of double precision polynomial
 *    coefficients (in ascending order), and the independent
 *    variable.
 *
 *  USAGE
 *
 *    double _XcmsPolynomial (order, coeffs, x)
 *    int order;
 *    double *coeffs;
 *    double x;
 *
 *  PROGRAMMER
 *
 *    Fred Fish
 *
 *  INTERNALS
 *
 *    Evalates the polynomial using recursion and the form:
 *
 *        P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
 *
 */

static double _XcmsPolynomial (order, coeffs, x)
register int order;
double *coeffs;
double x;
{
    auto double rtn_value;

#if 0
    auto double curr_coeff;
    if (order <= 0) {
    rtn_value = *coeffs;
    } else {
    curr_coeff = *coeffs;    /* Bug in Unisoft's compiler.  Does not */
    coeffs++;        /* generate good code for *coeffs++ */
    rtn_value = curr_coeff + x * _XcmsPolynomial (--order, coeffs, x);
    }
#else /* ++jrb -- removed tail recursion */
    coeffs += order;
    rtn_value = *coeffs--;
    while(order-- > 0)
    rtn_value = *coeffs-- + (x * rtn_value);
#endif

    return(rtn_value);
}


/*
 *  FUNCTION
 *
 *    _XcmsSine    double precision sine
 *
 *  KEY WORDS
 *
 *    sin
 *    machine independent routines
 *    trigonometric functions
 *    math libraries
 *
 *  DESCRIPTION
 *
 *    Returns double precision sine of double precision
 *    floating point argument.
 *
 *  USAGE
 *
 *    double _XcmsSine (x)
 *    double x;
 *
 *  REFERENCES
 *
 *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *    1968, pp. 112-120.
 *
 *  RESTRICTIONS
 *
 *    The sin and cos routines are interactive in the sense that
 *    in the process of reducing the argument to the range -PI/4
 *    to PI/4, each may call the other.  Ultimately one or the
 *    other uses a polynomial approximation on the reduced
 *    argument.  The sin approximation has a maximum relative error
 *    of 10**(-17.59) and the cos approximation has a maximum
 *    relative error of 10**(-16.18).
 *
 *    These error bounds assume exact arithmetic
 *    in the polynomial evaluation.  Additional rounding and
 *    truncation errors may occur as the argument is reduced
 *    to the range over which the polynomial approximation
 *    is valid, and as the polynomial is evaluated using
 *    finite-precision arithmetic.
 *
 *  PROGRAMMER
 *
 *    Fred Fish
 *
 *  INTERNALS
 *
 *    Computes sin(x) from:
 *
 *      (1)    Reduce argument x to range -PI to PI.
 *
 *      (2)    If x > PI/2 then call sin recursively
 *        using relation sin(x) = -sin(x - PI).
 *
 *      (3)    If x < -PI/2 then call sin recursively
 *        using relation sin(x) = -sin(x + PI).
 *
 *      (4)    If x > PI/4 then call cos using
 *        relation sin(x) = cos(PI/2 - x).
 *
 *      (5)    If x < -PI/4 then call cos using
 *        relation sin(x) = -cos(PI/2 + x).
 *
 *      (6)    If x is small enough that polynomial
 *        evaluation would cause underflow
 *        then return x, since sin(x)
 *        approaches x as x approaches zero.
 *
 *      (7)    By now x has been reduced to range
 *        -PI/4 to PI/4 and the approximation
 *        from HART pg. 118 can be used:
 *
 *        sin(x) = y * ( p(y) / q(y) )
 *        Where:
 *
 *        y    =  x * (4/PI)
 *
 *        p(y) =  SUM [ Pj * (y**(2*j)) ]
 *        over j = {0,1,2,3}
 *
 *        q(y) =  SUM [ Qj * (y**(2*j)) ]
 *        over j = {0,1,2,3}
 *
 *        P0   =  0.206643433369958582409167054e+7
 *        P1   =  -0.18160398797407332550219213e+6
 *        P2   =  0.359993069496361883172836e+4
 *        P3   =  -0.2010748329458861571949e+2
 *        Q0   =  0.263106591026476989637710307e+7
 *        Q1   =  0.3927024277464900030883986e+5
 *        Q2   =  0.27811919481083844087953e+3
 *        Q3   =  1.0000...
 *        (coefficients from HART table #3063 pg 234)
 *
 *
 *    **** NOTE ****      The range reduction relations used in
 *    this routine depend on the final approximation being valid
 *    over the negative argument range in addition to the positive
 *    argument range.  The particular approximation chosen from
 *    HART satisfies this requirement, although not explicitly
 *    stated in the text.  This may not be true of other
 *    approximations given in the reference.
 *
 */

double
_XcmsSine (x)
double x;
{
    double y;
    double yt2;
    double retval;

    if (x < -XCMS_PI || x > XCMS_PI) {
    x = _XcmsModulo (x, XCMS_TWOPI);
    if (x > XCMS_PI) {
        x = x - XCMS_TWOPI;
    } else if (x < -XCMS_PI) {
        x = x + XCMS_TWOPI;
    }
    }
    if (x > XCMS_HALFPI) {
    retval = -(_XcmsSine (x - XCMS_PI));
    } else if (x < -XCMS_HALFPI) {
    retval = -(_XcmsSine (x + XCMS_PI));
    } else if (x > XCMS_FOURTHPI) {
    retval = _XcmsCosine (XCMS_HALFPI - x);
    } else if (x < -XCMS_FOURTHPI) {
    retval = -(_XcmsCosine (XCMS_HALFPI + x));
    } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
    retval = x;
    } else {
    y = x / XCMS_FOURTHPI;
    yt2 = y * y;
    retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
    }
    return(retval);
}


/*
 *    NAME
 *        _XcmsArcTangent
 *
 *    SYNOPSIS
 */
double
_XcmsArcTangent(x)
    double x;
/*
 *    DESCRIPTION
 *        Computes the arctangent.
 *        This is an implementation of the Gauss algorithm as
 *        described in:
 *            Forman S. Acton, Numerical Methods That Work,
 *            New York, NY, Harper & Row, 1970.
 *
 *    RETURNS
 *        Returns the arctangent
 */
{
    double ai, a1, bi, b1, l, d;
    double maxerror;
    int i;

    if (x == 0.0)  {
    return (0.0);
    }
    if (x < 1.0) {
    maxerror = x * XCMS_MAXERROR;
    } else {
    maxerror = XCMS_MAXERROR;
    }
    ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
    bi = 1.0;
    for (i = 0; i < XCMS_MAXITER; i++) {
    a1 = (ai + bi) / 2.0;
    b1 = _XcmsSquareRoot((a1 * bi));
    if (a1 == b1)
        break;
    d = XCMS_FABS(a1 - b1);
    if (d < maxerror)
        break;
    ai = a1;
    bi = b1;
    }

    l = ((a1 > b1) ? b1 : a1);

    a1 = _XcmsSquareRoot(1 + (x * x));
    return (x / (a1 * l));
}