diff options
Diffstat (limited to 'libm/double/mod2pi.c')
-rw-r--r-- | libm/double/mod2pi.c | 122 |
1 files changed, 0 insertions, 122 deletions
diff --git a/libm/double/mod2pi.c b/libm/double/mod2pi.c deleted file mode 100644 index 057954a9b..000000000 --- a/libm/double/mod2pi.c +++ /dev/null @@ -1,122 +0,0 @@ -/* Program to test range reduction of trigonometry functions - * - * -- Steve Moshier - */ - -#include <math.h> -#ifdef ANSIPROT -extern double floor ( double ); -extern double ldexp ( double, int ); -extern double sin ( double ); -#else -double floor(), ldexp(), sin(); -#endif - -#define TPI 6.283185307179586476925 - -main() -{ -char s[40]; -double a, n, t, x, y, z; -int lflg; - -x = TPI/4.0; -t = 1.0; - -loop: - -t = 2.0 * t; - -/* Stop testing at a point beyond which the integer part of - * x/2pi cannot be represented exactly by a double precision number. - * The library trigonometry functions will probably give up long before - * this point is reached. - */ -if( t > 1.0e16 ) - exit(0); - -/* Adjust the following to choose a nontrivial x - * where test function(x) has a slope of about 1 or more. - */ -x = TPI * t + 0.5; - -z = x; -lflg = 0; - -inlup: - -/* floor() returns the largest integer less than its argument. - * If you do not have this, or AINT(), then you may convert x/TPI - * to a long integer and then back to double; but in that case - * x will be limited to the largest value that will fit into a - * long integer. - */ -n = floor( z/TPI ); - -/* Carefully subtract 2 pi n from x. - * This is done by subtracting n * 2**k in such a way that there - * is no arithmetic cancellation error at any step. The k are the - * bits in the number 2 pi. - * - * If you do not have ldexp(), then you may multiply or - * divide n by an appropriate power of 2 after each step. - * For example: - * a = z - 4*n; - * a -= 2*n; - * n /= 4; - * a -= n; n/4 - * n /= 8; - * a -= n; n/32 - * etc. - * This will only work if division by a power of 2 is exact. - */ - -a = z - ldexp(n, 2); /* 4n */ -a -= ldexp( n, 1); /* 2n */ -a -= ldexp( n, -2 ); /* n/4 */ -a -= ldexp( n, -5 ); /* n/32 */ -a -= ldexp( n, -9 ); /* n/512 */ -a += ldexp( n, -15 ); /* add n/32768 */ -a -= ldexp( n, -17 ); /* n/131072 */ -a -= ldexp( n, -18 ); -a -= ldexp( n, -20 ); -a -= ldexp( n, -22 ); -a -= ldexp( n, -24 ); -a -= ldexp( n, -28 ); -a -= ldexp( n, -32 ); -a -= ldexp( n, -37 ); -a -= ldexp( n, -39 ); -a -= ldexp( n, -40 ); -a -= ldexp( n, -42 ); -a -= ldexp( n, -46 ); -a -= ldexp( n, -47 ); - -/* Subtract what is left of 2 pi n after all the above reductions. - */ -a -= 2.44929359829470635445e-16 * n; - -/* If the test is extended too far, it is possible - * to have chosen the wrong value of n. The following - * will fix that, but at some reduction in accuracy. - */ -if( (a > TPI) || (a < -1e-11) ) - { - z = a; - lflg += 1; - printf( "Warning! Reduction failed on first try.\n" ); - goto inlup; - } -if( a < 0.0 ) - { - printf( "Warning! Reduced value < 0\n" ); - a += TPI; - } - -/* Compute the test function at x and at a = x mod 2 pi. - */ -y = sin(x); -z = sin(a); -printf( "sin(%.15e) error = %.3e\n", x, y-z ); -goto loop; -} - |