X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;ds=sidebyside;f=src%2Fmath%2F__rem_pio2.c;h=a40db9fc05fdd028d90f11d8f9ce8f8db6992fcf;hb=b72cd07f176b876aa51864d93aa8101477b1d732;hp=a7d779e0292ed48a2a12a0e49f3930bfd1b35503;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd;p=musl

diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index a7d779e0..a40db9fc 100644
--- a/src/math/__rem_pio2.c
+++ b/src/math/__rem_pio2.c
@@ -19,6 +19,12 @@
 
 #include "libm.h"
 
+#if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
+#define EPS DBL_EPSILON
+#elif FLT_EVAL_METHOD==2
+#define EPS LDBL_EPSILON
+#endif
+
 /*
  * invpio2:  53 bits of 2/pi
  * pio2_1:   first  33 bit of pi/2
@@ -29,8 +35,7 @@
  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  */
 static const double
-zero    = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-two24   = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+toint   = 1.5/EPS,
 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
 pio2_1  = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
@@ -42,18 +47,19 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 /* caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
 int __rem_pio2(double x, double *y)
 {
-	double z,w,t,r,fn;
+	union {double f; uint64_t i;} u = {x};
+	double_t z,w,t,r,fn;
 	double tx[3],ty[2];
-	int32_t e0,i,j,nx,n,ix,hx;
-	uint32_t low;
+	uint32_t ix;
+	int sign, n, ex, ey, i;
 
-	GET_HIGH_WORD(hx,x);
-	ix = hx & 0x7fffffff;
+	sign = u.i>>63;
+	ix = u.i>>32 & 0x7fffffff;
 	if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
 		if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
 			goto medium;  /* cancellation -- use medium case */
 		if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
-			if (hx > 0) {
+			if (!sign) {
 				z = x - pio2_1;  /* one round good to 85 bits */
 				y[0] = z - pio2_1t;
 				y[1] = (z-y[0]) - pio2_1t;
@@ -65,7 +71,7 @@ int __rem_pio2(double x, double *y)
 				return -1;
 			}
 		} else {
-			if (hx > 0) {
+			if (!sign) {
 				z = x - 2*pio2_1;
 				y[0] = z - 2*pio2_1t;
 				y[1] = (z-y[0]) - 2*pio2_1t;
@@ -82,7 +88,7 @@ int __rem_pio2(double x, double *y)
 		if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
 			if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
 				goto medium;
-			if (hx > 0) {
+			if (!sign) {
 				z = x - 3*pio2_1;
 				y[0] = z - 3*pio2_1t;
 				y[1] = (z-y[0]) - 3*pio2_1t;
@@ -96,7 +102,7 @@ int __rem_pio2(double x, double *y)
 		} else {
 			if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
 				goto medium;
-			if (hx > 0) {
+			if (!sign) {
 				z = x - 4*pio2_1;
 				y[0] = z - 4*pio2_1t;
 				y[1] = (z-y[0]) - 4*pio2_1t;
@@ -110,32 +116,25 @@ int __rem_pio2(double x, double *y)
 		}
 	}
 	if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
-		uint32_t high;
 medium:
-		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-		STRICT_ASSIGN(double, fn, x*invpio2 + 0x1.8p52);
-		fn = fn - 0x1.8p52;
-// FIXME
-#ifdef HAVE_EFFICIENT_IRINT
-		n = irint(fn);
-#else
+		/* rint(x/(pi/2)), Assume round-to-nearest. */
+		fn = x*invpio2 + toint - toint;
 		n = (int32_t)fn;
-#endif
 		r = x - fn*pio2_1;
 		w = fn*pio2_1t;  /* 1st round, good to 85 bits */
-		j = ix>>20;
 		y[0] = r - w;
-		GET_HIGH_WORD(high,y[0]);
-		i = j - ((high>>20)&0x7ff);
-		if (i > 16) {  /* 2nd round, good to 118 bits */
+		u.f = y[0];
+		ey = u.i>>52 & 0x7ff;
+		ex = ix>>20;
+		if (ex - ey > 16) { /* 2nd round, good to 118 bits */
 			t = r;
 			w = fn*pio2_2;
 			r = t - w;
 			w = fn*pio2_2t - ((t-r)-w);
 			y[0] = r - w;
-			GET_HIGH_WORD(high,y[0]);
-			i = j - ((high>>20)&0x7ff);
-			if (i > 49) {  /* 3rd round, good to 151 bits, covers all cases */
+			u.f = y[0];
+			ey = u.i>>52 & 0x7ff;
+			if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
 				t = r;
 				w = fn*pio2_3;
 				r = t - w;
@@ -143,7 +142,7 @@ medium:
 				y[0] = r - w;
 			}
 		}
-		y[1] = (r-y[0]) - w;
+		y[1] = (r - y[0]) - w;
 		return n;
 	}
 	/*
@@ -153,19 +152,21 @@ medium:
 		y[0] = y[1] = x - x;
 		return 0;
 	}
-	/* set z = scalbn(|x|,ilogb(x)-23) */
-	GET_LOW_WORD(low,x);
-	e0 = (ix>>20) - 1046;  /* e0 = ilogb(z)-23; */
-	INSERT_WORDS(z, ix - ((int32_t)(e0<<20)), low);
-	for (i=0; i<2; i++) {
-		tx[i] = (double)((int32_t)(z));
-		z = (z-tx[i])*two24;
+	/* set z = scalbn(|x|,-ilogb(x)+23) */
+	u.f = x;
+	u.i &= (uint64_t)-1>>12;
+	u.i |= (uint64_t)(0x3ff + 23)<<52;
+	z = u.f;
+	for (i=0; i < 2; i++) {
+		tx[i] = (double)(int32_t)z;
+		z     = (z-tx[i])*0x1p24;
 	}
-	tx[2] = z;
-	nx = 3;
-	while (tx[nx-1] == zero) nx--;  /* skip zero term */
-	n = __rem_pio2_large(tx,ty,e0,nx,1);
-	if (hx < 0) {
+	tx[i] = z;
+	/* skip zero terms, first term is non-zero */
+	while (tx[i] == 0.0)
+		i--;
+	n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
+	if (sign) {
 		y[0] = -ty[0];
 		y[1] = -ty[1];
 		return -n;