|
 |
ce426f |
# commit 62a728aeff93507ce5975f245a5f1d2046fb4503
|
|
|
ce426f |
# Author: Alan Modra <amodra@gmail.com>
|
|
|
ce426f |
# Date: Sat Aug 17 18:27:19 2013 +0930
|
|
|
ce426f |
#
|
|
|
ce426f |
# PowerPC floating point little-endian [6 of 15]
|
|
|
ce426f |
# http://sourceware.org/ml/libc-alpha/2013-07/msg00197.html
|
|
|
ce426f |
#
|
|
|
ce426f |
# A rewrite to make this code correct for little-endian.
|
|
|
ce426f |
#
|
|
|
ce426f |
# * sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c (mynumber): Replace
|
|
|
ce426f |
# union 32-bit int array member with 64-bit int array.
|
|
|
ce426f |
# (t515, tm256): Double rather than long double.
|
|
|
ce426f |
# (__ieee754_sqrtl): Rewrite using 64-bit arithmetic.
|
|
|
ce426f |
#
|
|
|
ce426f |
diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c
|
|
|
ce426f |
--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c 2014-05-27 22:20:12.000000000 -0500
|
|
|
ce426f |
+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c 2014-05-27 22:21:39.000000000 -0500
|
|
|
ce426f |
@@ -34,15 +34,13 @@
|
|
|
ce426f |
|
|
|
ce426f |
#include <math_private.h>
|
|
|
ce426f |
|
|
|
ce426f |
-typedef unsigned int int4;
|
|
|
ce426f |
-typedef union {int4 i[4]; long double x; double d[2]; } mynumber;
|
|
|
ce426f |
+typedef union {int64_t i[2]; long double x; double d[2]; } mynumber;
|
|
|
ce426f |
|
|
|
ce426f |
-static const mynumber
|
|
|
ce426f |
- t512 = {{0x5ff00000, 0x00000000, 0x00000000, 0x00000000 }}, /* 2^512 */
|
|
|
ce426f |
- tm256 = {{0x2ff00000, 0x00000000, 0x00000000, 0x00000000 }}; /* 2^-256 */
|
|
|
ce426f |
static const double
|
|
|
ce426f |
-two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */
|
|
|
ce426f |
-twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */
|
|
|
ce426f |
+ t512 = 0x1p512,
|
|
|
ce426f |
+ tm256 = 0x1p-256,
|
|
|
ce426f |
+ two54 = 0x1p54, /* 0x4350000000000000 */
|
|
|
ce426f |
+ twom54 = 0x1p-54; /* 0x3C90000000000000 */
|
|
|
ce426f |
|
|
|
ce426f |
/*********************************************************************/
|
|
|
ce426f |
/* An ultimate sqrt routine. Given an IEEE double machine number x */
|
|
|
ce426f |
@@ -54,56 +52,53 @@
|
|
|
ce426f |
static const long double big = 134217728.0, big1 = 134217729.0;
|
|
|
ce426f |
long double t,s,i;
|
|
|
ce426f |
mynumber a,c;
|
|
|
ce426f |
- int4 k, l, m;
|
|
|
ce426f |
- int n;
|
|
|
ce426f |
+ uint64_t k, l;
|
|
|
ce426f |
+ int64_t m, n;
|
|
|
ce426f |
double d;
|
|
|
ce426f |
|
|
|
ce426f |
a.x=x;
|
|
|
ce426f |
- k=a.i[0] & 0x7fffffff;
|
|
|
ce426f |
+ k=a.i[0] & INT64_C(0x7fffffffffffffff);
|
|
|
ce426f |
/*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/
|
|
|
ce426f |
- if (k>0x000fffff && k<0x7ff00000) {
|
|
|
ce426f |
+ if (k>INT64_C(0x000fffff00000000) && k
|
|
|
ce426f |
if (x < 0) return (big1-big1)/(big-big);
|
|
|
ce426f |
- l = (k&0x001fffff)|0x3fe00000;
|
|
|
ce426f |
- if (((a.i[2] & 0x7fffffff) | a.i[3]) != 0) {
|
|
|
ce426f |
- n = (int) ((l - k) * 2) >> 21;
|
|
|
ce426f |
- m = (a.i[2] >> 20) & 0x7ff;
|
|
|
ce426f |
+ l = (k&INT64_C(0x001fffffffffffff))|INT64_C(0x3fe0000000000000);
|
|
|
ce426f |
+ if ((a.i[1] & INT64_C(0x7fffffffffffffff)) != 0) {
|
|
|
ce426f |
+ n = (int64_t) ((l - k) * 2) >> 53;
|
|
|
ce426f |
+ m = (a.i[1] >> 52) & 0x7ff;
|
|
|
ce426f |
if (m == 0) {
|
|
|
ce426f |
a.d[1] *= two54;
|
|
|
ce426f |
- m = ((a.i[2] >> 20) & 0x7ff) - 54;
|
|
|
ce426f |
+ m = ((a.i[1] >> 52) & 0x7ff) - 54;
|
|
|
ce426f |
}
|
|
|
ce426f |
m += n;
|
|
|
ce426f |
- if ((int) m > 0)
|
|
|
ce426f |
- a.i[2] = (a.i[2] & 0x800fffff) | (m << 20);
|
|
|
ce426f |
- else if ((int) m <= -54) {
|
|
|
ce426f |
- a.i[2] &= 0x80000000;
|
|
|
ce426f |
- a.i[3] = 0;
|
|
|
ce426f |
+ if (m > 0)
|
|
|
ce426f |
+ a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
|
|
|
ce426f |
+ else if (m <= -54) {
|
|
|
ce426f |
+ a.i[1] &= INT64_C(0x8000000000000000);
|
|
|
ce426f |
} else {
|
|
|
ce426f |
m += 54;
|
|
|
ce426f |
- a.i[2] = (a.i[2] & 0x800fffff) | (m << 20);
|
|
|
ce426f |
+ a.i[1] = (a.i[1] & INT64_C(0x800fffffffffffff)) | (m << 52);
|
|
|
ce426f |
a.d[1] *= twom54;
|
|
|
ce426f |
}
|
|
|
ce426f |
}
|
|
|
ce426f |
a.i[0] = l;
|
|
|
ce426f |
s = a.x;
|
|
|
ce426f |
d = __ieee754_sqrt (a.d[0]);
|
|
|
ce426f |
- c.i[0] = 0x20000000+((k&0x7fe00000)>>1);
|
|
|
ce426f |
+ c.i[0] = INT64_C(0x2000000000000000)+((k&INT64_C(0x7fe0000000000000))>>1);
|
|
|
ce426f |
c.i[1] = 0;
|
|
|
ce426f |
- c.i[2] = 0;
|
|
|
ce426f |
- c.i[3] = 0;
|
|
|
ce426f |
i = d;
|
|
|
ce426f |
t = 0.5L * (i + s / i);
|
|
|
ce426f |
i = 0.5L * (t + s / t);
|
|
|
ce426f |
return c.x * i;
|
|
|
ce426f |
}
|
|
|
ce426f |
else {
|
|
|
ce426f |
- if (k>=0x7ff00000) {
|
|
|
ce426f |
- if (a.i[0] == 0xfff00000 && a.i[1] == 0)
|
|
|
ce426f |
+ if (k>=INT64_C(0x7ff0000000000000)) {
|
|
|
ce426f |
+ if (a.i[0] == INT64_C(0xfff0000000000000))
|
|
|
ce426f |
return (big1-big1)/(big-big); /* sqrt (-Inf) = NaN. */
|
|
|
ce426f |
return x; /* sqrt (NaN) = NaN, sqrt (+Inf) = +Inf. */
|
|
|
ce426f |
}
|
|
|
ce426f |
if (x == 0) return x;
|
|
|
ce426f |
if (x < 0) return (big1-big1)/(big-big);
|
|
|
ce426f |
- return tm256.x*__ieee754_sqrtl(x*t512.x);
|
|
|
ce426f |
+ return tm256*__ieee754_sqrtl(x*t512);
|
|
|
ce426f |
}
|
|
|
ce426f |
}
|
|
|
ce426f |
strong_alias (__ieee754_sqrtl, __sqrtl_finite)
|