if (hx >= 0) { /* x > 0 */
if (hx > hy || (hx == hy && lx > ly)) { /* x > y, x -= ulp */
if (lx == 0)
- hx -= 1;
- lx -= 1;
+ hx--;
+ lx--;
} else { /* x < y, x += ulp */
- lx += 1;
+ lx++;
if (lx == 0)
- hx += 1;
+ hx++;
}
} else { /* x < 0 */
if (hy >= 0 || hx > hy || (hx == hy && lx > ly)) { /* x < y, x -= ulp */
if (lx == 0)
- hx -= 1;
- lx -= 1;
+ hx--;
+ lx--;
} else { /* x > y, x += ulp */
- lx += 1;
+ lx++;
if (lx == 0)
- hx += 1;
+ hx++;
}
}
hy = hx & 0x7ff00000;
}
if (hx >= 0) { /* x > 0 */
if (hx > hy) { /* x > y, x -= ulp */
- hx -= 1;
+ hx--;
} else { /* x < y, x += ulp */
- hx += 1;
+ hx++;
}
} else { /* x < 0 */
if (hy >= 0 || hx > hy) { /* x < y, x -= ulp */
- hx -= 1;
+ hx--;
} else { /* x > y, x += ulp */
- hx += 1;
+ hx++;
}
}
hy = hx & 0x7f800000;
if(x > 0.0 ^ x < y) { /* x -= ulp */
if (ux.bits.manl == 0) {
if ((ux.bits.manh&~LDBL_NBIT) == 0)
- ux.bits.exp -= 1;
+ ux.bits.exp--;
ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT);
}
- ux.bits.manl -= 1;
+ ux.bits.manl--;
} else { /* x += ulp */
- ux.bits.manl += 1;
+ ux.bits.manl++;
if (ux.bits.manl == 0) {
ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT);
if ((ux.bits.manh&~LDBL_NBIT)==0)
- ux.bits.exp += 1;
+ ux.bits.exp++;
}
}
if (ux.bits.exp == 0x7fff) /* overflow */
return x;
}
if (hx >= 0 ^ x < y) /* x -= ulp */
- hx -= 1;
+ hx--;
else /* x += ulp */
- hx += 1;
+ hx++;
ix = hx & 0x7f800000;
if (ix >= 0x7f800000) /* overflow */
return x+x;
if (y == 1.0L)
return x;
- // FIXME: this is wrong, see pow special cases in posix2008
+ // FIXME: this is wrong, see pow special cases in c99 F.9.4.4
if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
return y - y; /* +-1**inf is NaN */
if (x == 1.0L)
return (x*p)/(x*p);
if (hx >= 0x7ff00000 || /* x not finite */
(hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */
+ // FIXME: why long double?
return ((long double)x*p)/((long double)x*p);
if (hp <= 0x7fdfffff)
if (hp == 0) /* p = 0 */
return (x*p)/(x*p);
if (hx >= 0x7f800000 || hp > 0x7f800000) /* x not finite, p is NaN */
+ // FIXME: why long double?
return ((long double)x*p)/((long double)x*p);
if (hp <= 0x7effffff)
/* determine ix = ilogb(x) */
if (hx < 0x00100000) { /* subnormal x */
if (hx == 0) {
- for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+ for (ix = -1043, i=lx; i>0; i<<=1) ix--;
} else {
- for (ix = -1022, i=hx<<11; i>0; i<<=1) ix -=1;
+ for (ix = -1022, i=hx<<11; i>0; i<<=1) ix--;
}
} else
ix = (hx>>20) - 1023;
/* determine iy = ilogb(y) */
if (hy < 0x00100000) { /* subnormal y */
if (hy == 0) {
- for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+ for (iy = -1043, i=ly; i>0; i<<=1) iy--;
} else {
- for (iy = -1022, i=hy<<11; i>0; i<<=1) iy -=1;
+ for (iy = -1022, i=hy<<11; i>0; i<<=1) iy--;
}
} else
iy = (hy>>20) - 1023;
while (hx < 0x00100000) { /* normalize x */
hx = hx + hx + (lx>>31);
lx = lx + lx;
- iy -= 1;
+ iy--;
}
if (iy >= -1022) { /* normalize output */
hx = (hx-0x00100000)|((iy+1023)<<20);
if (y < 0x1p-1021) {
if (x + x > y || (x + x == y && (q & 1))) {
q++;
- x-=y;
+ x -= y;
}
} else if (x > 0.5*y || (x == 0.5*y && (q & 1))) {
q++;
- x-=y;
+ x -= y;
}
GET_HIGH_WORD(hx, x);
SET_HIGH_WORD(x, hx ^ sx);
/* determine ix = ilogb(x) */
if (hx < 0x00800000) { /* subnormal x */
- for (ix = -126, i=hx<<8; i>0; i<<=1) ix -=1;
+ for (ix = -126, i=hx<<8; i>0; i<<=1) ix--;
} else
ix = (hx>>23) - 127;
/* determine iy = ilogb(y) */
if (hy < 0x00800000) { /* subnormal y */
- for (iy = -126, i=hy<<8; i>0; i<<=1) iy -=1;
+ for (iy = -126, i=hy<<8; i>0; i<<=1) iy--;
} else
iy = (hy>>23) - 127;
}
while (hx < 0x00800000) { /* normalize x */
hx <<= 1;
- iy -= 1;
+ iy--;
}
if (iy >= -126) { /* normalize output */
hx = (hx-0x00800000)|((iy+127)<<23);
else if (j0 == 18)
i1 = 0x80000000;
else
- i0 = (i0&(~i))|((0x20000)>>j0);
+ i0 = (i0 & ~i)|(0x20000>>j0);
}
}
} else if (j0 > 51) {
return x; /* x is integral */
i >>= 1;
if ((i1&i) != 0)
- i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+ i1 = (i1 & ~i)|(0x40000000>>(j0-20));
}
INSERT_WORDS(x, i0, i1);
STRICT_ASSIGN(double, w, TWO52[sx] + x);
if (k == 0) { /* 0 or subnormal x */
if ((u.bits.manh|u.bits.manl) == 0) /* +-0 */
return x;
- u.e *= 0x1p+128;
+ u.e *= 0x1p128;
k = u.bits.exp - 128;
if (n < -50000)
return tiny*x; /*underflow*/
- }
+ }
if (k == 0x7fff) /* NaN or Inf */
return x + x;
k = k + n;
i1 = 0;
}
} else {
- i = (0x000fffff)>>j0;
+ i = 0x000fffff>>j0;
if (((i0&i)|i1) == 0)
return x; /* x is integral */
/* raise inexact */
return x + x; /* inf or NaN */
return x; /* x is integral */
} else {
- i = ((uint32_t)(0xffffffff))>>(j0-20);
+ i = (uint32_t)0xffffffff>>(j0-20);
if ((i1&i) == 0)
return x; /* x is integral */
/* raise inexact */