OpenE2K
/
glibc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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.
This commit is contained in:
Alan Modra 2013-08-17 18:27:19 +09:30
parent 32c301dfc9
commit 62a728aeff
2 changed files with 30 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
ChangeLog
View File

@ -1,3 +1,10 @@
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):

51
sysdeps/ieee754/ldbl-128ibm/e_sqrtl.c
View File

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