Compute sin(pi*x) without actually doing the pi*x multiplication.
sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is the precision of x. The new argument reduction is an optimization compared to the old code, and it removes a chunk of dead code. Accuracy tests in the intervals (-21,-20), (-20,-19), ... (-1,0) show no differences between the old and new code. Obtained from: bde
This commit is contained in:
parent
f0d2731dd8
commit
795b92049d
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=270893
@ -156,37 +156,35 @@ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
|
||||
|
||||
static const double zero= 0.00000000000000000000e+00;
|
||||
|
||||
static double sin_pi(double x)
|
||||
/*
|
||||
* Compute sin(pi*x) without actually doing the pi*x multiplication.
|
||||
* sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
|
||||
* the precision of x.
|
||||
*/
|
||||
static double
|
||||
sin_pi(double x)
|
||||
{
|
||||
volatile double vz;
|
||||
double y,z;
|
||||
int n,ix;
|
||||
int n;
|
||||
|
||||
GET_HIGH_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
y = -x;
|
||||
|
||||
if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
|
||||
y = -x; /* x is assume negative */
|
||||
vz = y+0x1p52; /* depend on 0 <= y < 0x1p52 */
|
||||
z = vz-0x1p52; /* rint(y) for the above range */
|
||||
if (z == y)
|
||||
return (zero);
|
||||
|
||||
vz = y+0x1p50;
|
||||
GET_LOW_WORD(n,vz); /* bits for rounded y (units 0.25) */
|
||||
z = vz-0x1p50; /* y rounded to a multiple of 0.25 */
|
||||
if (z > y) {
|
||||
z -= 0.25; /* adjust to round down */
|
||||
n--;
|
||||
}
|
||||
n &= 7; /* octant of y mod 2 */
|
||||
y = y - z + n * 0.25; /* y mod 2 */
|
||||
|
||||
/*
|
||||
* argument reduction, make sure inexact flag not raised if input
|
||||
* is an integer
|
||||
*/
|
||||
z = floor(y);
|
||||
if(z!=y) { /* inexact anyway */
|
||||
y *= 0.5;
|
||||
y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */
|
||||
n = (int) (y*4.0);
|
||||
} else {
|
||||
if(ix>=0x43400000) {
|
||||
y = zero; n = 0; /* y must be even */
|
||||
} else {
|
||||
if(ix<0x43300000) z = y+two52; /* exact */
|
||||
GET_LOW_WORD(n,z);
|
||||
n &= 1;
|
||||
y = n;
|
||||
n<<= 2;
|
||||
}
|
||||
}
|
||||
switch (n) {
|
||||
case 0: y = __kernel_sin(pi*y,zero,0); break;
|
||||
case 1:
|
||||
|
@ -89,37 +89,35 @@ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
|
||||
|
||||
static const float zero= 0.0000000000e+00;
|
||||
|
||||
static float sin_pif(float x)
|
||||
/*
|
||||
* Compute sin(pi*x) without actually doing the pi*x multiplication.
|
||||
* sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
|
||||
* the precision of x.
|
||||
*/
|
||||
static float
|
||||
sin_pi(float x)
|
||||
{
|
||||
volatile float vz;
|
||||
float y,z;
|
||||
int n,ix;
|
||||
int n;
|
||||
|
||||
GET_FLOAT_WORD(ix,x);
|
||||
ix &= 0x7fffffff;
|
||||
y = -x;
|
||||
|
||||
if(ix<0x3e800000) return __kernel_sindf(pi*x);
|
||||
y = -x; /* x is assume negative */
|
||||
vz = y+0x1p23F; /* depend on 0 <= y < 0x1p23 */
|
||||
z = vz-0x1p23F; /* rintf(y) for the above range */
|
||||
if (z == y)
|
||||
return (zero);
|
||||
|
||||
vz = y+0x1p21F;
|
||||
GET_FLOAT_WORD(n,vz); /* bits for rounded y (units 0.25) */
|
||||
z = vz-0x1p21F; /* y rounded to a multiple of 0.25 */
|
||||
if (z > y) {
|
||||
z -= 0.25F; /* adjust to round down */
|
||||
n--;
|
||||
}
|
||||
n &= 7; /* octant of y mod 2 */
|
||||
y = y - z + n * 0.25F; /* y mod 2 */
|
||||
|
||||
/*
|
||||
* argument reduction, make sure inexact flag not raised if input
|
||||
* is an integer
|
||||
*/
|
||||
z = floorf(y);
|
||||
if(z!=y) { /* inexact anyway */
|
||||
y *= (float)0.5;
|
||||
y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */
|
||||
n = (int) (y*(float)4.0);
|
||||
} else {
|
||||
if(ix>=0x4b800000) {
|
||||
y = zero; n = 0; /* y must be even */
|
||||
} else {
|
||||
if(ix<0x4b000000) z = y+two23; /* exact */
|
||||
GET_FLOAT_WORD(n,z);
|
||||
n &= 1;
|
||||
y = n;
|
||||
n<<= 2;
|
||||
}
|
||||
}
|
||||
switch (n) {
|
||||
case 0: y = __kernel_sindf(pi*y); break;
|
||||
case 1:
|
||||
@ -157,7 +155,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
|
||||
if(hx<0) {
|
||||
if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
|
||||
return one/zero;
|
||||
t = sin_pif(x);
|
||||
t = sin_pi(x);
|
||||
if(t==zero) return one/zero; /* -integer */
|
||||
nadj = __ieee754_logf(pi/fabsf(t*x));
|
||||
if(t<zero) *signgamp = -1;
|
||||
|
Loading…
Reference in New Issue
Block a user