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
3.9 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)