mirror of
https://github.com/Stichting-MINIX-Research-Foundation/xsrc.git
synced 2025-09-09 12:49:44 -04:00
595 lines
15 KiB
C
595 lines
15 KiB
C
/* $Xorg: cmsTrig.c,v 1.3 2000/08/17 19:45:10 cpqbld 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.
|
||
*/
|
||
|
||
/* $XFree86: xc/lib/X11/cmsTrig.c,v 3.9 2003/04/13 19:22:20 dawes Exp $ */
|
||
/*
|
||
* 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.
|
||
*
|
||
*/
|
||
|
||
#include "Xcmsint.h"
|
||
|
||
/* forward/static */
|
||
static double _XcmsModulo(double value, double base);
|
||
static double _XcmsPolynomial(
|
||
register int order,
|
||
double const *coeffs,
|
||
double x);
|
||
static double
|
||
_XcmsModuloF(
|
||
double val,
|
||
register double *dp);
|
||
|
||
/*
|
||
* 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 */
|
||
#ifdef _CRAY
|
||
#define XCMS_DMAXPOWTWO ((double)(1 < 47))
|
||
#else
|
||
#define XCMS_DMAXPOWTWO ((double)(XCMS_LONG_MAX) * \
|
||
(1L << ((XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(int) + 1)))
|
||
#endif
|
||
|
||
/*
|
||
* LOCAL VARIABLES
|
||
*/
|
||
|
||
static double const cos_pcoeffs[] = {
|
||
0.12905394659037374438e7,
|
||
-0.37456703915723204710e6,
|
||
0.13432300986539084285e5,
|
||
-0.11231450823340933092e3
|
||
};
|
||
|
||
static double const cos_qcoeffs[] = {
|
||
0.12905394659037373590e7,
|
||
0.23467773107245835052e5,
|
||
0.20969518196726306286e3,
|
||
1.0
|
||
};
|
||
|
||
static double const sin_pcoeffs[] = {
|
||
0.20664343336995858240e7,
|
||
-0.18160398797407332550e6,
|
||
0.35999306949636188317e4,
|
||
-0.20107483294588615719e2
|
||
};
|
||
|
||
static double const 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(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(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(
|
||
double val,
|
||
register double *dp)
|
||
{
|
||
register double abs;
|
||
/*
|
||
* Don't use a register for this. The extra precision this results
|
||
* in on some systems causes problems.
|
||
*/
|
||
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(
|
||
register int order,
|
||
double const *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 = 0.0, bi, b1 = 0.0, 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));
|
||
}
|