Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,29 ****
--- 20,30 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
*** 32,42 ****
#include "libm.h"
#include "fma.h"
#include "fenv_inlines.h"
#if defined(__sparc)
-
static const union {
unsigned i[2];
double d;
} C[] = {
{ 0x3fe00000u, 0 },
--- 33,42 ----
*** 78,96 ****
/*
* fmal for SPARC: 128-bit quad precision, big-endian
*/
long double
! __fmal(long double x, long double y, long double z) {
union {
unsigned int i[4];
long double q;
} xx, yy, zz;
union {
unsigned int i[2];
double d;
} u;
double dx[5], dy[5], dxy[9], c, s;
unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
unsigned int rm, sticky;
unsigned int fsr;
--- 78,98 ----
/*
* fmal for SPARC: 128-bit quad precision, big-endian
*/
long double
! __fmal(long double x, long double y, long double z)
! {
union {
unsigned int i[4];
long double q;
} xx, yy, zz;
union {
unsigned int i[2];
double d;
} u;
+
double dx[5], dy[5], dxy[9], c, s;
unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
unsigned int rm, sticky;
unsigned int fsr;
*** 117,166 ****
dummy = snan;
dummy += snan;
xx.i[0] |= 0x8000;
return (xx.q);
}
cx = 3; /* quiet nan */
! } else
cx = 2; /* inf */
} else if (hx == 0) {
cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
/* subnormal or zero */
! } else
cx = 1; /* finite nonzero */
if (hy >= 0x7fff0000) {
if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
if (!(hy & 0x8000)) {
dummy = snan;
dummy += snan;
yy.i[0] |= 0x8000;
return (yy.q);
}
cy = 3;
! } else
cy = 2;
} else if (hy == 0) {
cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
! } else
cy = 1;
if (hz >= 0x7fff0000) {
if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
if (!(hz & 0x8000)) {
dummy = snan;
dummy += snan;
zz.i[0] |= 0x8000;
return (zz.q);
}
cz = 3;
! } else
cz = 2;
} else if (hz == 0) {
cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
! } else
cz = 1;
/* get the fsr and clear current exceptions */
__fenv_getfsr32(&fsr);
fsr &= ~FSR_CEXC;
--- 119,177 ----
dummy = snan;
dummy += snan;
xx.i[0] |= 0x8000;
return (xx.q);
}
+
cx = 3; /* quiet nan */
! } else {
cx = 2; /* inf */
+ }
} else if (hx == 0) {
cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
/* subnormal or zero */
! } else {
cx = 1; /* finite nonzero */
+ }
if (hy >= 0x7fff0000) {
if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
if (!(hy & 0x8000)) {
dummy = snan;
dummy += snan;
yy.i[0] |= 0x8000;
return (yy.q);
}
+
cy = 3;
! } else {
cy = 2;
+ }
} else if (hy == 0) {
cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
! } else {
cy = 1;
+ }
if (hz >= 0x7fff0000) {
if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
if (!(hz & 0x8000)) {
dummy = snan;
dummy += snan;
zz.i[0] |= 0x8000;
return (zz.q);
}
+
cz = 3;
! } else {
cz = 2;
+ }
} else if (hz == 0) {
cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
! } else {
cz = 1;
+ }
/* get the fsr and clear current exceptions */
__fenv_getfsr32(&fsr);
fsr &= ~FSR_CEXC;
*** 169,178 ****
--- 180,190 ----
/* if x or y is a quiet nan, return it */
if (cx == 3) {
__fenv_setfsr32(&fsr);
return (x);
}
+
if (cy == 3) {
__fenv_setfsr32(&fsr);
return (y);
}
*** 199,216 ****
/*
* if z is also inf, either we have inf-inf or
* the result is the same as z depending on signs
*/
if (cz == 2) {
! if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
dummy = inf;
dummy -= inf;
zz.i[0] = 0x7fffffff;
zz.i[1] = zz.i[2] = zz.i[3] =
0xffffffff;
return (zz.q);
}
__fenv_setfsr32(&fsr);
return (z);
}
/* otherwise the result is inf with appropriate sign */
--- 211,229 ----
/*
* if z is also inf, either we have inf-inf or
* the result is the same as z depending on signs
*/
if (cz == 2) {
! if ((int)((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
dummy = inf;
dummy -= inf;
zz.i[0] = 0x7fffffff;
zz.i[1] = zz.i[2] = zz.i[3] =
0xffffffff;
return (zz.q);
}
+
__fenv_setfsr32(&fsr);
return (z);
}
/* otherwise the result is inf with appropriate sign */
*** 231,247 ****
* now x, y, and z are all finite; handle cases where x or y
* is zero
*/
if (cx == 0 || cy == 0) {
/* either we have 0-0 or the result is the same as z */
! if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
0) {
zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
0;
__fenv_setfsr32(&fsr);
return (zz.q);
}
__fenv_setfsr32(&fsr);
return (z);
}
/* if we get here, x and y are nonzero finite, z must be zero */
--- 244,261 ----
* now x, y, and z are all finite; handle cases where x or y
* is zero
*/
if (cx == 0 || cy == 0) {
/* either we have 0-0 or the result is the same as z */
! if (cz == 0 && (int)((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
0) {
zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
0;
__fenv_setfsr32(&fsr);
return (zz.q);
}
+
__fenv_setfsr32(&fsr);
return (z);
}
/* if we get here, x and y are nonzero finite, z must be zero */
*** 259,268 ****
--- 273,283 ----
* of x and y
*/
sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000;
ex = hx >> 16;
hx &= 0xffff;
+
if (!ex) {
if (hx | (xx.i[1] & 0xfffe0000)) {
ex = 1;
} else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) {
hx = xx.i[1];
*** 278,298 ****
} else {
hx = xx.i[3];
xx.i[1] = xx.i[2] = xx.i[3] = 0;
ex = -95;
}
while ((hx & 0x10000) == 0) {
hx = (hx << 1) | (xx.i[1] >> 31);
xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
xx.i[3] <<= 1;
ex--;
}
! } else
hx |= 0x10000;
ey = hy >> 16;
hy &= 0xffff;
if (!ey) {
if (hy | (yy.i[1] & 0xfffe0000)) {
ey = 1;
} else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
hy = yy.i[1];
--- 293,317 ----
} else {
hx = xx.i[3];
xx.i[1] = xx.i[2] = xx.i[3] = 0;
ex = -95;
}
+
while ((hx & 0x10000) == 0) {
hx = (hx << 1) | (xx.i[1] >> 31);
xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
xx.i[3] <<= 1;
ex--;
}
! } else {
hx |= 0x10000;
+ }
+
ey = hy >> 16;
hy &= 0xffff;
+
if (!ey) {
if (hy | (yy.i[1] & 0xfffe0000)) {
ey = 1;
} else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
hy = yy.i[1];
*** 308,363 ****
} else {
hy = yy.i[3];
yy.i[1] = yy.i[2] = yy.i[3] = 0;
ey = -95;
}
while ((hy & 0x10000) == 0) {
hy = (hy << 1) | (yy.i[1] >> 31);
yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
yy.i[3] <<= 1;
ey--;
}
! } else
hy |= 0x10000;
exy = ex + ey - 0x3fff;
/* convert the significands of x and y to doubles */
c = twom16;
! dx[0] = (double) ((int) hx) * c;
! dy[0] = (double) ((int) hy) * c;
c *= twom24;
! dx[1] = (double) ((int) (xx.i[1] >> 8)) * c;
! dy[1] = (double) ((int) (yy.i[1] >> 8)) * c;
c *= twom24;
! dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) &
0xffffff)) * c;
! dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) &
0xffffff)) * c;
c *= twom24;
! dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) &
! 0xffffff)) * c;
! dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) &
! 0xffffff)) * c;
c *= twom24;
! dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c;
! dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c;
/* form the "digits" of the product */
dxy[0] = dx[0] * dy[0];
dxy[1] = dx[0] * dy[1] + dx[1] * dy[0];
dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
! dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] +
! dx[3] * dy[0];
! dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] +
! dx[3] * dy[1] + dx[4] * dy[0];
! dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] +
! dx[4] * dy[1];
dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
dxy[7] = dx[3] * dy[4] + dx[4] * dy[3];
dxy[8] = dx[4] * dy[4];
/* split odd-numbered terms and combine into even-numbered terms */
--- 327,383 ----
} else {
hy = yy.i[3];
yy.i[1] = yy.i[2] = yy.i[3] = 0;
ey = -95;
}
+
while ((hy & 0x10000) == 0) {
hy = (hy << 1) | (yy.i[1] >> 31);
yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
yy.i[3] <<= 1;
ey--;
}
! } else {
hy |= 0x10000;
+ }
+
exy = ex + ey - 0x3fff;
/* convert the significands of x and y to doubles */
c = twom16;
! dx[0] = (double)((int)hx) * c;
! dy[0] = (double)((int)hy) * c;
c *= twom24;
! dx[1] = (double)((int)(xx.i[1] >> 8)) * c;
! dy[1] = (double)((int)(yy.i[1] >> 8)) * c;
c *= twom24;
! dx[2] = (double)((int)(((xx.i[1] << 16) | (xx.i[2] >> 16)) &
0xffffff)) * c;
! dy[2] = (double)((int)(((yy.i[1] << 16) | (yy.i[2] >> 16)) &
0xffffff)) * c;
c *= twom24;
! dx[3] = (double)((int)(((xx.i[2] << 8) | (xx.i[3] >> 24)) & 0xffffff)) *
! c;
! dy[3] = (double)((int)(((yy.i[2] << 8) | (yy.i[3] >> 24)) & 0xffffff)) *
! c;
c *= twom24;
! dx[4] = (double)((int)(xx.i[3] & 0xffffff)) * c;
! dy[4] = (double)((int)(yy.i[3] & 0xffffff)) * c;
/* form the "digits" of the product */
dxy[0] = dx[0] * dy[0];
dxy[1] = dx[0] * dy[1] + dx[1] * dy[0];
dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
! dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] + dx[3] * dy[0];
! dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] + dx[3] * dy[1] +
! dx[4] * dy[0];
! dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] + dx[4] * dy[1];
dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
dxy[7] = dx[3] * dy[4] + dx[4] * dy[3];
dxy[8] = dx[4] * dy[4];
/* split odd-numbered terms and combine into even-numbered terms */
*** 377,386 ****
--- 397,407 ----
/* propagate carries, adjusting the exponent if need be */
dxy[7] = dxy[6] + dxy[8];
dxy[5] = dxy[4] + dxy[7];
dxy[3] = dxy[2] + dxy[5];
dxy[1] = dxy[0] + dxy[3];
+
if (dxy[1] >= two) {
dxy[0] *= half;
dxy[1] *= half;
dxy[2] *= half;
dxy[3] *= half;
*** 446,455 ****
--- 467,477 ----
/* extract the sign, exponent, and significand of z */
sz = zz.i[0] & 0x80000000;
ez = hz >> 16;
z0 = hz & 0xffff;
+
if (!ez) {
if (z0 | (zz.i[1] & 0xfffe0000)) {
z1 = zz.i[1];
z2 = zz.i[2];
z3 = zz.i[3];
*** 468,477 ****
--- 490,500 ----
} else {
z0 = zz.i[3];
z1 = z2 = z3 = 0;
ez = -95;
}
+
while ((z0 & 0x10000) == 0) {
z0 = (z0 << 1) | (z1 >> 31);
z1 = (z1 << 1) | (z2 >> 31);
z2 = (z2 << 1) | (z3 >> 31);
z3 <<= 1;
*** 481,587 ****
z0 |= 0x10000;
z1 = zz.i[1];
z2 = zz.i[2];
z3 = zz.i[3];
}
z4 = z5 = z6 = z7 = 0;
/*
* now x*y is represented by sxy, exy, and xy[0-7], and z is
* represented likewise; swap if need be so |xy| <= |z|
*/
if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
! (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 ||
! (xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
! e = sxy; sxy = sz; sz = e;
! e = exy; exy = ez; ez = e;
! e = xy0; xy0 = z0; z0 = e;
! e = xy1; xy1 = z1; z1 = e;
! e = xy2; xy2 = z2; z2 = e;
! e = xy3; xy3 = z3; z3 = e;
! z4 = xy4; xy4 = 0;
! z5 = xy5; xy5 = 0;
! z6 = xy6; xy6 = 0;
! z7 = xy7; xy7 = 0;
}
/* shift the significand of xy keeping a sticky bit */
e = ez - exy;
if (e > 236) {
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
xy7 = 1;
} else if (e >= 224) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 |
! ((xy0 << 1) << (255 - e));
xy7 = xy0 >> (e - 224);
if (sticky)
xy7 |= 1;
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
} else if (e >= 192) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 |
! ((xy1 << 1) << (223 - e));
xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
if (sticky)
xy7 |= 1;
xy6 = xy0 >> (e - 192);
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0;
} else if (e >= 160) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 |
! ((xy2 << 1) << (191 - e));
xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
xy5 = xy0 >> (e - 160);
xy0 = xy1 = xy2 = xy3 = xy4 = 0;
} else if (e >= 128) {
sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
xy4 = xy0 >> (e - 128);
xy0 = xy1 = xy2 = xy3 = 0;
} else if (e >= 96) {
sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
xy3 = xy0 >> (e - 96);
xy0 = xy1 = xy2 = 0;
} else if (e >= 64) {
sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
xy2 = xy0 >> (e - 64);
xy0 = xy1 = 0;
} else if (e >= 32) {
sticky = xy7 | ((xy6 << 1) << (63 - e));
xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
xy1 = xy0 >> (e - 32);
xy0 = 0;
} else if (e) {
sticky = (xy7 << 1) << (31 - e);
xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
if (sticky)
xy7 |= 1;
xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
--- 504,644 ----
z0 |= 0x10000;
z1 = zz.i[1];
z2 = zz.i[2];
z3 = zz.i[3];
}
+
z4 = z5 = z6 = z7 = 0;
/*
* now x*y is represented by sxy, exy, and xy[0-7], and z is
* represented likewise; swap if need be so |xy| <= |z|
*/
if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
! (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 || (xy3 == z3 &&
! (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
! e = sxy;
! sxy = sz;
! sz = e;
! e = exy;
! exy = ez;
! ez = e;
! e = xy0;
! xy0 = z0;
! z0 = e;
! e = xy1;
! xy1 = z1;
! z1 = e;
! e = xy2;
! xy2 = z2;
! z2 = e;
! e = xy3;
! xy3 = z3;
! z3 = e;
! z4 = xy4;
! xy4 = 0;
! z5 = xy5;
! xy5 = 0;
! z6 = xy6;
! xy6 = 0;
! z7 = xy7;
! xy7 = 0;
}
/* shift the significand of xy keeping a sticky bit */
e = ez - exy;
+
if (e > 236) {
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
xy7 = 1;
} else if (e >= 224) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 | ((xy0 <<
! 1) << (255 - e));
xy7 = xy0 >> (e - 224);
+
if (sticky)
xy7 |= 1;
+
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
} else if (e >= 192) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | ((xy1 << 1) <<
! (223 - e));
xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = xy0 >> (e - 192);
xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0;
} else if (e >= 160) {
! sticky = xy7 | xy6 | xy5 | xy4 | xy3 | ((xy2 << 1) << (191 -
! e));
xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
xy5 = xy0 >> (e - 160);
xy0 = xy1 = xy2 = xy3 = xy4 = 0;
} else if (e >= 128) {
sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
xy4 = xy0 >> (e - 128);
xy0 = xy1 = xy2 = xy3 = 0;
} else if (e >= 96) {
sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
xy3 = xy0 >> (e - 96);
xy0 = xy1 = xy2 = 0;
} else if (e >= 64) {
sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
xy2 = xy0 >> (e - 64);
xy0 = xy1 = 0;
} else if (e >= 32) {
sticky = xy7 | ((xy6 << 1) << (63 - e));
xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
xy1 = xy0 >> (e - 32);
xy0 = 0;
} else if (e) {
sticky = (xy7 << 1) << (31 - e);
xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
+
if (sticky)
xy7 |= 1;
+
xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
*** 597,606 ****
--- 654,664 ----
xy3 = ~xy3;
xy4 = ~xy4;
xy5 = ~xy5;
xy6 = ~xy6;
xy7 = -xy7;
+
if (xy7 == 0)
if (++xy6 == 0)
if (++xy5 == 0)
if (++xy4 == 0)
if (++xy3 == 0)
*** 611,717 ****
/* add, propagating carries */
z7 += xy7;
e = (z7 < xy7);
z6 += xy6;
if (e) {
z6++;
e = (z6 <= xy6);
! } else
e = (z6 < xy6);
z5 += xy5;
if (e) {
z5++;
e = (z5 <= xy5);
! } else
e = (z5 < xy5);
z4 += xy4;
if (e) {
z4++;
e = (z4 <= xy4);
! } else
e = (z4 < xy4);
z3 += xy3;
if (e) {
z3++;
e = (z3 <= xy3);
! } else
e = (z3 < xy3);
z2 += xy2;
if (e) {
z2++;
e = (z2 <= xy2);
! } else
e = (z2 < xy2);
z1 += xy1;
if (e) {
z1++;
e = (z1 <= xy1);
! } else
e = (z1 < xy1);
z0 += xy0;
if (e)
z0++;
/* postnormalize and collect rounding information into z4 */
if (ez < 1) {
/* result is tiny; shift right until exponent is within range */
e = 1 - ez;
if (e > 116) {
z4 = 1; /* result can't be exactly zero */
z0 = z1 = z2 = z3 = 0;
} else if (e >= 96) {
! sticky = z7 | z6 | z5 | z4 | z3 | z2 |
! ((z1 << 1) << (127 - e));
z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
if (sticky)
z4 |= 1;
z3 = z0 >> (e - 96);
z0 = z1 = z2 = 0;
} else if (e >= 64) {
! sticky = z7 | z6 | z5 | z4 | z3 |
! ((z2 << 1) << (95 - e));
z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
if (sticky)
z4 |= 1;
z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
z2 = z0 >> (e - 64);
z0 = z1 = 0;
} else if (e >= 32) {
sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
if (sticky)
z4 |= 1;
z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
z1 = z0 >> (e - 32);
z0 = 0;
} else {
sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
z4 = (z4 >> e) | ((z3 << 1) << (31 - e));
if (sticky)
z4 |= 1;
z3 = (z3 >> e) | ((z2 << 1) << (31 - e));
z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
z0 >>= e;
}
ez = 1;
} else if (z0 >= 0x20000) {
/* carry out; shift right by one */
sticky = (z4 & 1) | z5 | z6 | z7;
z4 = (z4 >> 1) | (z3 << 31);
if (sticky)
z4 |= 1;
z3 = (z3 >> 1) | (z2 << 31);
z2 = (z2 >> 1) | (z1 << 31);
z1 = (z1 >> 1) | (z0 << 31);
z0 >>= 1;
ez++;
} else {
! if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7)
! != 0) {
/*
* borrow/cancellation; shift left as much as
* exponent allows
*/
while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
--- 669,806 ----
/* add, propagating carries */
z7 += xy7;
e = (z7 < xy7);
z6 += xy6;
+
if (e) {
z6++;
e = (z6 <= xy6);
! } else {
e = (z6 < xy6);
+ }
+
z5 += xy5;
+
if (e) {
z5++;
e = (z5 <= xy5);
! } else {
e = (z5 < xy5);
+ }
+
z4 += xy4;
+
if (e) {
z4++;
e = (z4 <= xy4);
! } else {
e = (z4 < xy4);
+ }
+
z3 += xy3;
+
if (e) {
z3++;
e = (z3 <= xy3);
! } else {
e = (z3 < xy3);
+ }
+
z2 += xy2;
+
if (e) {
z2++;
e = (z2 <= xy2);
! } else {
e = (z2 < xy2);
+ }
+
z1 += xy1;
+
if (e) {
z1++;
e = (z1 <= xy1);
! } else {
e = (z1 < xy1);
+ }
+
z0 += xy0;
+
if (e)
z0++;
/* postnormalize and collect rounding information into z4 */
if (ez < 1) {
/* result is tiny; shift right until exponent is within range */
e = 1 - ez;
+
if (e > 116) {
z4 = 1; /* result can't be exactly zero */
z0 = z1 = z2 = z3 = 0;
} else if (e >= 96) {
! sticky = z7 | z6 | z5 | z4 | z3 | z2 | ((z1 << 1) <<
! (127 - e));
z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
+
if (sticky)
z4 |= 1;
+
z3 = z0 >> (e - 96);
z0 = z1 = z2 = 0;
} else if (e >= 64) {
! sticky = z7 | z6 | z5 | z4 | z3 | ((z2 << 1) << (95 -
! e));
z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
+
if (sticky)
z4 |= 1;
+
z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
z2 = z0 >> (e - 64);
z0 = z1 = 0;
} else if (e >= 32) {
sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
+
if (sticky)
z4 |= 1;
+
z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
z1 = z0 >> (e - 32);
z0 = 0;
} else {
sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
z4 = (z4 >> e) | ((z3 << 1) << (31 - e));
+
if (sticky)
z4 |= 1;
+
z3 = (z3 >> e) | ((z2 << 1) << (31 - e));
z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
z0 >>= e;
}
+
ez = 1;
} else if (z0 >= 0x20000) {
/* carry out; shift right by one */
sticky = (z4 & 1) | z5 | z6 | z7;
z4 = (z4 >> 1) | (z3 << 31);
+
if (sticky)
z4 |= 1;
+
z3 = (z3 >> 1) | (z2 << 31);
z2 = (z2 >> 1) | (z1 << 31);
z1 = (z1 >> 1) | (z0 << 31);
z0 >>= 1;
ez++;
} else {
! if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7) !=
! 0) {
/*
* borrow/cancellation; shift left as much as
* exponent allows
*/
while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
*** 723,732 ****
--- 812,822 ----
z5 = z6;
z6 = z7;
z7 = 0;
ez -= 32;
}
+
while (z0 < 0x10000 && ez > 1) {
z0 = (z0 << 1) | (z1 >> 31);
z1 = (z1 << 1) | (z2 >> 31);
z2 = (z2 << 1) | (z3 >> 31);
z3 = (z3 << 1) | (z4 >> 31);
*** 735,757 ****
z6 = (z6 << 1) | (z7 >> 31);
z7 <<= 1;
ez--;
}
}
if (z5 | z6 | z7)
z4 |= 1;
}
/* get the rounding mode */
rm = fsr >> 30;
/* strip off the integer bit, if there is one */
ibit = z0 & 0x10000;
! if (ibit)
z0 -= 0x10000;
! else {
ez = 0;
if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
zz.i[0] = rm == FSR_RM ? 0x80000000 : 0;
zz.i[1] = zz.i[2] = zz.i[3] = 0;
__fenv_setfsr32(&fsr);
return (zz.q);
--- 825,850 ----
z6 = (z6 << 1) | (z7 >> 31);
z7 <<= 1;
ez--;
}
}
+
if (z5 | z6 | z7)
z4 |= 1;
}
/* get the rounding mode */
rm = fsr >> 30;
/* strip off the integer bit, if there is one */
ibit = z0 & 0x10000;
!
! if (ibit) {
z0 -= 0x10000;
! } else {
ez = 0;
+
if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
zz.i[0] = rm == FSR_RM ? 0x80000000 : 0;
zz.i[1] = zz.i[2] = zz.i[3] = 0;
__fenv_setfsr32(&fsr);
return (zz.q);
*** 768,779 ****
/* round and raise exceptions */
if (z4) {
fsr |= FSR_NXC;
/* decide whether to round the fraction up */
! if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u ||
! (z4 == 0x80000000u && (z3 & 1))))) {
/* round up and renormalize if necessary */
if (++z3 == 0)
if (++z2 == 0)
if (++z1 == 0)
if (++z0 == 0x10000) {
--- 861,872 ----
/* round and raise exceptions */
if (z4) {
fsr |= FSR_NXC;
/* decide whether to round the fraction up */
! if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u || (z4 ==
! 0x80000000u && (z3 & 1))))) {
/* round up and renormalize if necessary */
if (++z3 == 0)
if (++z2 == 0)
if (++z1 == 0)
if (++z0 == 0x10000) {
*** 790,799 ****
--- 883,893 ----
zz.i[1] = zz.i[2] = zz.i[3] = 0;
} else {
zz.i[0] = sz | 0x7ffeffff;
zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
}
+
fsr |= FSR_OFC | FSR_NXC;
} else {
zz.i[0] = sz | (ez << 16) | z0;
zz.i[1] = z1;
zz.i[2] = z2;
*** 812,826 ****
--- 906,922 ----
}
/* restore the fsr and emulate exceptions as needed */
if ((fsr & FSR_CEXC) & (fsr >> 23)) {
__fenv_setfsr32(&fsr);
+
if (fsr & FSR_OFC) {
dummy = huge;
dummy *= huge;
} else if (fsr & FSR_UFC) {
dummy = tiny;
+
if (fsr & FSR_NXC)
dummy *= tiny;
else
dummy -= tiny2;
} else {
*** 829,843 ****
}
} else {
fsr |= (fsr & 0x1f) << 5;
__fenv_setfsr32(&fsr);
}
return (zz.q);
}
-
#elif defined(__x86)
-
static const union {
unsigned i[2];
double d;
} C[] = {
{ 0, 0x3fe00000u },
--- 925,938 ----
}
} else {
fsr |= (fsr & 0x1f) << 5;
__fenv_setfsr32(&fsr);
}
+
return (zz.q);
}
#elif defined(__x86)
static const union {
unsigned i[2];
double d;
} C[] = {
{ 0, 0x3fe00000u },
*** 869,883 ****
/*
* fmal for x86: 80-bit extended double precision, little-endian
*/
long double
! __fmal(long double x, long double y, long double z) {
union {
unsigned i[NI];
long double e;
} xx, yy, zz;
long double xhi, yhi, xlo, ylo, t;
unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
unsigned oldcwsw, cwsw, rm, sticky, carry;
int ex, ey, ez, exy, sxy, sz, e, tinyafter;
volatile double dummy;
--- 964,980 ----
/*
* fmal for x86: 80-bit extended double precision, little-endian
*/
long double
! __fmal(long double x, long double y, long double z)
! {
union {
unsigned i[NI];
long double e;
} xx, yy, zz;
+
long double xhi, yhi, xlo, ylo, t;
unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
unsigned oldcwsw, cwsw, rm, sticky, carry;
int ex, ey, ez, exy, sxy, sz, e, tinyafter;
volatile double dummy;
*** 910,927 ****
--- 1007,1027 ----
* normalize the significands (this will raise the denormal operand
* exception if need be)
*/
sxy = (xx.i[2] ^ yy.i[2]) & 0x8000;
sz = zz.i[2] & 0x8000;
+
if (!ex) {
xx.e = x * two63;
ex = (xx.i[2] & 0x7fff) - 63;
}
+
if (!ey) {
yy.e = y * two63;
ey = (yy.i[2] & 0x7fff) - 63;
}
+
if (!ez) {
zz.e = z * two63;
ez = (zz.i[2] & 0x7fff) - 63;
}
*** 943,952 ****
--- 1043,1053 ----
yhi = ((y + twom32) + two32) - two32;
xlo = x - xhi;
ylo = y - yhi;
x *= y;
y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo;
+
if (x >= two) {
x *= half;
y *= half;
exy++;
}
*** 966,1016 ****
/*
* now x*y is represented by sxy, exy, and xy[0-4], and z is
* represented likewise; swap if need be so |xy| <= |z|
*/
! if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 &&
! (xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) {
! e = sxy; sxy = sz; sz = e;
! e = exy; exy = ez; ez = e;
! e = xy0; xy0 = z0; z0 = e;
! e = xy1; xy1 = z1; z1 = e;
! z2 = xy2; xy2 = 0;
! z3 = xy3; xy3 = 0;
}
/* shift the significand of xy keeping a sticky bit */
e = ez - exy;
if (e > 130) {
xy0 = xy1 = xy2 = xy3 = 0;
xy4 = 1;
} else if (e >= 128) {
sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e));
xy4 = xy0 >> (e - 128);
if (sticky)
xy4 |= 1;
xy0 = xy1 = xy2 = xy3 = 0;
} else if (e >= 96) {
sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e));
xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
if (sticky)
xy4 |= 1;
xy3 = xy0 >> (e - 96);
xy0 = xy1 = xy2 = 0;
} else if (e >= 64) {
sticky = xy3 | ((xy2 << 1) << (95 - e));
xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
if (sticky)
xy4 |= 1;
xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
xy2 = xy0 >> (e - 64);
xy0 = xy1 = 0;
} else if (e >= 32) {
sticky = (xy3 << 1) << (63 - e);
xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
if (sticky)
xy4 |= 1;
xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
xy1 = xy0 >> (e - 32);
xy0 = 0;
} else if (e) {
--- 1067,1136 ----
/*
* now x*y is represented by sxy, exy, and xy[0-4], and z is
* represented likewise; swap if need be so |xy| <= |z|
*/
! if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
! (xy1 == z1 && (xy2 | xy3) != 0)))))) {
! e = sxy;
! sxy = sz;
! sz = e;
! e = exy;
! exy = ez;
! ez = e;
! e = xy0;
! xy0 = z0;
! z0 = e;
! e = xy1;
! xy1 = z1;
! z1 = e;
! z2 = xy2;
! xy2 = 0;
! z3 = xy3;
! xy3 = 0;
}
/* shift the significand of xy keeping a sticky bit */
e = ez - exy;
+
if (e > 130) {
xy0 = xy1 = xy2 = xy3 = 0;
xy4 = 1;
} else if (e >= 128) {
sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e));
xy4 = xy0 >> (e - 128);
+
if (sticky)
xy4 |= 1;
+
xy0 = xy1 = xy2 = xy3 = 0;
} else if (e >= 96) {
sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e));
xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
+
if (sticky)
xy4 |= 1;
+
xy3 = xy0 >> (e - 96);
xy0 = xy1 = xy2 = 0;
} else if (e >= 64) {
sticky = xy3 | ((xy2 << 1) << (95 - e));
xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
+
if (sticky)
xy4 |= 1;
+
xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
xy2 = xy0 >> (e - 64);
xy0 = xy1 = 0;
} else if (e >= 32) {
sticky = (xy3 << 1) << (63 - e);
xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
+
if (sticky)
xy4 |= 1;
+
xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
xy1 = xy0 >> (e - 32);
xy0 = 0;
} else if (e) {
*** 1026,1035 ****
--- 1146,1156 ----
xy0 = ~xy0;
xy1 = ~xy1;
xy2 = ~xy2;
xy3 = ~xy3;
xy4 = -xy4;
+
if (xy4 == 0)
if (++xy3 == 0)
if (++xy2 == 0)
if (++xy1 == 0)
xy0++;
*** 1037,1110 ****
/* add, propagating carries */
z4 += xy4;
carry = (z4 < xy4);
z3 += xy3;
if (carry) {
z3++;
carry = (z3 <= xy3);
! } else
carry = (z3 < xy3);
z2 += xy2;
if (carry) {
z2++;
carry = (z2 <= xy2);
! } else
carry = (z2 < xy2);
z1 += xy1;
if (carry) {
z1++;
carry = (z1 <= xy1);
! } else
carry = (z1 < xy1);
z0 += xy0;
if (carry) {
z0++;
carry = (z0 <= xy0);
! } else
carry = (z0 < xy0);
/* for a magnitude subtract, ignore the last carry out */
if (sxy ^ sz)
carry = 0;
/* postnormalize and collect rounding information into z2 */
if (ez < 1) {
/* result is tiny; shift right until exponent is within range */
e = 1 - ez;
if (e > 67) {
z2 = 1; /* result can't be exactly zero */
z0 = z1 = 0;
} else if (e >= 64) {
sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e));
z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e));
if (sticky)
z2 |= 1;
z1 = carry >> (e - 64);
z0 = 0;
} else if (e >= 32) {
sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e));
z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
if (sticky)
z2 |= 1;
z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e));
z0 = carry >> (e - 32);
} else {
sticky = z4 | z3 | (z2 << 1) << (31 - e);
z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
if (sticky)
z2 |= 1;
z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
z0 = (z0 >> e) | ((carry << 1) << (31 - e));
}
ez = 1;
} else if (carry) {
/* carry out; shift right by one */
sticky = (z2 & 1) | z3 | z4;
z2 = (z2 >> 1) | (z1 << 31);
if (sticky)
z2 |= 1;
z1 = (z1 >> 1) | (z0 << 31);
z0 = (z0 >> 1) | 0x80000000;
ez++;
} else {
if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) {
--- 1158,1252 ----
/* add, propagating carries */
z4 += xy4;
carry = (z4 < xy4);
z3 += xy3;
+
if (carry) {
z3++;
carry = (z3 <= xy3);
! } else {
carry = (z3 < xy3);
+ }
+
z2 += xy2;
+
if (carry) {
z2++;
carry = (z2 <= xy2);
! } else {
carry = (z2 < xy2);
+ }
+
z1 += xy1;
+
if (carry) {
z1++;
carry = (z1 <= xy1);
! } else {
carry = (z1 < xy1);
+ }
+
z0 += xy0;
+
if (carry) {
z0++;
carry = (z0 <= xy0);
! } else {
carry = (z0 < xy0);
+ }
/* for a magnitude subtract, ignore the last carry out */
if (sxy ^ sz)
carry = 0;
/* postnormalize and collect rounding information into z2 */
if (ez < 1) {
/* result is tiny; shift right until exponent is within range */
e = 1 - ez;
+
if (e > 67) {
z2 = 1; /* result can't be exactly zero */
z0 = z1 = 0;
} else if (e >= 64) {
sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e));
z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e));
+
if (sticky)
z2 |= 1;
+
z1 = carry >> (e - 64);
z0 = 0;
} else if (e >= 32) {
sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e));
z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
+
if (sticky)
z2 |= 1;
+
z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e));
z0 = carry >> (e - 32);
} else {
sticky = z4 | z3 | (z2 << 1) << (31 - e);
z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
+
if (sticky)
z2 |= 1;
+
z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
z0 = (z0 >> e) | ((carry << 1) << (31 - e));
}
+
ez = 1;
} else if (carry) {
/* carry out; shift right by one */
sticky = (z2 & 1) | z3 | z4;
z2 = (z2 >> 1) | (z1 << 31);
+
if (sticky)
z2 |= 1;
+
z1 = (z1 >> 1) | (z0 << 31);
z0 = (z0 >> 1) | 0x80000000;
ez++;
} else {
if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) {
*** 1118,1148 ****
--- 1260,1294 ----
z2 = z3;
z3 = z4;
z4 = 0;
ez -= 32;
}
+
while (z0 < 0x80000000u && ez > 1) {
z0 = (z0 << 1) | (z1 >> 31);
z1 = (z1 << 1) | (z2 >> 31);
z2 = (z2 << 1) | (z3 >> 31);
z3 = (z3 << 1) | (z4 >> 31);
z4 <<= 1;
ez--;
}
}
+
if (z3 | z4)
z2 |= 1;
}
/* get the rounding mode */
rm = oldcwsw & 0x0c000000;
/* adjust exponent if result is subnormal */
tinyafter = 0;
+
if (!(z0 & 0x80000000)) {
ez = 0;
tinyafter = 1;
+
if (!(z0 | z1 | z2)) { /* exact zero */
zz.i[2] = rm == FCW_RM ? 0x8000 : 0;
zz.i[1] = zz.i[0] = 0;
__fenv_setcwsw(&oldcwsw);
return (zz.e);
*** 1156,1175 ****
if (sz && (rm == FCW_RM || rm == FCW_RP))
rm = (FCW_RM + FCW_RP) - rm;
/* round */
if (z2) {
! if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u ||
! (z2 == 0x80000000u && (z1 & 1))))) {
/* round up and renormalize if necessary */
if (++z1 == 0) {
if (++z0 == 0) {
z0 = 0x80000000;
ez++;
} else if (z0 == 0x80000000) {
/* rounded up to smallest normal */
ez = 1;
if ((rm == FCW_RP && z2 >
0x80000000u) || (rm == FCW_RN &&
z2 >= 0xc0000000u))
/*
* would have rounded up to
--- 1302,1322 ----
if (sz && (rm == FCW_RM || rm == FCW_RP))
rm = (FCW_RM + FCW_RP) - rm;
/* round */
if (z2) {
! if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u || (z2 ==
! 0x80000000u && (z1 & 1))))) {
/* round up and renormalize if necessary */
if (++z1 == 0) {
if (++z0 == 0) {
z0 = 0x80000000;
ez++;
} else if (z0 == 0x80000000) {
/* rounded up to smallest normal */
ez = 1;
+
if ((rm == FCW_RP && z2 >
0x80000000u) || (rm == FCW_RN &&
z2 >= 0xc0000000u))
/*
* would have rounded up to
*** 1182,1201 ****
--- 1329,1350 ----
}
}
/* restore the control and status words, check for over/underflow */
__fenv_setcwsw(&oldcwsw);
+
if (ez >= 0x7fff) {
if (rm == FCW_RN || rm == FCW_RP) {
zz.i[2] = sz | 0x7fff;
zz.i[1] = 0x80000000;
zz.i[0] = 0;
} else {
zz.i[2] = sz | 0x7ffe;
zz.i[1] = 0xffffffff;
zz.i[0] = 0xffffffff;
}
+
dummy = huge;
dummy *= huge;
} else {
zz.i[2] = sz | ez;
zz.i[1] = z0;
*** 1205,1214 ****
--- 1354,1364 ----
* tinyafter => result rounded w/ unbounded range would be tiny,
* z2 nonzero => result delivered is inexact
*/
if (tinyafter) {
dummy = tiny;
+
if (z2)
dummy *= tiny;
else
dummy -= tiny2;
} else if (z2) {
*** 1217,1225 ****
}
}
return (zz.e);
}
-
#else
#error Unknown architecture
#endif
--- 1367,1374 ----