home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / evbl0627.zip / everblue_20010627.zip / x11 / Xcms_Trig.c < prev    next >
C/C++ Source or Header  |  1999-11-02  |  16KB  |  611 lines
  1. /* $XConsortium: cmsTrig.c,v 1.7 95/06/08 23:20:39 gildea Exp $" */
  2. /* $XFree86: xc/lib/X11/cmsTrig.c,v 3.2 1996/05/06 05:54:11 dawes Exp $ */
  3.  
  4. /*
  5.  * Code and supporting documentation (c) Copyright 1990 1991 Tektronix, Inc.
  6.  *     All Rights Reserved
  7.  * 
  8.  * This file is a component of an X Window System-specific implementation
  9.  * of Xcms based on the TekColor Color Management System.  Permission is
  10.  * hereby granted to use, copy, modify, sell, and otherwise distribute this
  11.  * software and its documentation for any purpose and without fee, provided
  12.  * that this copyright, permission, and disclaimer notice is reproduced in
  13.  * all copies of this software and in supporting documentation.  TekColor
  14.  * is a trademark of Tektronix, Inc.
  15.  * 
  16.  * Tektronix makes no representation about the suitability of this software
  17.  * for any purpose.  It is provided "as is" and with all faults.
  18.  * 
  19.  * TEKTRONIX DISCLAIMS ALL WARRANTIES APPLICABLE TO THIS SOFTWARE,
  20.  * INCLUDING THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
  21.  * PARTICULAR PURPOSE.  IN NO EVENT SHALL TEKTRONIX BE LIABLE FOR ANY
  22.  * SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER
  23.  * RESULTING FROM LOSS OF USE, DATA, OR PROFITS, WHETHER IN AN ACTION OF
  24.  * CONTRACT, NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
  25.  * CONNECTION WITH THE USE OR THE PERFORMANCE OF THIS SOFTWARE.
  26.  */
  27.  
  28. /*
  29.  *    It should be pointed out that for simplicity's sake, the
  30.  *    environment parameters are defined as floating point constants,
  31.  *    rather than octal or hexadecimal initializations of allocated
  32.  *    storage areas.  This means that the range of allowed numbers
  33.  *    may not exactly match the hardware's capabilities.  For example,
  34.  *    if the maximum positive double precision floating point number
  35.  *    is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
  36.  *    defined to be 1.11E100 then the numbers between 1.11E100 and
  37.  *    1.11...E100 are considered to be undefined.  For most
  38.  *    applications, this will cause no problems.
  39.  *
  40.  *    An alternate method is to allocate a global static "double" variable,
  41.  *    say "maxdouble", and use a union declaration and initialization
  42.  *    to initialize it with the proper bits for the EXACT maximum value.
  43.  *    This was not done because the only compilers available to the
  44.  *    author did not fully support union initialization features.
  45.  *
  46.  */
  47.  
  48. /*
  49.  *    EXTERNS
  50.  */
  51. extern double _XcmsSquareRoot();
  52.  
  53. /*
  54.  *    FORWARD DECLARATIONS
  55.  */
  56. double _XcmsCosine();
  57. static double _XcmsModulo();
  58. static double _XcmsModuloF();
  59. static double _XcmsPolynomial();
  60. double _XcmsSine();
  61. double _XcmsArcTangent();
  62.  
  63. /*
  64.  *    DEFINES
  65.  */
  66.  
  67. #ifdef __STDC__
  68. #define Const const
  69. #else
  70. #define Const /**/
  71. #endif
  72. #define XCMS_MAXERROR           0.000001
  73. #define XCMS_MAXITER           10000
  74. #define XCMS_PI               3.14159265358979323846264338327950
  75. #define XCMS_TWOPI         6.28318530717958620
  76. #define XCMS_HALFPI        1.57079632679489660
  77. #define XCMS_FOURTHPI        0.785398163397448280
  78. #define XCMS_SIXTHPI        0.523598775598298820
  79. #define XCMS_RADIANS(d)        ((d) * XCMS_PI / 180.0)
  80. #define XCMS_DEGREES(r)        ((r) * 180.0 / XCMS_PI)
  81. #define XCMS_X6_UNDERFLOWS    (4.209340e-52)    /* X**6 almost underflows */
  82. #define XCMS_X16_UNDERFLOWS    (5.421010e-20)    /* X**16 almost underflows*/
  83. #define XCMS_CHAR_BIT        8
  84. #define XCMS_LONG_MAX        0x7FFFFFFF
  85. #define XCMS_DEXPLEN        11
  86. #define XCMS_NBITS(type)    (XCMS_CHAR_BIT * (int)sizeof(type))
  87. #define XCMS_FABS(x)        ((x) < 0.0 ? -(x) : (x))
  88.  
  89. /* XCMS_DMAXPOWTWO - largest power of two exactly representable as a double */
  90. #ifdef _CRAY
  91. #define XCMS_DMAXPOWTWO    ((double)(1 < 47))
  92. #else
  93. #ifdef __alpha__
  94. #define XCMS_DMAXPOWTWO    ((double)(XCMS_LONG_MAX) * \
  95.         (1L << (XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(int) + 1))
  96. #else
  97. #define XCMS_DMAXPOWTWO    ((double)(XCMS_LONG_MAX) * \
  98.         (1L << (XCMS_NBITS(double)-XCMS_DEXPLEN) - XCMS_NBITS(long) + 1))
  99. #endif
  100. #endif
  101.  
  102. /*
  103.  *    LOCAL VARIABLES
  104.  */
  105.  
  106. static double Const cos_pcoeffs[] = {
  107.     0.12905394659037374438e7,
  108.    -0.37456703915723204710e6,
  109.     0.13432300986539084285e5,
  110.    -0.11231450823340933092e3
  111. };
  112.  
  113. static double Const cos_qcoeffs[] = {
  114.     0.12905394659037373590e7,
  115.     0.23467773107245835052e5,
  116.     0.20969518196726306286e3,
  117.     1.0
  118. };
  119.  
  120. static double Const sin_pcoeffs[] = {
  121.     0.20664343336995858240e7,
  122.    -0.18160398797407332550e6,
  123.     0.35999306949636188317e4,
  124.    -0.20107483294588615719e2
  125. };
  126.  
  127. static double Const sin_qcoeffs[] = {
  128.     0.26310659102647698963e7,
  129.     0.39270242774649000308e5,
  130.     0.27811919481083844087e3,
  131.     1.0
  132. };
  133.  
  134. /*
  135.  *
  136.  *  FUNCTION
  137.  *
  138.  *    _XcmsCosine   double precision cosine
  139.  *
  140.  *  KEY WORDS
  141.  *
  142.  *    cos
  143.  *    machine independent routines
  144.  *    trigonometric functions
  145.  *    math libraries
  146.  *
  147.  *  DESCRIPTION
  148.  *
  149.  *    Returns double precision cosine of double precision
  150.  *    floating point argument.
  151.  *
  152.  *  USAGE
  153.  *
  154.  *    double _XcmsCosine (x)
  155.  *    double x;
  156.  *
  157.  *  REFERENCES
  158.  *
  159.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  160.  *    1968, pp. 112-120.
  161.  *
  162.  *  RESTRICTIONS
  163.  *
  164.  *    The sin and cos routines are interactive in the sense that
  165.  *    in the process of reducing the argument to the range -PI/4
  166.  *    to PI/4, each may call the other.  Ultimately one or the
  167.  *    other uses a polynomial approximation on the reduced
  168.  *    argument.  The sin approximation has a maximum relative error
  169.  *    of 10**(-17.59) and the cos approximation has a maximum
  170.  *    relative error of 10**(-16.18).
  171.  *
  172.  *    These error bounds assume exact arithmetic
  173.  *    in the polynomial evaluation.  Additional rounding and
  174.  *    truncation errors may occur as the argument is reduced
  175.  *    to the range over which the polynomial approximation
  176.  *    is valid, and as the polynomial is evaluated using
  177.  *    finite-precision arithmetic.
  178.  *    
  179.  *  PROGRAMMER
  180.  *
  181.  *    Fred Fish
  182.  *
  183.  *  INTERNALS
  184.  *
  185.  *    Computes cos(x) from:
  186.  *
  187.  *        (1)    Reduce argument x to range -PI to PI.
  188.  *
  189.  *        (2)    If x > PI/2 then call cos recursively
  190.  *            using relation cos(x) = -cos(x - PI).
  191.  *
  192.  *        (3)    If x < -PI/2 then call cos recursively
  193.  *            using relation cos(x) = -cos(x + PI).
  194.  *
  195.  *        (4)    If x > PI/4 then call sin using
  196.  *            relation cos(x) = sin(PI/2 - x).
  197.  *
  198.  *        (5)    If x < -PI/4 then call cos using
  199.  *            relation cos(x) = sin(PI/2 + x).
  200.  *
  201.  *        (6)    If x would cause underflow in approx
  202.  *            evaluation arithmetic then return
  203.  *            sqrt(1.0 - x**2).
  204.  *
  205.  *        (7)    By now x has been reduced to range
  206.  *            -PI/4 to PI/4 and the approximation
  207.  *            from HART pg. 119 can be used:
  208.  *
  209.  *            cos(x) = ( p(y) / q(y) )
  210.  *            Where:
  211.  *
  212.  *            y    =    x * (4/PI)
  213.  *
  214.  *            p(y) =  SUM [ Pj * (y**(2*j)) ]
  215.  *            over j = {0,1,2,3}
  216.  *
  217.  *            q(y) =  SUM [ Qj * (y**(2*j)) ]
  218.  *            over j = {0,1,2,3}
  219.  *
  220.  *            P0   =    0.12905394659037374438571854e+7
  221.  *            P1   =    -0.3745670391572320471032359e+6
  222.  *            P2   =    0.134323009865390842853673e+5
  223.  *            P3   =    -0.112314508233409330923e+3
  224.  *            Q0   =    0.12905394659037373590295914e+7
  225.  *            Q1   =    0.234677731072458350524124e+5
  226.  *            Q2   =    0.2096951819672630628621e+3
  227.  *            Q3   =    1.0000...
  228.  *            (coefficients from HART table #3843 pg 244)
  229.  *
  230.  *
  231.  *    **** NOTE ****    The range reduction relations used in
  232.  *    this routine depend on the final approximation being valid
  233.  *    over the negative argument range in addition to the positive
  234.  *    argument range.  The particular approximation chosen from
  235.  *    HART satisfies this requirement, although not explicitly
  236.  *    stated in the text.  This may not be true of other
  237.  *    approximations given in the reference.
  238.  *            
  239.  */
  240.  
  241. double _XcmsCosine (x)
  242. double x;
  243. {
  244.     auto double y;
  245.     auto double yt2;
  246.     double retval;
  247.  
  248.     if (x < -XCMS_PI || x > XCMS_PI) {
  249.     x = _XcmsModulo (x, XCMS_TWOPI);
  250.         if (x > XCMS_PI) {
  251.         x = x - XCMS_TWOPI;
  252.         } else if (x < -XCMS_PI) {
  253.         x = x + XCMS_TWOPI;
  254.         }
  255.     }
  256.     if (x > XCMS_HALFPI) {
  257.     retval = -(_XcmsCosine (x - XCMS_PI));
  258.     } else if (x < -XCMS_HALFPI) {
  259.     retval = -(_XcmsCosine (x + XCMS_PI));
  260.     } else if (x > XCMS_FOURTHPI) {
  261.     retval = _XcmsSine (XCMS_HALFPI - x);
  262.     } else if (x < -XCMS_FOURTHPI) {
  263.     retval = _XcmsSine (XCMS_HALFPI + x);
  264.     } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
  265.     retval = _XcmsSquareRoot (1.0 - (x * x));
  266.     } else {
  267.     y = x / XCMS_FOURTHPI;
  268.     yt2 = y * y;
  269.     retval = _XcmsPolynomial (3, cos_pcoeffs, yt2) / _XcmsPolynomial (3, cos_qcoeffs, yt2);
  270.     }
  271.     return (retval);
  272. }
  273.  
  274.  
  275. /*
  276.  *  FUNCTION
  277.  *
  278.  *    _XcmsModulo   double precision modulo
  279.  *
  280.  *  KEY WORDS
  281.  *
  282.  *    _XcmsModulo
  283.  *    machine independent routines
  284.  *    math libraries
  285.  *
  286.  *  DESCRIPTION
  287.  *
  288.  *    Returns double precision modulo of two double
  289.  *    precision arguments.
  290.  *
  291.  *  USAGE
  292.  *
  293.  *    double _XcmsModulo (value, base)
  294.  *    double value;
  295.  *    double base;
  296.  *
  297.  *  PROGRAMMER
  298.  *
  299.  *    Fred Fish
  300.  *
  301.  */
  302. static double _XcmsModulo (value, base)
  303. double value;
  304. double base;
  305. {
  306.     auto double intpart;
  307.  
  308.     value /= base;
  309.     value = _XcmsModuloF (value, &intpart);
  310.     value *= base;
  311.     return(value);
  312. }
  313.  
  314.  
  315. /*
  316.  * frac = (double) _XcmsModuloF(double val, double *dp)
  317.  *    return fractional part of 'val'
  318.  *    set *dp to integer part of 'val'
  319.  *
  320.  * Note  -> only compiled for the CA or KA.  For the KB/MC,
  321.  * "math.c" instantiates a copy of the inline function
  322.  * defined in "math.h".
  323.  */
  324. static double
  325. _XcmsModuloF(val, dp)
  326. double val;
  327. register double *dp;
  328. {
  329.     register double abs;
  330.     /*
  331.      * Don't use a register for this.  The extra precision this results
  332.      * in on some systems causes problems.
  333.      */
  334.     double ip;
  335.  
  336.     /* should check for illegal values here - nan, inf, etc */
  337.     abs = XCMS_FABS(val);
  338.     if (abs >= XCMS_DMAXPOWTWO) {
  339.         ip = val;
  340.     } else {
  341.         ip = abs + XCMS_DMAXPOWTWO;    /* dump fraction */
  342.         ip -= XCMS_DMAXPOWTWO;    /* restore w/o frac */
  343.         if (ip > abs)        /* if it rounds up */
  344.             ip -= 1.0;    /* fix it */
  345.         ip = XCMS_FABS(ip);
  346.     }
  347.     *dp = ip;
  348.     return (val - ip); /* signed fractional part */
  349. }
  350.  
  351.  
  352. /*
  353.  *  FUNCTION
  354.  *
  355.  *    _XcmsPolynomial   double precision polynomial evaluation
  356.  *
  357.  *  KEY WORDS
  358.  *
  359.  *    poly
  360.  *    machine independent routines
  361.  *    math libraries
  362.  *
  363.  *  DESCRIPTION
  364.  *
  365.  *    Evaluates a polynomial and returns double precision
  366.  *    result.  Is passed a the order of the polynomial,
  367.  *    a pointer to an array of double precision polynomial
  368.  *    coefficients (in ascending order), and the independent
  369.  *    variable.
  370.  *
  371.  *  USAGE
  372.  *
  373.  *    double _XcmsPolynomial (order, coeffs, x)
  374.  *    int order;
  375.  *    double *coeffs;
  376.  *    double x;
  377.  *
  378.  *  PROGRAMMER
  379.  *
  380.  *    Fred Fish
  381.  *
  382.  *  INTERNALS
  383.  *
  384.  *    Evalates the polynomial using recursion and the form:
  385.  *
  386.  *        P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
  387.  *
  388.  */
  389.  
  390. static double _XcmsPolynomial (order, coeffs, x)
  391. register int order;
  392. double Const *coeffs;
  393. double x;
  394. {
  395.     auto double rtn_value;
  396.  
  397. #if 0
  398.     auto double curr_coeff;
  399.     if (order <= 0) {
  400.     rtn_value = *coeffs;
  401.     } else {
  402.     curr_coeff = *coeffs;    /* Bug in Unisoft's compiler.  Does not */
  403.     coeffs++;        /* generate good code for *coeffs++ */
  404.     rtn_value = curr_coeff + x * _XcmsPolynomial (--order, coeffs, x);
  405.     }
  406. #else /* ++jrb -- removed tail recursion */
  407.     coeffs += order;
  408.     rtn_value = *coeffs--;
  409.     while(order-- > 0)
  410.     rtn_value = *coeffs-- + (x * rtn_value);
  411. #endif
  412.  
  413.     return(rtn_value);
  414. }
  415.  
  416.  
  417. /*
  418.  *  FUNCTION
  419.  *
  420.  *    _XcmsSine    double precision sine
  421.  *
  422.  *  KEY WORDS
  423.  *
  424.  *    sin
  425.  *    machine independent routines
  426.  *    trigonometric functions
  427.  *    math libraries
  428.  *
  429.  *  DESCRIPTION
  430.  *
  431.  *    Returns double precision sine of double precision
  432.  *    floating point argument.
  433.  *
  434.  *  USAGE
  435.  *
  436.  *    double _XcmsSine (x)
  437.  *    double x;
  438.  *
  439.  *  REFERENCES
  440.  *
  441.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  442.  *    1968, pp. 112-120.
  443.  *
  444.  *  RESTRICTIONS
  445.  *
  446.  *    The sin and cos routines are interactive in the sense that
  447.  *    in the process of reducing the argument to the range -PI/4
  448.  *    to PI/4, each may call the other.  Ultimately one or the
  449.  *    other uses a polynomial approximation on the reduced
  450.  *    argument.  The sin approximation has a maximum relative error
  451.  *    of 10**(-17.59) and the cos approximation has a maximum
  452.  *    relative error of 10**(-16.18).
  453.  *
  454.  *    These error bounds assume exact arithmetic
  455.  *    in the polynomial evaluation.  Additional rounding and
  456.  *    truncation errors may occur as the argument is reduced
  457.  *    to the range over which the polynomial approximation
  458.  *    is valid, and as the polynomial is evaluated using
  459.  *    finite-precision arithmetic.
  460.  *    
  461.  *  PROGRAMMER
  462.  *
  463.  *    Fred Fish
  464.  *
  465.  *  INTERNALS
  466.  *
  467.  *    Computes sin(x) from:
  468.  *
  469.  *      (1)    Reduce argument x to range -PI to PI.
  470.  *
  471.  *      (2)    If x > PI/2 then call sin recursively
  472.  *        using relation sin(x) = -sin(x - PI).
  473.  *
  474.  *      (3)    If x < -PI/2 then call sin recursively
  475.  *        using relation sin(x) = -sin(x + PI).
  476.  *
  477.  *      (4)    If x > PI/4 then call cos using
  478.  *        relation sin(x) = cos(PI/2 - x).
  479.  *
  480.  *      (5)    If x < -PI/4 then call cos using
  481.  *        relation sin(x) = -cos(PI/2 + x).
  482.  *
  483.  *      (6)    If x is small enough that polynomial
  484.  *        evaluation would cause underflow
  485.  *        then return x, since sin(x)
  486.  *        approaches x as x approaches zero.
  487.  *
  488.  *      (7)    By now x has been reduced to range
  489.  *        -PI/4 to PI/4 and the approximation
  490.  *        from HART pg. 118 can be used:
  491.  *
  492.  *        sin(x) = y * ( p(y) / q(y) )
  493.  *        Where:
  494.  *
  495.  *        y    =  x * (4/PI)
  496.  *
  497.  *        p(y) =  SUM [ Pj * (y**(2*j)) ]
  498.  *        over j = {0,1,2,3}
  499.  *
  500.  *        q(y) =  SUM [ Qj * (y**(2*j)) ]
  501.  *        over j = {0,1,2,3}
  502.  *
  503.  *        P0   =  0.206643433369958582409167054e+7
  504.  *        P1   =  -0.18160398797407332550219213e+6
  505.  *        P2   =  0.359993069496361883172836e+4
  506.  *        P3   =  -0.2010748329458861571949e+2
  507.  *        Q0   =  0.263106591026476989637710307e+7
  508.  *        Q1   =  0.3927024277464900030883986e+5
  509.  *        Q2   =  0.27811919481083844087953e+3
  510.  *        Q3   =  1.0000...
  511.  *        (coefficients from HART table #3063 pg 234)
  512.  *
  513.  *
  514.  *    **** NOTE ****      The range reduction relations used in
  515.  *    this routine depend on the final approximation being valid
  516.  *    over the negative argument range in addition to the positive
  517.  *    argument range.  The particular approximation chosen from
  518.  *    HART satisfies this requirement, although not explicitly
  519.  *    stated in the text.  This may not be true of other
  520.  *    approximations given in the reference.
  521.  *            
  522.  */
  523.  
  524. double
  525. _XcmsSine (x)
  526. double x;
  527. {
  528.     double y;
  529.     double yt2;
  530.     double retval;
  531.  
  532.     if (x < -XCMS_PI || x > XCMS_PI) {
  533.     x = _XcmsModulo (x, XCMS_TWOPI);
  534.     if (x > XCMS_PI) {
  535.         x = x - XCMS_TWOPI;
  536.     } else if (x < -XCMS_PI) {
  537.         x = x + XCMS_TWOPI;
  538.     }
  539.     }
  540.     if (x > XCMS_HALFPI) {
  541.     retval = -(_XcmsSine (x - XCMS_PI));
  542.     } else if (x < -XCMS_HALFPI) {
  543.     retval = -(_XcmsSine (x + XCMS_PI));
  544.     } else if (x > XCMS_FOURTHPI) {
  545.     retval = _XcmsCosine (XCMS_HALFPI - x);
  546.     } else if (x < -XCMS_FOURTHPI) {
  547.     retval = -(_XcmsCosine (XCMS_HALFPI + x));
  548.     } else if (x < XCMS_X6_UNDERFLOWS && x > -XCMS_X6_UNDERFLOWS) {
  549.     retval = x;
  550.     } else {
  551.     y = x / XCMS_FOURTHPI;
  552.     yt2 = y * y;
  553.     retval = y * (_XcmsPolynomial (3, sin_pcoeffs, yt2) / _XcmsPolynomial(3, sin_qcoeffs, yt2));
  554.     }
  555.     return(retval);
  556. }
  557.  
  558.  
  559. /*
  560.  *    NAME
  561.  *        _XcmsArcTangent
  562.  *
  563.  *    SYNOPSIS
  564.  */
  565. double
  566. _XcmsArcTangent(x)
  567.     double x;
  568. /*
  569.  *    DESCRIPTION
  570.  *        Computes the arctangent.
  571.  *        This is an implementation of the Gauss algorithm as
  572.  *        described in:
  573.  *            Forman S. Acton, Numerical Methods That Work,
  574.  *            New York, NY, Harper & Row, 1970.
  575.  *
  576.  *    RETURNS
  577.  *        Returns the arctangent 
  578.  */
  579. {
  580.     double ai, a1, bi, b1, l, d;
  581.     double maxerror;
  582.     int i;
  583.  
  584.     if (x == 0.0)  {
  585.     return (0.0);
  586.     }
  587.     if (x < 1.0) {
  588.     maxerror = x * XCMS_MAXERROR;
  589.     } else {
  590.     maxerror = XCMS_MAXERROR;
  591.     }
  592.     ai = _XcmsSquareRoot( 1.0 / (1.0 + (x * x)) );
  593.     bi = 1.0;
  594.     for (i = 0; i < XCMS_MAXITER; i++) {
  595.     a1 = (ai + bi) / 2.0;
  596.     b1 = _XcmsSquareRoot((a1 * bi));
  597.     if (a1 == b1)
  598.         break;
  599.     d = XCMS_FABS(a1 - b1);
  600.     if (d < maxerror)
  601.         break;
  602.     ai = a1;
  603.     bi = b1;
  604.     }
  605.  
  606.     l = ((a1 > b1) ? b1 : a1);
  607.  
  608.     a1 = _XcmsSquareRoot(1 + (x * x));
  609.     return (x / (a1 * l));
  610. }
  611.