summaryrefslogtreecommitdiff
path: root/src/math/remquol.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/math/remquol.c')
-rw-r--r--src/math/remquol.c270
1 files changed, 100 insertions, 170 deletions
diff --git a/src/math/remquol.c b/src/math/remquol.c
index a2e11728..9b065c00 100644
--- a/src/math/remquol.c
+++ b/src/math/remquol.c
@@ -1,15 +1,3 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_remquol.c */
-/*-
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunSoft, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -18,177 +6,119 @@ long double remquol(long double x, long double y, int *quo)
return remquo(x, y, quo);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-
-#define BIAS (LDBL_MAX_EXP - 1)
-
-#if LDBL_MANL_SIZE > 32
-typedef uint64_t manl_t;
-#else
-typedef uint32_t manl_t;
-#endif
-
-#if LDBL_MANH_SIZE > 32
-typedef uint64_t manh_t;
-#else
-typedef uint32_t manh_t;
-#endif
-
-/*
- * These macros add and remove an explicit integer bit in front of the
- * fractional mantissa, if the architecture doesn't have such a bit by
- * default already.
- */
-#ifdef LDBL_IMPLICIT_NBIT
-#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE))
-#define HFRAC_BITS LDBL_MANH_SIZE
-#else
-#define SET_NBIT(hx) (hx)
-#define HFRAC_BITS (LDBL_MANH_SIZE - 1)
-#endif
-
-#define MANL_SHIFT (LDBL_MANL_SIZE - 1)
-
-static const long double Zero[] = {0.0, -0.0};
-
-/*
- * Return the IEEE remainder and set *quo to the last n bits of the
- * quotient, rounded to the nearest integer. We choose n=31 because
- * we wind up computing all the integer bits of the quotient anyway as
- * a side-effect of computing the remainder by the shift and subtract
- * method. In practice, this is far more bits than are needed to use
- * remquo in reduction algorithms.
- *
- * Assumptions:
- * - The low part of the mantissa fits in a manl_t exactly.
- * - The high part of the mantissa fits in an int64_t with enough room
- * for an explicit integer bit in front of the fractional bits.
- */
long double remquol(long double x, long double y, int *quo)
{
- union IEEEl2bits ux, uy;
- int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */
- manh_t hy;
- manl_t lx,ly,lz;
- int ix,iy,n,q,sx,sxy;
-
- ux.e = x;
- uy.e = y;
- sx = ux.bits.sign;
- sxy = sx ^ uy.bits.sign;
- ux.bits.sign = 0; /* |x| */
- uy.bits.sign = 0; /* |y| */
- x = ux.e;
-
- /* purge off exception values */
- if ((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */
- (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
- (uy.bits.exp == BIAS + LDBL_MAX_EXP &&
- ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */
+ union ldshape ux = {x}, uy = {y};
+ int ex = ux.i.se & 0x7fff;
+ int ey = uy.i.se & 0x7fff;
+ int sx = ux.i.se >> 15;
+ int sy = uy.i.se >> 15;
+ uint32_t q;
+
+ *quo = 0;
+ if (y == 0 || isnan(y) || ex == 0x7fff)
return (x*y)/(x*y);
- if (ux.bits.exp <= uy.bits.exp) {
- if ((ux.bits.exp < uy.bits.exp) ||
- (ux.bits.manh <= uy.bits.manh &&
- (ux.bits.manh < uy.bits.manh ||
- ux.bits.manl < uy.bits.manl))) {
- q = 0;
- goto fixup; /* |x|<|y| return x or x-y */
- }
- if (ux.bits.manh == uy.bits.manh && ux.bits.manl == uy.bits.manl) {
- *quo = sxy ? -1 : 1;
- return Zero[sx]; /* |x|=|y| return x*0*/
- }
- }
-
- /* determine ix = ilogb(x) */
- if (ux.bits.exp == 0) { /* subnormal x */
- ux.e *= 0x1.0p512;
- ix = ux.bits.exp - (BIAS + 512);
- } else {
- ix = ux.bits.exp - BIAS;
+ if (x == 0)
+ return x;
+
+ /* normalize x and y */
+ if (!ex) {
+ ux.i.se = ex;
+ ux.f *= 0x1p120f;
+ ex = ux.i.se - 120;
}
-
- /* determine iy = ilogb(y) */
- if (uy.bits.exp == 0) { /* subnormal y */
- uy.e *= 0x1.0p512;
- iy = uy.bits.exp - (BIAS + 512);
- } else {
- iy = uy.bits.exp - BIAS;
+ if (!ey) {
+ uy.i.se = ey;
+ uy.f *= 0x1p120f;
+ ey = uy.i.se - 120;
}
- /* set up {hx,lx}, {hy,ly} and align y to x */
- hx = SET_NBIT(ux.bits.manh);
- hy = SET_NBIT(uy.bits.manh);
- lx = ux.bits.manl;
- ly = uy.bits.manl;
-
- /* fix point fmod */
- n = ix - iy;
q = 0;
-
- while (n--) {
- hz = hx - hy;
- lz = lx - ly;
- if (lx < ly)
- hz -= 1;
- if (hz < 0) {
- hx = hx + hx + (lx>>MANL_SHIFT);
- lx = lx + lx;
- } else {
- hx = hz + hz + (lz>>MANL_SHIFT);
- lx = lz + lz;
+ if (ex >= ey) {
+ /* x mod y */
+#if LDBL_MANT_DIG == 64
+ uint64_t i, mx, my;
+ mx = ux.i.m;
+ my = uy.i.m;
+ for (; ex > ey; ex--) {
+ i = mx - my;
+ if (mx >= my) {
+ mx = 2*i;
+ q++;
+ q <<= 1;
+ } else if (2*mx < mx) {
+ mx = 2*mx - my;
+ q <<= 1;
+ q++;
+ } else {
+ mx = 2*mx;
+ q <<= 1;
+ }
+ }
+ i = mx - my;
+ if (mx >= my) {
+ mx = i;
q++;
}
- q <<= 1;
- }
- hz = hx - hy;
- lz = lx - ly;
- if (lx < ly)
- hz -= 1;
- if (hz >= 0) {
- hx = hz;
- lx = lz;
- q++;
- }
-
- /* convert back to floating value and restore the sign */
- if ((hx|lx) == 0) { /* return sign(x)*0 */
- q &= 0x7fffffff;
- *quo = sxy ? -q : q;
- return Zero[sx];
- }
- while (hx < (1ULL<<HFRAC_BITS)) { /* normalize x */
- hx = hx + hx + (lx>>MANL_SHIFT);
- lx = lx + lx;
- iy -= 1;
- }
- ux.bits.manh = hx; /* The integer bit is truncated here if needed. */
- ux.bits.manl = lx;
- if (iy < LDBL_MIN_EXP) {
- ux.bits.exp = iy + (BIAS + 512);
- ux.e *= 0x1p-512;
- } else {
- ux.bits.exp = iy + BIAS;
- }
- ux.bits.sign = 0;
- x = ux.e;
-fixup:
- y = fabsl(y);
- if (y < LDBL_MIN * 2) {
- if (x + x > y || (x + x == y && (q & 1))) {
+ if (mx == 0)
+ ex = -120;
+ else
+ for (; mx >> 63 == 0; mx *= 2, ex--);
+ ux.i.m = mx;
+#elif LDBL_MANT_DIG == 113
+ uint64_t hi, lo, xhi, xlo, yhi, ylo;
+ xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
+ yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
+ xlo = ux.i2.lo;
+ ylo = ux.i2.lo;
+ for (; ex > ey; ex--) {
+ hi = xhi - yhi;
+ lo = xlo - ylo;
+ if (xlo < ylo)
+ hi -= 1;
+ if (hi >> 63 == 0) {
+ xhi = 2*hi + (lo>>63);
+ xlo = 2*lo;
+ q++;
+ } else {
+ xhi = 2*xhi + (xlo>>63);
+ xlo = 2*xlo;
+ }
+ q <<= 1;
+ }
+ hi = xhi - yhi;
+ lo = xlo - ylo;
+ if (xlo < ylo)
+ hi -= 1;
+ if (hi >> 63 == 0) {
+ xhi = hi;
+ xlo = lo;
q++;
- x-=y;
}
- } else if (x > 0.5*y || (x == 0.5*y && (q & 1))) {
- q++;
- x-=y;
+ if ((xhi|xlo) == 0)
+ ex = -120;
+ else
+ for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
+ ux.i2.hi = xhi;
+ ux.i2.lo = xlo;
+#endif
}
- ux.e = x;
- ux.bits.sign ^= sx;
- x = ux.e;
-
+ /* scale result and decide between |x| and |x|-|y| */
+ if (ex <= 0) {
+ ux.i.se = ex + 120;
+ ux.f *= 0x1p-120f;
+ } else
+ ux.i.se = ex;
+ x = ux.f;
+ if (sy)
+ y = -y;
+ if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
+ x -= y;
+ q++;
+ }
q &= 0x7fffffff;
- *quo = sxy ? -q : q;
- return x;
+ *quo = sx^sy ? -(int)q : (int)q;
+ return sx ? -x : x;
}
#endif
generated by cgit v1.2.1 (git 2.18.0) at 2024-04-24 06:23:15 +0000