Used synchronise script to update files

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));

+}