You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
113 lines
4.0 KiB
113 lines
4.0 KiB
# commit 62a728aeff93507ce5975f245a5f1d2046fb4503 |
|
# Author: Alan Modra <amodra@gmail.com> |
|
# Date: Sat Aug 17 18:27:19 2013 +0930 |
|
# |
|
# PowerPC floating point little-endian [6 of 15] |
|
# http://sourceware.org/ml/libc-alpha/2013-07/msg00197.html |
|
# |
|
# A rewrite to make this code correct for little-endian. |
|
# |
|
# * 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. |
|
# |
|
diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c |
|
--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c 2014-05-27 22:20:12.000000000 -0500 |
|
+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c 2014-05-27 22:21:39.000000000 -0500 |
|
@@ -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 @@ |
|
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)
|
|
|