diff options
Diffstat (limited to 'libm/float/sindgf.c')
| -rw-r--r-- | libm/float/sindgf.c | 232 |
1 files changed, 232 insertions, 0 deletions
diff --git a/libm/float/sindgf.c b/libm/float/sindgf.c new file mode 100644 index 000000000..a3f5851c8 --- /dev/null +++ b/libm/float/sindgf.c @@ -0,0 +1,232 @@ +/* sindgf.c + * + * Circular sine of angle in degrees + * + * + * + * SYNOPSIS: + * + * float x, y, sindgf(); + * + * y = sindgf( 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 Q(x**2). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-3600 100,000 1.2e-7 3.0e-8 + * + * ERROR MESSAGES: + * + * message condition value returned + * sin total loss x > 2^24 0.0 + * + */ + +/* cosdgf.c + * + * Circular cosine of angle in degrees + * + * + * + * SYNOPSIS: + * + * float x, y, cosdgf(); + * + * y = cosdgf( 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 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 -8192,+8192 100,000 3.0e-7 3.0e-8 + */ + +/* +Cephes Math Library Release 2.2: June, 1992 +Copyright 1985, 1987, 1988, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +/* Single precision circular sine + * test interval: [-pi/4, +pi/4] + * trials: 10000 + * peak relative error: 6.8e-8 + * rms relative error: 2.6e-8 + */ +#include <math.h> + + +/*static float FOPI = 1.27323954473516;*/ + +extern float PIO4F; + +/* These are for a 24-bit significand: */ +static float T24M1 = 16777215.; + +static float PI180 = 0.0174532925199432957692; /* pi/180 */ + +float sindgf( float xx ) +{ +float x, y, z; +long j; +int sign; + +sign = 1; +x = xx; +if( xx < 0 ) + { + sign = -1; + x = -xx; + } +if( x > T24M1 ) + { + mtherr( "sindgf", TLOSS ); + return(0.0); + } +j = 0.022222222222222222222 * x; /* integer part of x/45 */ +y = j; +/* map zeros to origin */ +if( j & 1 ) + { + j += 1; + y += 1.0; + } +j &= 7; /* octant modulo 360 degrees */ +/* reflect in x axis */ +if( j > 3) + { + sign = -sign; + j -= 4; + } + +x = x - y * 45.0; +x *= PI180; /* multiply by pi/180 to convert to radians */ + +z = x * x; +if( (j==1) || (j==2) ) + { +/* + y = ((( 2.4462803166E-5 * z + - 1.3887580023E-3) * z + + 4.1666650433E-2) * z + - 4.9999999968E-1) * z + + 1.0; +*/ + +/* measured relative error in +/- pi/4 is 7.8e-8 */ + y = (( 2.443315711809948E-005 * z + - 1.388731625493765E-003) * z + + 4.166664568298827E-002) * z * z; + y -= 0.5 * z; + y += 1.0; + } +else + { +/* Theoretical relative error = 3.8e-9 in [-pi/4, +pi/4] */ + y = ((-1.9515295891E-4 * z + + 8.3321608736E-3) * z + - 1.6666654611E-1) * z * x; + y += x; + } + +if(sign < 0) + y = -y; +return( y); +} + + +/* Single precision circular cosine + * test interval: [-pi/4, +pi/4] + * trials: 10000 + * peak relative error: 8.3e-8 + * rms relative error: 2.2e-8 + */ + +float cosdgf( float xx ) +{ +register float x, y, z; +int j, sign; + +/* make argument positive */ +sign = 1; +x = xx; +if( x < 0 ) + x = -x; + +if( x > T24M1 ) + { + mtherr( "cosdgf", TLOSS ); + return(0.0); + } + +j = 0.02222222222222222222222 * x; /* integer part of x/PIO4 */ +y = j; +/* integer and fractional part modulo one octant */ +if( j & 1 ) /* map zeros to origin */ + { + j += 1; + y += 1.0; + } +j &= 7; +if( j > 3) + { + j -=4; + sign = -sign; + } + +if( j > 1 ) + sign = -sign; + +x = x - y * 45.0; /* x mod 45 degrees */ +x *= PI180; /* multiply by pi/180 to convert to radians */ + +z = x * x; + +if( (j==1) || (j==2) ) + { + y = (((-1.9515295891E-4 * z + + 8.3321608736E-3) * z + - 1.6666654611E-1) * z * x) + + x; + } +else + { + y = (( 2.443315711809948E-005 * z + - 1.388731625493765E-003) * z + + 4.166664568298827E-002) * z * z; + y -= 0.5 * z; + y += 1.0; + } +if(sign < 0) + y = -y; +return( y ); +} + |