309 lines
6.0 KiB
C
309 lines
6.0 KiB
C
|
/* sindg.c
|
|||
|
|
*
|
|||
|
|
* Circular sine of angle in degrees
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* double x, y, sindg();
|
|||
|
|
*
|
|||
|
|
* y = sindg( x );
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* Range reduction is into intervals of 45 degrees.
|
|||
|
|
*
|
|||
|
|
* 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 P(x**2).
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* ACCURACY:
|
|||
|
|
*
|
|||
|
|
* Relative error:
|
|||
|
|
* arithmetic domain # trials peak rms
|
|||
|
|
* DEC +-1000 3100 3.3e-17 9.0e-18
|
|||
|
|
* IEEE +-1000 30000 2.3e-16 5.6e-17
|
|||
|
|
*
|
|||
|
|
* ERROR MESSAGES:
|
|||
|
|
*
|
|||
|
|
* message condition value returned
|
|||
|
|
* sindg total loss x > 8.0e14 (DEC) 0.0
|
|||
|
|
* x > 1.0e14 (IEEE)
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
/* cosdg.c
|
|||
|
|
*
|
|||
|
|
* Circular cosine of angle in degrees
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* double x, y, cosdg();
|
|||
|
|
*
|
|||
|
|
* y = cosdg( x );
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* Range reduction is into intervals of 45 degrees.
|
|||
|
|
*
|
|||
|
|
* Two polynomial approximating functions are employed.
|
|||
|
|
* Between 0 and pi/4 the cosine is approximated by
|
|||
|
|
* 1 - x**2 P(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
|
|||
|
|
* DEC +-1000 3400 3.5e-17 9.1e-18
|
|||
|
|
* IEEE +-1000 30000 2.1e-16 5.7e-17
|
|||
|
|
* See also sin().
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/* Cephes Math Library Release 2.0: April, 1987
|
|||
|
|
* Copyright 1985, 1987 by Stephen L. Moshier
|
|||
|
|
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */
|
|||
|
|
|
|||
|
|
#include "mconf.h"
|
|||
|
|
|
|||
|
|
#ifdef UNK
|
|||
|
|
const static double sincof[] = {
|
|||
|
|
1.58962301572218447952E-10,
|
|||
|
|
-2.50507477628503540135E-8,
|
|||
|
|
2.75573136213856773549E-6,
|
|||
|
|
-1.98412698295895384658E-4,
|
|||
|
|
8.33333333332211858862E-3,
|
|||
|
|
-1.66666666666666307295E-1
|
|||
|
|
};
|
|||
|
|
const static double coscof[] = {
|
|||
|
|
1.13678171382044553091E-11,
|
|||
|
|
-2.08758833757683644217E-9,
|
|||
|
|
2.75573155429816611547E-7,
|
|||
|
|
-2.48015872936186303776E-5,
|
|||
|
|
1.38888888888806666760E-3,
|
|||
|
|
-4.16666666666666348141E-2,
|
|||
|
|
4.99999999999999999798E-1
|
|||
|
|
};
|
|||
|
|
const static double PI180 = 1.74532925199432957692E-2; /* pi/180 */
|
|||
|
|
const static double lossth = 1.0e14;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef DEC
|
|||
|
|
static unsigned short sincof[] = {
|
|||
|
|
0030056,0143750,0177170,0073013,
|
|||
|
|
0131727,0027455,0044510,0132205,
|
|||
|
|
0033470,0167432,0131752,0042263,
|
|||
|
|
0135120,0006400,0146776,0174027,
|
|||
|
|
0036410,0104210,0104207,0137202,
|
|||
|
|
0137452,0125252,0125252,0125103
|
|||
|
|
};
|
|||
|
|
static unsigned short coscof[] = {
|
|||
|
|
0027107,0176030,0153315,0110312,
|
|||
|
|
0131017,0072476,0007450,0123243,
|
|||
|
|
0032623,0171174,0070066,0146445,
|
|||
|
|
0134320,0006400,0147355,0163313,
|
|||
|
|
0035666,0005540,0133012,0165067,
|
|||
|
|
0137052,0125252,0125252,0125206,
|
|||
|
|
0040000,0000000,0000000,0000000
|
|||
|
|
};
|
|||
|
|
static unsigned short P1[] = {0036616,0175065,0011224,0164711};
|
|||
|
|
#define PI180 *(double *)P1
|
|||
|
|
const static double lossth = 8.0e14;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef IBMPC
|
|||
|
|
static unsigned short sincof[] = {
|
|||
|
|
0x0ec1,0x1fcf,0xd8fd,0x3de5,
|
|||
|
|
0x1691,0xa929,0xe5e5,0xbe5a,
|
|||
|
|
0x4896,0x567d,0x1de3,0x3ec7,
|
|||
|
|
0xdf03,0x19bf,0x01a0,0xbf2a,
|
|||
|
|
0xf7d0,0x1110,0x1111,0x3f81,
|
|||
|
|
0x5548,0x5555,0x5555,0xbfc5
|
|||
|
|
};
|
|||
|
|
static unsigned short coscof[] = {
|
|||
|
|
0xb219,0x1ad9,0xff83,0x3da8,
|
|||
|
|
0x14d4,0xc1e5,0xeea7,0xbe21,
|
|||
|
|
0xd9a5,0x8e06,0x7e4f,0x3e92,
|
|||
|
|
0xbcd9,0x19dd,0x01a0,0xbefa,
|
|||
|
|
0x5d47,0x16c1,0xc16c,0x3f56,
|
|||
|
|
0x5551,0x5555,0x5555,0xbfa5,
|
|||
|
|
0x0000,0x0000,0x0000,0x3fe0
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
static unsigned short P1[] = {0x9d39,0xa252,0xdf46,0x3f91};
|
|||
|
|
#define PI180 *(double *)P1
|
|||
|
|
const static double lossth = 1.0e14;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef MIEEE
|
|||
|
|
static unsigned short sincof[] = {
|
|||
|
|
0x3de5,0xd8fd,0x1fcf,0x0ec1,
|
|||
|
|
0xbe5a,0xe5e5,0xa929,0x1691,
|
|||
|
|
0x3ec7,0x1de3,0x567d,0x4896,
|
|||
|
|
0xbf2a,0x01a0,0x19bf,0xdf03,
|
|||
|
|
0x3f81,0x1111,0x1110,0xf7d0,
|
|||
|
|
0xbfc5,0x5555,0x5555,0x5548
|
|||
|
|
};
|
|||
|
|
static unsigned short coscof[] = {
|
|||
|
|
0x3da8,0xff83,0x1ad9,0xb219,
|
|||
|
|
0xbe21,0xeea7,0xc1e5,0x14d4,
|
|||
|
|
0x3e92,0x7e4f,0x8e06,0xd9a5,
|
|||
|
|
0xbefa,0x01a0,0x19dd,0xbcd9,
|
|||
|
|
0x3f56,0xc16c,0x16c1,0x5d47,
|
|||
|
|
0xbfa5,0x5555,0x5555,0x5551,
|
|||
|
|
0x3fe0,0x0000,0x0000,0x0000
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
static unsigned short P1[] = {
|
|||
|
|
0x3f91,0xdf46,0xa252,0x9d39
|
|||
|
|
};
|
|||
|
|
#define PI180 *(double *)P1
|
|||
|
|
const static double lossth = 1.0e14;
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
#ifdef ANSIPROT
|
|||
|
|
extern double polevl ( double, void *, int );
|
|||
|
|
extern double floor ( double );
|
|||
|
|
extern double ldexp ( double, int );
|
|||
|
|
#else
|
|||
|
|
double polevl(), floor(), ldexp();
|
|||
|
|
#endif
|
|||
|
|
extern double PIO4;
|
|||
|
|
|
|||
|
|
double sindg(x)
|
|||
|
|
double x;
|
|||
|
|
{
|
|||
|
|
double y, z, zz;
|
|||
|
|
int j, sign;
|
|||
|
|
|
|||
|
|
/* make argument positive but save the sign */
|
|||
|
|
sign = 1;
|
|||
|
|
if( x < 0 )
|
|||
|
|
{
|
|||
|
|
x = -x;
|
|||
|
|
sign = -1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if( x > lossth )
|
|||
|
|
{
|
|||
|
|
mtherr( "sindg", TLOSS );
|
|||
|
|
return(0.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
y = floor( x/45.0 ); /* 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;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
z = x - y * 45.0; /* x mod 45 degrees */
|
|||
|
|
z *= PI180; /* multiply by pi/180 to convert to radians */
|
|||
|
|
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 cosdg(x)
|
|||
|
|
double x;
|
|||
|
|
{
|
|||
|
|
double y, z, zz;
|
|||
|
|
int j, sign;
|
|||
|
|
|
|||
|
|
/* make argument positive */
|
|||
|
|
sign = 1;
|
|||
|
|
if( x < 0 )
|
|||
|
|
x = -x;
|
|||
|
|
|
|||
|
|
if( x > lossth )
|
|||
|
|
{
|
|||
|
|
mtherr( "cosdg", TLOSS );
|
|||
|
|
return(0.0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
y = floor( x/45.0 );
|
|||
|
|
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 */
|
|||
|
|
j = z;
|
|||
|
|
if( j & 1 ) /* map zeros to origin */
|
|||
|
|
{
|
|||
|
|
j += 1;
|
|||
|
|
y += 1.0;
|
|||
|
|
}
|
|||
|
|
j = j & 07;
|
|||
|
|
if( j > 3)
|
|||
|
|
{
|
|||
|
|
j -=4;
|
|||
|
|
sign = -sign;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if( j > 1 )
|
|||
|
|
sign = -sign;
|
|||
|
|
|
|||
|
|
z = x - y * 45.0; /* x mod 45 degrees */
|
|||
|
|
z *= PI180; /* multiply by pi/180 to convert to radians */
|
|||
|
|
|
|||
|
|
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);
|
|||
|
|
}
|