4.4BSD/usr/src/contrib/X11R5-lib/lib/X/XcmsTrig.c
/* $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));
}