Output of file : MOD2PI.C contained in archive : CEPHES22.ZIP
/* Program to test range reduction of trigonometry functions
*
* -- Steve Moshier
*/

#define TPI 6.283185307179586476925

main()
{
char s[40];
double a, n, t, x, y, z;
int lflg;
double floor(), ldexp(), sin();

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

### 3 Responses to “Category : C Source CodeArchive   : CEPHES22.ZIPFilename : MOD2PI.C”

1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/