diff options
Diffstat (limited to 'libX11/src/xcms/cmsTrig.c')
-rw-r--r-- | libX11/src/xcms/cmsTrig.c | 1186 |
1 files changed, 593 insertions, 593 deletions
diff --git a/libX11/src/xcms/cmsTrig.c b/libX11/src/xcms/cmsTrig.c index 5a01a56c8..26ae05f75 100644 --- a/libX11/src/xcms/cmsTrig.c +++ b/libX11/src/xcms/cmsTrig.c @@ -1,593 +1,593 @@ - -/* - * 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. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif -#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 (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(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)); -} +
+/*
+ * 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.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#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 (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(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));
+}
|