summaryrefslogtreecommitdiff
path: root/libm/double/sin.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/sin.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/sin.c')
-rw-r--r--libm/double/sin.c387
1 files changed, 0 insertions, 387 deletions
diff --git a/libm/double/sin.c b/libm/double/sin.c
deleted file mode 100644
index 24746d79d..000000000
--- a/libm/double/sin.c
+++ /dev/null
@@ -1,387 +0,0 @@
-/* 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 150000 3.0e-17 7.8e-18
- * IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-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 130000 2.1e-16 5.4e-17
- * DEC 0,+1.07e9 17000 3.0e-17 7.2e-18
- */
-
-/* sin.c */
-
-/*
-Cephes Math Library Release 2.8: June, 2000
-Copyright 1985, 1995, 2000 by Stephen L. Moshier
-*/
-
-#include <math.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 coscof[6] = {
--1.13585365213876817300E-11,
- 2.08757008419747316778E-9,
--2.75573141792967388112E-7,
- 2.48015872888517045348E-5,
--1.38888888888730564116E-3,
- 4.16666666666665929218E-2,
-};
-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 unsigned 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,
-};
-static unsigned short coscof[24] = {
-0127107,0151115,0002060,0152325,
-0031017,0072353,0155161,0174053,
-0132623,0171173,0172542,0057056,
-0034320,0006400,0147102,0023652,
-0135666,0005540,0133012,0076213,
-0037052,0125252,0125252,0125126,
-};
-/* 7.853981629014015197753906250000E-1 */
-static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
-/* 4.960467869796758577649598009884E-10 */
-static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
-/* 2.860594363054915898381331279295E-18 */
-static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
-#define DP1 *(double *)P1
-#define DP2 *(double *)P2
-#define DP3 *(double *)P3
-#endif
-
-#ifdef IBMPC
-static unsigned 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,
-};
-static unsigned short coscof[24] = {
-0x1a9b,0xa086,0xfa49,0xbda8,
-0x3f05,0x7b4e,0xee9d,0x3e21,
-0x4bc6,0x7eac,0x7e4f,0xbe92,
-0x44f5,0x19c8,0x01a0,0x3efa,
-0x4f91,0x16c1,0xc16c,0xbf56,
-0x554b,0x5555,0x5555,0x3fa5,
-};
-/*
- 7.85398125648498535156E-1,
- 3.77489470793079817668E-8,
- 2.69515142907905952645E-15,
-*/
-static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
-static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
-static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
-#define DP1 *(double *)P1
-#define DP2 *(double *)P2
-#define DP3 *(double *)P3
-#endif
-
-#ifdef MIEEE
-static unsigned 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 unsigned short coscof[24] = {
-0xbda8,0xfa49,0xa086,0x1a9b,
-0x3e21,0xee9d,0x7b4e,0x3f05,
-0xbe92,0x7e4f,0x7eac,0x4bc6,
-0x3efa,0x01a0,0x19c8,0x44f5,
-0xbf56,0xc16c,0x16c1,0x4f91,
-0x3fa5,0x5555,0x5555,0x554b,
-};
-static unsigned short P1[] = {0x3fe9,0x21fb,0x4000,0x0000};
-static unsigned short P2[] = {0x3e64,0x442d,0x0000,0x0000};
-static unsigned short P3[] = {0x3ce8,0x4698,0x98cc,0x5170};
-#define DP1 *(double *)P1
-#define DP2 *(double *)P2
-#define DP3 *(double *)P3
-#endif
-
-#ifdef ANSIPROT
-extern double polevl ( double, void *, int );
-extern double p1evl ( double, void *, int );
-extern double floor ( double );
-extern double ldexp ( double, int );
-extern int isnan ( double );
-extern int isfinite ( double );
-#else
-double polevl(), floor(), ldexp();
-int isnan(), isfinite();
-#endif
-extern double PIO4;
-static double lossth = 1.073741824e9;
-#ifdef NANS
-extern double NAN;
-#endif
-#ifdef INFINITIES
-extern double INFINITY;
-#endif
-
-
-double sin(x)
-double x;
-{
-double y, z, zz;
-int j, sign;
-
-#ifdef MINUSZERO
-if( x == 0.0 )
- return(x);
-#endif
-#ifdef NANS
-if( isnan(x) )
- return(x);
-if( !isfinite(x) )
- {
- mtherr( "sin", DOMAIN );
- return(NAN);
- }
-#endif
-/* 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 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
- }
-else
- {
-/* y = z + z * (zz * polevl( zz, sincof, 5 ));*/
- y = z + z * z * z * 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;
-
-#ifdef NANS
-if( isnan(x) )
- return(x);
-if( !isfinite(x) )
- {
- mtherr( "cos", DOMAIN );
- return(NAN);
- }
-#endif
-
-/* 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 ));*/
- y = z + z * z * z * polevl( zz, sincof, 5 );
- }
-else
- {
- y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
- }
-
-if(sign < 0)
- y = -y;
-
-return(y);
-}
-
-
-
-
-
-/* Degrees, minutes, seconds to radians: */
-
-/* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
-#ifdef DEC
-static unsigned 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 );
-}