/* sin.c * * Circular sine * * * * SYNOPSIS: * * double x, y, sin(); * * y = sin( x ); * * * * DESCRIPTION: * * Range reduction is into intervals of pi/4. The reduction * error is nearly eliminated by contriving an extended precision * modular arithmetic. * * Two polynomial approximating functions are employed. * Between 0 and pi/4 the sine is approximated by * x + x**3 P(x**2). * Between pi/4 and pi/2 the cosine is represented as * 1 - x**2 Q(x**2). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC 0, 10 20000 2.5e-17 7.1e-18 * DEC 0, 1.07e9 10000 2.8e-17 7.2e-18 * IEEE -1.07e9,+1.07e9 40000 2.2e-16 5.6e-17 * * ERROR MESSAGES: * * message condition value returned * sin total loss x > 1.073741824e9 0.0 * * Partial loss of accuracy begins to occur at x = 2**30 * = 1.074e9. The loss is not gradual, but jumps suddenly to * about 1 part in 10e7. Results may be meaningless for * x > 2**49 = 5.6e14. The routine as implemented flags a * TLOSS error for x > 2**30 and returns 0.0. */ /* cos.c * * Circular cosine * * * * SYNOPSIS: * * double x, y, cos(); * * y = cos( x ); * * * * DESCRIPTION: * * Range reduction is into intervals of pi/4. The reduction * error is nearly eliminated by contriving an extended precision * modular arithmetic. * * Two polynomial approximating functions are employed. * Between 0 and pi/4 the cosine is approximated by * 1 - x**2 Q(x**2). * Between pi/4 and pi/2 the sine is represented as * x + x**3 P(x**2). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -1.07e9,+1.07e9 30000 2.0e-16 5.6e-17 * DEC 0,+1.07e9 17000 3.0e-17 7.2e-18 */ /* sin.c */ /* Cephes Math Library Release 2.1: October, 1988 Copyright 1985, 1987, 1988 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ #include "mconf.h" #ifdef UNK static double sincof[] = { 1.58962301576546568060E-10, -2.50507477628578072866E-8, 2.75573136213857245213E-6, -1.98412698295895385996E-4, 8.33333333332211858878E-3, -1.66666666666666307295E-1, }; static double DP1 = 7.85398125648498535156E-1; static double DP2 = 3.77489470793079817668E-8; static double DP3 = 2.69515142907905952645E-15; static double lossth = 1.073741824e9; #endif #ifdef DEC static short sincof[] = { 0030056,0143750,0177214,0163153, 0131727,0027455,0044510,0175352, 0033470,0167432,0131752,0042414, 0135120,0006400,0146776,0174027, 0036410,0104210,0104207,0137202, 0137452,0125252,0125252,0125103, }; /* 7.853981629014015197753906250000E-1 */ static short P1[] = {0040111,0007732,0120000,0000000,}; /* 4.960467869796758577649598009884E-10 */ static short P2[] = {0030410,0055060,0100000,0000000,}; /* 2.860594363054915898381331279295E-18 */ static short P3[] = {0021523,0011431,0105056,0001560,}; #define DP1 *(double *)P1 #define DP2 *(double *)P2 #define DP3 *(double *)P3 static double lossth = 1.073741824e9; #endif #ifdef IBMPC static short sincof[] = { 0x9ccd,0x1fd1,0xd8fd,0x3de5, 0x1f5d,0xa929,0xe5e5,0xbe5a, 0x48a1,0x567d,0x1de3,0x3ec7, 0xdf03,0x19bf,0x01a0,0xbf2a, 0xf7d0,0x1110,0x1111,0x3f81, 0x5548,0x5555,0x5555,0xbfc5, }; /* 7.85398125648498535156E-1, 3.77489470793079817668E-8, 2.69515142907905952645E-15, */ static short P1[] = {0x0000,0x4000,0x21fb,0x3fe9}; static short P2[] = {0x0000,0x0000,0x442d,0x3e64}; static short P3[] = {0x5170,0x98cc,0x4698,0x3ce8}; #define DP1 *(double *)P1 #define DP2 *(double *)P2 #define DP3 *(double *)P3 static double lossth = 1.073741824e9; #endif #ifdef MIEEE static short sincof[] = { 0x3de5,0xd8fd,0x1fd1,0x9ccd, 0xbe5a,0xe5e5,0xa929,0x1f5d, 0x3ec7,0x1de3,0x567d,0x48a1, 0xbf2a,0x01a0,0x19bf,0xdf03, 0x3f81,0x1111,0x1110,0xf7d0, 0xbfc5,0x5555,0x5555,0x5548, }; static short P1[] = { 0x3fe9,0x21fb,0x4000,0x0000 }; static short P2[] = { 0x3e64,0x442d,0x0000,0x0000 }; static short P3[] = { 0x3ce8,0x4698,0x98cc,0x5170 }; #define DP1 *(double *)P1 #define DP2 *(double *)P2 #define DP3 *(double *)P3 static double lossth = 1.073741824e9; #endif #ifdef UNK static double coscof[7] = { 1.13678171380010505367E-11, -2.08758833757646780967E-9, 2.75573155429816368675E-7, -2.48015872936186303093E-5, 1.38888888888806666760E-3, -4.16666666666666348141E-2, 4.99999999999999999798E-1, }; #endif #ifdef DEC static short coscof[28] = { 0027107,0176030,0153276,0031137, 0131017,0072476,0007450,0105310, 0032623,0171174,0070066,0146400, 0134320,0006400,0147355,0163313, 0035666,0005540,0133012,0165067, 0137052,0125252,0125252,0125206, 0040000,0000000,0000000,0000000, }; #endif #ifdef IBMPC static short coscof[28] = { 0xc64c,0x1ad7,0xff83,0x3da8, 0x1159,0xc1e5,0xeea7,0xbe21, 0xd9a0,0x8e06,0x7e4f,0x3e92, 0xbcd9,0x19dd,0x01a0,0xbefa, 0x5d47,0x16c1,0xc16c,0x3f56, 0x5551,0x5555,0x5555,0xbfa5, 0x0000,0x0000,0x0000,0x3fe0, }; #endif #ifdef MIEEE static short coscof[28] = { 0x3da8,0xff83,0x1ad7,0xc64c, 0xbe21,0xeea7,0xc1e5,0x1159, 0x3e92,0x7e4f,0x8e06,0xd9a0, 0xbefa,0x01a0,0x19dd,0xbcd9, 0x3f56,0xc16c,0x16c1,0x5d47, 0xbfa5,0x5555,0x5555,0x5551, 0x3fe0,0x0000,0x0000,0x0000, }; #endif extern double PIO4; double sin(x) double x; { double y, z, zz; int rflg, j, sign; double polevl(), floor(), ldexp(); /* make argument positive but save the sign */ sign = 1; if( x < 0 ) { x = -x; sign = -1; } if( x > lossth ) { mtherr( "sin", TLOSS ); return(0.0); } y = floor( x/PIO4 ); /* integer part of x/PIO4 */ /* strip high bits of integer part to prevent integer overflow */ z = ldexp( y, -4 ); z = floor(z); /* integer part of y/8 */ z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */ j = z; /* convert to integer for tests on the phase angle */ /* map zeros to origin */ if( j & 1 ) { j += 1; y += 1.0; } j = j & 07; /* octant modulo 360 degrees */ /* reflect in x axis */ if( j > 3) { sign = -sign; j -= 4; } /* Extended precision modular arithmetic */ z = ((x - y * DP1) - y * DP2) - y * DP3; zz = z * z; if( (j==1) || (j==2) ) { y = 1.0 - zz * polevl( zz, coscof, 6 ); } else { y = z + z * (zz * polevl( zz, sincof, 5 )); } if(sign < 0) y = -y; return(y); } double cos(x) double x; { double y, z, zz; long i; int j, sign, refl; double polevl(), floor(), ldexp(); /* make argument positive */ sign = 1; if( x < 0 ) x = -x; if( x > lossth ) { mtherr( "cos", TLOSS ); return(0.0); } y = floor( x/PIO4 ); z = ldexp( y, -4 ); z = floor(z); /* integer part of y/8 */ z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */ /* integer and fractional part modulo one octant */ i = z; if( i & 1 ) /* map zeros to origin */ { i += 1; y += 1.0; } j = i & 07; if( j > 3) { j -=4; sign = -sign; } if( j > 1 ) sign = -sign; /* Extended precision modular arithmetic */ z = ((x - y * DP1) - y * DP2) - y * DP3; zz = z * z; if( (j==1) || (j==2) ) { y = z + z * (zz * polevl( zz, sincof, 5 )); } else { y = 1.0 - zz * polevl( zz, coscof, 6 ); } if(sign < 0) y = -y; return(y); } /* Degrees, minutes, seconds to radians: */ /* 1 arc second, in radians = 4.8481368110953599358991410e-5 */ #ifdef DEC static short P648[] = {034513,054170,0176773,0116043,}; #define P64800 *(double *)P648 #else static double P64800 = 4.8481368110953599358991410e-5; #endif double radian(d,m,s) double d,m,s; { return( ((d*60.0 + m)*60.0 + s)*P64800 ); }