diff options
author | marha <marha@users.sourceforge.net> | 2012-11-19 10:38:33 +0100 |
---|---|---|
committer | marha <marha@users.sourceforge.net> | 2012-11-19 10:38:33 +0100 |
commit | 24635abae6008bef13e30d798b3f33abab412770 (patch) | |
tree | e799fbde24e0fd935af76b0bc48d30ef69f75d54 /libX11/src/xcms/cmsTrig.c | |
parent | e0844ae8b5ef87049537a7e0ebff81acc2695256 (diff) | |
parent | 6ce1d8f0f8c23e186175a7c84c21d7bfbe168dc5 (diff) | |
download | vcxsrv-24635abae6008bef13e30d798b3f33abab412770.tar.gz vcxsrv-24635abae6008bef13e30d798b3f33abab412770.tar.bz2 vcxsrv-24635abae6008bef13e30d798b3f33abab412770.zip |
Merge remote-tracking branch 'origin/released'
* origin/released:
Changed file permissions
dos -> unix
Conflicts:
libX11/include/X11/Xregion.h
libX11/src/ConvSel.c
libX11/src/CrGlCur.c
libX11/src/CrWindow.c
libX11/src/GetDflt.c
libX11/src/StrKeysym.c
libX11/src/Window.c
libX11/src/xkb/XKBBind.c
libX11/src/xkb/XKBGetMap.c
libX11/src/xkb/XKBSetGeom.c
libX11/src/xkb/XKBUse.c
libX11/src/xlibi18n/XimProto.h
libX11/src/xlibi18n/lcDynamic.c
libXdmcp/Key.c
libXdmcp/Write.c
libxcb/src/xcb_windefs.h
xkbcomp/keycodes.c
xkbcomp/xkbpath.c
xorg-server/hw/xwin/glx/winpriv.h
xorg-server/xkeyboard-config/rules/bin/ml1_s.sh
xorg-server/xkeyboard-config/rules/bin/ml1v1_s.sh
xorg-server/xkeyboard-config/rules/bin/ml1v_s.sh
xorg-server/xkeyboard-config/rules/bin/ml_s.sh
xorg-server/xkeyboard-config/rules/bin/mln_s.sh
xorg-server/xkeyboard-config/rules/bin/mlnvn_s.sh
xorg-server/xkeyboard-config/rules/bin/mlv_s.sh
xorg-server/xkeyboard-config/rules/compat/.gitignore
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 26ae05f75..5a01a56c8 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)); +} |