Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : SIN.C

 
Output of file : SIN.C contained in archive : CEPHES22.ZIP
/* 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 );
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : SIN.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/