home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / mod2pi.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.5 KB  |  115 lines
  1. /* Program to test range reduction of trigonometry functions
  2.  *
  3.  * -- Steve Moshier
  4.  */
  5.  
  6. #define TPI 6.283185307179586476925
  7.  
  8. main()
  9. {
  10. char s[40];
  11. double a, n, t, x, y, z;
  12. int lflg;
  13. double floor(), ldexp(), sin();
  14.  
  15. x = TPI/4.0;
  16. t = 1.0;
  17.  
  18. loop:
  19.  
  20. t = 2.0 * t;
  21.  
  22. /* Stop testing at a point beyond which the integer part of
  23.  * x/2pi cannot be represented exactly by a double precision number.
  24.  * The library trigonometry functions will probably give up long before
  25.  * this point is reached.
  26.  */
  27. if( t > 1.0e16 )
  28.     exit(0);
  29.  
  30. /* Adjust the following to choose a nontrivial x
  31.  * where test function(x) has a slope of about 1 or more.
  32.  */
  33. x = TPI * t  + 0.5;
  34.  
  35. z = x;
  36. lflg = 0;
  37.  
  38. inlup:
  39.  
  40. /* floor() returns the largest integer less than its argument.
  41.  * If you do not have this, or AINT(), then you may convert x/TPI
  42.  * to a long integer and then back to double; but in that case
  43.  * x will be limited to the largest value that will fit into a
  44.  * long integer.
  45.  */
  46. n = floor( z/TPI );
  47.  
  48. /* Carefully subtract 2 pi n from x.
  49.  * This is done by subtracting n * 2**k in such a way that there
  50.  * is no arithmetic cancellation error at any step.  The k are the
  51.  * bits in the number 2 pi.
  52.  *
  53.  * If you do not have ldexp(), then you may multiply or
  54.  * divide n by an appropriate power of 2 after each step.
  55.  * For example:
  56.  *  a = z - 4*n;
  57.  *  a -= 2*n;
  58.  *  n /= 4;
  59.  *  a -= n;   n/4
  60.  *  n /= 8;
  61.  *  a -= n;   n/32
  62.  * etc.
  63.  * This will only work if division by a power of 2 is exact.
  64.  */
  65.  
  66. a = z - ldexp(n, 2);    /* 4n */
  67. a -= ldexp( n, 1);    /* 2n */
  68. a -= ldexp( n, -2 );    /* n/4 */
  69. a -= ldexp( n, -5 );    /* n/32 */
  70. a -= ldexp( n, -9 );    /* n/512 */
  71. a += ldexp( n, -15 );    /* add n/32768 */
  72. a -= ldexp( n, -17 );    /* n/131072 */
  73. a -= ldexp( n, -18 );
  74. a -= ldexp( n, -20 );
  75. a -= ldexp( n, -22 );
  76. a -= ldexp( n, -24 );
  77. a -= ldexp( n, -28 );
  78. a -= ldexp( n, -32 );
  79. a -= ldexp( n, -37 );
  80. a -= ldexp( n, -39 );
  81. a -= ldexp( n, -40 );
  82. a -= ldexp( n, -42 );
  83. a -= ldexp( n, -46 );
  84. a -= ldexp( n, -47 );
  85.  
  86. /* Subtract what is left of 2 pi n after all the above reductions.
  87.  */
  88. a -= 2.44929359829470635445e-16 * n;
  89.  
  90. /* If the test is extended too far, it is possible
  91.  * to have chosen the wrong value of n.  The following
  92.  * will fix that, but at some reduction in accuracy.
  93.  */
  94. if( (a > TPI) || (a < -1e-11) )
  95.     {
  96.     z = a;
  97.     lflg += 1;
  98.     printf( "Warning! Reduction failed on first try.\n" );
  99.     goto inlup;
  100.     }
  101. if( a < 0.0 )
  102.     {
  103.     printf( "Warning! Reduced value < 0\n" );
  104.     a += TPI;
  105.     }
  106.  
  107. /* Compute the test function at x and at a = x mod 2 pi.
  108.  */
  109. y = sin(x);
  110. z = sin(a);
  111. printf( "sin(%.15e) error = %.3e\n", x, y-z );
  112. goto loop;
  113. }
  114.  
  115.