diff options
Diffstat (limited to 'libm/float/sinf.c')
| -rw-r--r-- | libm/float/sinf.c | 283 |
1 files changed, 283 insertions, 0 deletions
diff --git a/libm/float/sinf.c b/libm/float/sinf.c new file mode 100644 index 000000000..2f1bb45b8 --- /dev/null +++ b/libm/float/sinf.c @@ -0,0 +1,283 @@ +/* sinf.c + * + * Circular sine + * + * + * + * SYNOPSIS: + * + * float x, y, sinf(); + * + * y = sinf( 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 + * IEEE -4096,+4096 100,000 1.2e-7 3.0e-8 + * IEEE -8192,+8192 100,000 3.0e-7 3.0e-8 + * + * ERROR MESSAGES: + * + * message condition value returned + * sin total loss x > 2^24 0.0 + * + * Partial loss of accuracy begins to occur at x = 2^13 + * = 8192. Results may be meaningless for x >= 2^24 + * The routine as implemented flags a TLOSS error + * for x >= 2^24 and returns 0.0. + */ + +/* cosf.c + * + * Circular cosine + * + * + * + * SYNOPSIS: + * + * float x, y, cosf(); + * + * y = cosf( 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 -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; +/* Note, these constants are for a 32-bit significand: */ +/* +static float DP1 = 0.7853851318359375; +static float DP2 = 1.30315311253070831298828125e-5; +static float DP3 = 3.03855025325309630e-11; +static float lossth = 65536.; +*/ + +/* These are for a 24-bit significand: */ +static float DP1 = 0.78515625; +static float DP2 = 2.4187564849853515625e-4; +static float DP3 = 3.77489497744594108e-8; +static float lossth = 8192.; +static float T24M1 = 16777215.; + +static float sincof[] = { +-1.9515295891E-4, + 8.3321608736E-3, +-1.6666654611E-1 +}; +static float coscof[] = { + 2.443315711809948E-005, +-1.388731625493765E-003, + 4.166664568298827E-002 +}; + +float sinf( float xx ) +{ +float *p; +float x, y, z; +register unsigned long j; +register int sign; + +sign = 1; +x = xx; +if( xx < 0 ) + { + sign = -1; + x = -xx; + } +if( x > T24M1 ) + { + mtherr( "sinf", TLOSS ); + return(0.0); + } +j = FOPI * x; /* integer part of x/(PI/4) */ +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; + } + +if( x > lossth ) + { + mtherr( "sinf", PLOSS ); + x = x - y * PIO4F; + } +else + { +/* Extended precision modular arithmetic */ + x = ((x - y * DP1) - y * DP2) - y * DP3; + } +/*einits();*/ +z = x * x; +if( (j==1) || (j==2) ) + { +/* measured relative error in +/- pi/4 is 7.8e-8 */ +/* + y = (( 2.443315711809948E-005 * z + - 1.388731625493765E-003) * z + + 4.166664568298827E-002) * z * z; +*/ + p = coscof; + y = *p++; + y = y * z + *p++; + y = y * z + *p++; + y *= 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; +*/ + p = sincof; + y = *p++; + y = y * z + *p++; + y = y * z + *p++; + y *= z * x; + y += x; + } +/*einitd();*/ +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 cosf( float xx ) +{ +float x, y, z; +int j, sign; + +/* make argument positive */ +sign = 1; +x = xx; +if( x < 0 ) + x = -x; + +if( x > T24M1 ) + { + mtherr( "cosf", TLOSS ); + return(0.0); + } + +j = FOPI * 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; + +if( x > lossth ) + { + mtherr( "cosf", PLOSS ); + x = x - y * PIO4F; + } +else +/* Extended precision modular arithmetic */ + x = ((x - y * DP1) - y * DP2) - y * DP3; + +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 ); +} + |