diff options
| -rw-r--r-- | ChangeLog | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c | 51 |
2 files changed, 30 insertions, 28 deletions
@@ -1,5 +1,12 @@ 2013-10-04 Alan Modra <amodra@gmail.com> + * sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c (mynumber): Replace + union 32-bit int array member with 64-bit int array. + (t515, tm256): Double rather than long double. + (__ieee754_sqrtl): Rewrite using 64-bit arithmetic. + +2013-10-04 Alan Modra <amodra@gmail.com> + * sysdeps/ieee754/ldbl-128ibm/ieee754.h (union ieee854_long_double): Delete. (IEEE854_LONG_DOUBLE_BIAS): Delete. diff --git a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c index 2b0f7c6..61feb36 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c @@ -34,15 +34,13 @@ #include <math_private.h> -typedef unsigned int int4; -typedef union {int4 i[4]; long double x; double d[2]; } mynumber; +typedef union {int64_t i[2]; long double x; double d[2]; } mynumber; -static const mynumber - t512 = {{0x5ff00000, 0x00000000, 0x00000000, 0x00000000 }}, /* 2^512 */ - tm256 = {{0x2ff00000, 0x00000000, 0x00000000, 0x00000000 }}; /* 2^-256 */ static const double -two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */ -twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */ + t512 = 0x1p512, + tm256 = 0x1p-256, + two54 = 0x1p54, /* 0x4350000000000000 */ + twom54 = 0x1p-54; /* 0x3C90000000000000 */ /*********************************************************************/ /* An ultimate sqrt routine. Given an IEEE double machine number x */ @@ -54,56 +52,53 @@ long double __ieee754_sqrtl(long double x) static const long double big = 134217728.0, big1 = 134217729.0; long double t,s,i; mynumber a,c; - int4 k, l, m; - int n; + uint64_t k, l; + int64_t m, n; double d; a.x=x; - k=a.i[0] & 0x7fffffff; + k=a.i[0] & INT64_C(0x7fffffffffffffff); /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ - if (k>0x000fffff && k<0x7ff00000) { + if (k>INT64_C(0x000fffff00000000) && k<INT64_C(0x7ff0000000000000)) { if (x < 0) return (big1-big1)/(big-big); - l = (k&0x001fffff)|0x3fe00000; - if (((a.i[2] & 0x7fffffff) | a.i[3]) != 0) { - n = (int) ((l - k) * 2) >> 21; - m = (a.i[2] >> 20) & 0x7ff; + l = (k&INT64_C(0x001fffffffffffff))|INT64_C(0x3fe0000000000000); + if ((a.i[1] & INT64_C(0x7fffffffffffffff)) != 0) { + n = (int64_t) ((l - k) * 2) >> 53; + m = (a.i[1] >> 52) & 0x7ff; if (m == 0) { a.d[1] *= two54; - m = ((a.i[2] >> 20) & 0x7ff) - 54; + m = ((a.i[1] >> 52) & 0x7ff) - 54; } m += n; - if ((int) m > 0) - a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); - else if ((int) m <= -54) { - a.i[2] &= 0x80000000; - a.i[3] = 0; + if (m > 0) + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); + else if (m <= -54) { + a.i[1] &= INT64_C(0x8000000000000000); } else { m += 54; - a.i[2] = (a.i[2] & 0x800fffff) | (m << 20); + a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52); a.d[1] *= twom54; } } a.i[0] = l; s = a.x; d = __ieee754_sqrt (a.d[0]); - c.i[0] = 0x20000000+((k&0x7fe00000)>>1); + c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1); c.i[1] = 0; - c.i[2] = 0; - c.i[3] = 0; i = d; t = 0.5L * (i + s / i); i = 0.5L * (t + s / t); return c.x * i; } else { - if (k>=0x7ff00000) { - if (a.i[0] == 0xfff00000 && a.i[1] == 0) + if (k>=INT64_C(0x7ff0000000000000)) { + if (a.i[0] == INT64_C(0xfff0000000000000)) return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN. */ return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */ } if (x == 0) return x; if (x < 0) return (big1-big1)/(big-big); - return tm256.x*__ieee754_sqrtl(x*t512.x); + return tm256*__ieee754_sqrtl(x*t512); } } strong_alias (__ieee754_sqrtl, __sqrtl_finite) |