-/* origin: FreeBSD /usr/src/lib/msun/src/e_fmod.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.
- * ====================================================
- */
-/*
- * fmod(x,y)
- * Return x mod y in exact arithmetic
- * Method: shift and subtract
- */
-
-#include "libm.h"
-
-static const double Zero[] = {0.0, -0.0,};
+#include <math.h>
+#include <stdint.h>
double fmod(double x, double y)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- uint32_t lx,ly,lz;
+ union {double f; uint64_t i;} ux = {x}, uy = {y};
+ int ex = ux.i>>52 & 0x7ff;
+ int ey = uy.i>>52 & 0x7ff;
+ int sx = ux.i>>63;
+ uint64_t i;
- EXTRACT_WORDS(hx, lx, x);
- EXTRACT_WORDS(hy, ly, y);
- sx = hx & 0x80000000; /* sign of x */
- hx ^= sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
+ /* in the followings uxi should be ux.i, but then gcc wrongly adds */
+ /* float load/store to inner loops ruining performance and code size */
+ uint64_t uxi = ux.i;
- /* purge off exception values */
- if ((hy|ly) == 0 || hx >= 0x7ff00000 || /* y=0,or x not finite */
- (hy|((ly|-ly)>>31)) > 0x7ff00000) /* or y is NaN */
+ if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
return (x*y)/(x*y);
- if (hx <= hy) {
- if (hx < hy || lx < ly) /* |x| < |y| */
- return x;
- if (lx == ly) /* |x| = |y|, return x*0 */
- return Zero[(uint32_t)sx>>31];
+ if (uxi<<1 <= uy.i<<1) {
+ if (uxi<<1 == uy.i<<1)
+ return 0*x;
+ return x;
}
- /* determine ix = ilogb(x) */
- if (hx < 0x00100000) { /* subnormal x */
- if (hx == 0) {
- for (ix = -1043, i = lx; i > 0; i <<= 1)
- ix -= 1;
- } else {
- for (ix = -1022, i = hx<<11; i > 0; i <<= 1)
- ix -= 1;
- }
- } 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;
- } else {
- for (iy = -1022, i = hy<<11; i > 0; i <<= 1)
- iy -= 1;
- }
- } else
- iy = (hy>>20) - 1023;
-
- /* set up {hx,lx}, {hy,ly} and align y to x */
- if (ix >= -1022)
- hx = 0x00100000|(0x000fffff&hx);
- else { /* subnormal x, shift x to normal */
- n = -1022-ix;
- if (n <= 31) {
- hx = (hx<<n)|(lx>>(32-n));
- lx <<= n;
- } else {
- hx = lx<<(n-32);
- lx = 0;
- }
+ /* normalize x and y */
+ if (!ex) {
+ for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
+ uxi <<= -ex + 1;
+ } else {
+ uxi &= -1ULL >> 12;
+ uxi |= 1ULL << 52;
}
- if(iy >= -1022)
- hy = 0x00100000|(0x000fffff&hy);
- else { /* subnormal y, shift y to normal */
- n = -1022-iy;
- if (n <= 31) {
- hy = (hy<<n)|(ly>>(32-n));
- ly <<= n;
- } else {
- hy = ly<<(n-32);
- ly = 0;
- }
+ if (!ey) {
+ for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
+ uy.i <<= -ey + 1;
+ } else {
+ uy.i &= -1ULL >> 12;
+ uy.i |= 1ULL << 52;
}
- /* fix point fmod */
- n = ix - iy;
- while (n--) {
- hz = hx-hy;
- lz = lx-ly;
- if (lx < ly)
- hz -= 1;
- if (hz < 0) {
- hx = hx+hx+(lx>>31);
- lx = lx+lx;
- } else {
- if ((hz|lz) == 0) /* return sign(x)*0 */
- return Zero[(uint32_t)sx>>31];
- hx = hz+hz+(lz>>31);
- lx = lz+lz;
+ /* x mod y */
+ for (; ex > ey; ex--) {
+ i = uxi - uy.i;
+ if (i >> 63 == 0) {
+ if (i == 0)
+ return 0*x;
+ uxi = i;
}
+ uxi <<= 1;
}
- hz = hx-hy;
- lz = lx-ly;
- if (lx < ly)
- hz -= 1;
- if (hz >= 0) {
- hx = hz;
- lx = lz;
+ i = uxi - uy.i;
+ if (i >> 63 == 0) {
+ if (i == 0)
+ return 0*x;
+ uxi = i;
}
+ for (; uxi>>52 == 0; uxi <<= 1, ex--);
- /* convert back to floating value and restore the sign */
- if ((hx|lx) == 0) /* return sign(x)*0 */
- return Zero[(uint32_t)sx>>31];
- while (hx < 0x00100000) { /* normalize x */
- hx = hx+hx+(lx>>31);
- lx = lx+lx;
- iy -= 1;
- }
- if (iy >= -1022) { /* normalize output */
- hx = ((hx-0x00100000)|((iy+1023)<<20));
- INSERT_WORDS(x, hx|sx, lx);
- } else { /* subnormal output */
- n = -1022 - iy;
- if (n <= 20) {
- lx = (lx>>n)|((uint32_t)hx<<(32-n));
- hx >>= n;
- } else if (n <= 31) {
- lx = (hx<<(32-n))|(lx>>n);
- hx = sx;
- } else {
- lx = hx>>(n-32); hx = sx;
- }
- INSERT_WORDS(x, hx|sx, lx);
+ /* scale result */
+ if (ex > 0) {
+ uxi -= 1ULL << 52;
+ uxi |= (uint64_t)ex << 52;
+ } else {
+ uxi >>= -ex + 1;
}
- return x; /* exact output */
+ uxi |= (uint64_t)sx << 63;
+ ux.i = uxi;
+ return ux.f;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/e_fmodf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/*
- * fmodf(x,y)
- * Return x mod y in exact arithmetic
- * Method: shift and subtract
- */
-
-#include "libm.h"
-
-static const float Zero[] = {0.0, -0.0,};
+#include <math.h>
+#include <stdint.h>
float fmodf(float x, float y)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
+ union {float f; uint32_t i;} ux = {x}, uy = {y};
+ int ex = ux.i>>23 & 0xff;
+ int ey = uy.i>>23 & 0xff;
+ uint32_t sx = ux.i & 0x80000000;
+ uint32_t i;
+ uint32_t uxi = ux.i;
- GET_FLOAT_WORD(hx, x);
- GET_FLOAT_WORD(hy, y);
- sx = hx & 0x80000000; /* sign of x */
- hx ^= sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
-
- /* purge off exception values */
- if (hy == 0 || hx >= 0x7f800000 || /* y=0,or x not finite */
- hy > 0x7f800000) /* or y is NaN */
+ if (uy.i<<1 == 0 || isnan(y) || ex == 0xff)
return (x*y)/(x*y);
- if (hx < hy) /* |x| < |y| */
+ if (uxi<<1 <= uy.i<<1) {
+ if (uxi<<1 == uy.i<<1)
+ return 0*x;
return x;
- if (hx == hy) /* |x| = |y|, return x*0 */
- return Zero[(uint32_t)sx>>31];
-
- /* determine ix = ilogb(x) */
- if (hx < 0x00800000) { /* subnormal x */
- for (ix = -126, i = hx<<8; i > 0; i <<= 1)
- ix -= 1;
- } 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;
- } else
- iy = (hy>>23) - 127;
+ }
- /* set up {hx,lx}, {hy,ly} and align y to x */
- if (ix >= -126)
- hx = 0x00800000|(0x007fffff&hx);
- else { /* subnormal x, shift x to normal */
- n = -126-ix;
- hx = hx<<n;
+ /* normalize x and y */
+ if (!ex) {
+ for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1);
+ uxi <<= -ex + 1;
+ } else {
+ uxi &= -1U >> 9;
+ uxi |= 1U << 23;
}
- if (iy >= -126)
- hy = 0x00800000|(0x007fffff&hy);
- else { /* subnormal y, shift y to normal */
- n = -126-iy;
- hy = hy<<n;
+ if (!ey) {
+ for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1);
+ uy.i <<= -ey + 1;
+ } else {
+ uy.i &= -1U >> 9;
+ uy.i |= 1U << 23;
}
- /* fix point fmod */
- n = ix - iy;
- while (n--) {
- hz = hx-hy;
- if (hz<0)
- hx = hx+hx;
- else {
- if(hz == 0) /* return sign(x)*0 */
- return Zero[(uint32_t)sx>>31];
- hx = hz+hz;
+ /* x mod y */
+ for (; ex > ey; ex--) {
+ i = uxi - uy.i;
+ if (i >> 31 == 0) {
+ if (i == 0)
+ return 0*x;
+ uxi = i;
}
+ uxi <<= 1;
}
- hz = hx-hy;
- if (hz >= 0)
- hx = hz;
-
- /* convert back to floating value and restore the sign */
- if (hx == 0) /* return sign(x)*0 */
- return Zero[(uint32_t)sx>>31];
- while (hx < 0x00800000) { /* normalize x */
- hx = hx+hx;
- iy -= 1;
+ i = uxi - uy.i;
+ if (i >> 31 == 0) {
+ if (i == 0)
+ return 0*x;
+ uxi = i;
}
- if (iy >= -126) { /* normalize output */
- hx = ((hx-0x00800000)|((iy+127)<<23));
- SET_FLOAT_WORD(x, hx|sx);
- } else { /* subnormal output */
- n = -126 - iy;
- hx >>= n;
- SET_FLOAT_WORD(x, hx|sx);
+ for (; uxi>>23 == 0; uxi <<= 1, ex--);
+
+ /* scale result up */
+ if (ex > 0) {
+ uxi -= 1U << 23;
+ uxi |= (uint32_t)ex << 23;
+ } else {
+ uxi >>= -ex + 1;
}
- return x; /* exact output */
+ uxi |= sx;
+ ux.i = uxi;
+ return ux.f;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/e_fmodl.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
return fmod(x, y);
}
#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,};
-
-/*
- * fmodl(x,y)
- * Return x mod y in exact arithmetic
- * Method: shift and subtract
- *
- * 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 fmodl(long double x, long double y)
{
- 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,sx;
-
- ux.e = x;
- uy.e = y;
- sx = ux.bits.sign;
+ union ldshape ux = {x}, uy = {y};
+ int ex = ux.i.se & 0x7fff;
+ int ey = uy.i.se & 0x7fff;
+ int sx = ux.i.se & 0x8000;
- /* 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 */
+ 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))) /* |x|<|y| return x or x-y */
- return x;
- if (ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl)
- return Zero[sx]; /* |x| = |y| return x*0 */
+ ux.i.se = ex;
+ uy.i.se = ey;
+ if (ux.f <= uy.f) {
+ if (ux.f == uy.f)
+ return 0*x;
+ return x;
}
- /* 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;
+ /* normalize x and y */
+ if (!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.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;
-
- 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;
+ /* 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) {
+ if (i == 0)
+ return 0*x;
+ mx = 2*i;
+ } else if (2*mx < mx) {
+ mx = 2*mx - my;
} else {
- if ((hz|lz)==0) /* return sign(x)*0 */
- return Zero[sx];
- hx = hz+hz+(lz>>MANL_SHIFT);
- lx = lz+lz;
+ mx = 2*mx;
}
}
- hz = hx-hy;
- lz = lx-ly;
- if (lx < ly)
- hz -= 1;
- if (hz >= 0) {
- hx = hz;
- lx = lz;
+ i = mx - my;
+ if (mx >= my) {
+ if (i == 0)
+ return 0*x;
+ mx = i;
}
-
- /* convert back to floating value and restore the sign */
- if ((hx|lx) == 0) /* return sign(x)*0 */
- return Zero[sx];
- while (hx < (1ULL<<HFRAC_BITS)) { /* normalize x */
- hx = hx+hx+(lx>>MANL_SHIFT);
- lx = lx+lx;
- iy -= 1;
+ 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) {
+ if ((hi|lo) == 0)
+ return 0*x;
+ xhi = 2*hi + (lo>>63);
+ xlo = 2*lo;
+ } else {
+ xhi = 2*xhi + (xlo>>63);
+ xlo = 2*xlo;
+ }
}
- ux.bits.manh = hx; /* The mantissa 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;
+ hi = xhi - yhi;
+ lo = xlo - ylo;
+ if (xlo < ylo)
+ hi -= 1;
+ if (hi >> 63 == 0) {
+ if ((hi|lo) == 0)
+ return 0*x;
+ xhi = hi;
+ xlo = lo;
}
- return ux.e; /* exact output */
+ for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
+ ux.i2.hi = xhi;
+ ux.i2.lo = xlo;
+#endif
+
+ /* scale result */
+ if (ex <= 0) {
+ ux.i.se = (ex+120)|sx;
+ ux.f *= 0x1p-120f;
+ } else
+ ux.i.se = ex|sx;
+ return ux.f;
}
#endif
double modf(double x, double *iptr)
{
- union {double x; uint64_t n;} u = {x};
+ union {double f; uint64_t i;} u = {x};
uint64_t mask;
- int e;
-
- e = (int)(u.n>>52 & 0x7ff) - 0x3ff;
+ int e = (int)(u.i>>52 & 0x7ff) - 0x3ff;
/* no fractional part */
if (e >= 52) {
*iptr = x;
- if (e == 0x400 && u.n<<12 != 0) /* nan */
+ if (e == 0x400 && u.i<<12 != 0) /* nan */
return x;
- u.n &= (uint64_t)1<<63;
- return u.x;
+ u.i &= 1ULL<<63;
+ return u.f;
}
/* no integral part*/
if (e < 0) {
- u.n &= (uint64_t)1<<63;
- *iptr = u.x;
+ u.i &= 1ULL<<63;
+ *iptr = u.f;
return x;
}
- mask = (uint64_t)-1>>12 >> e;
- if ((u.n & mask) == 0) {
+ mask = -1ULL>>12>>e;
+ if ((u.i & mask) == 0) {
*iptr = x;
- u.n &= (uint64_t)1<<63;
- return u.x;
+ u.i &= 1ULL<<63;
+ return u.f;
}
- u.n &= ~mask;
- *iptr = u.x;
- STRICT_ASSIGN(double, x, x - *iptr);
- return x;
+ u.i &= ~mask;
+ *iptr = u.f;
+ return x - u.f;
}
float modff(float x, float *iptr)
{
- union {float x; uint32_t n;} u = {x};
+ union {float f; uint32_t i;} u = {x};
uint32_t mask;
- int e;
-
- e = (int)(u.n>>23 & 0xff) - 0x7f;
+ int e = (int)(u.i>>23 & 0xff) - 0x7f;
/* no fractional part */
if (e >= 23) {
*iptr = x;
- if (e == 0x80 && u.n<<9 != 0) { /* nan */
+ if (e == 0x80 && u.i<<9 != 0) { /* nan */
return x;
}
- u.n &= 0x80000000;
- return u.x;
+ u.i &= 0x80000000;
+ return u.f;
}
/* no integral part */
if (e < 0) {
- u.n &= 0x80000000;
- *iptr = u.x;
+ u.i &= 0x80000000;
+ *iptr = u.f;
return x;
}
mask = 0x007fffff>>e;
- if ((u.n & mask) == 0) {
+ if ((u.i & mask) == 0) {
*iptr = x;
- u.n &= 0x80000000;
- return u.x;
+ u.i &= 0x80000000;
+ return u.f;
}
- u.n &= ~mask;
- *iptr = u.x;
- STRICT_ASSIGN(float, x, x - *iptr);
- return x;
+ u.i &= ~mask;
+ *iptr = u.f;
+ return x - u.f;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/s_modfl.c */
-/*-
- * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * Derived from s_modf.c, which has the following Copyright:
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, 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
return modf(x, (double *)iptr);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
-
-#if LDBL_MANL_SIZE > 32
-#define MASK ((uint64_t)-1)
-#else
-#define MASK ((uint32_t)-1)
+#if LDBL_MANT_DIG == 64
+#define TOINT 0x1p63
+#elif LDBL_MANT_DIG == 113
+#define TOINT 0x1p112
#endif
-/* Return the last n bits of a word, representing the fractional part. */
-#define GETFRAC(bits, n) ((bits) & ~(MASK << (n)))
-/* The number of fraction bits in manh, not counting the integer bit */
-#define HIBITS (LDBL_MANT_DIG - LDBL_MANL_SIZE)
-
-static const long double zero[] = { 0.0, -0.0 };
-
long double modfl(long double x, long double *iptr)
{
- union IEEEl2bits ux;
- int e;
+ union ldshape u = {x};
+ uint64_t mask;
+ int e = (u.i.se & 0x7fff) - 0x3fff;
+ int s = u.i.se >> 15;
+ long double absx;
+ long double y;
- ux.e = x;
- e = ux.bits.exp - LDBL_MAX_EXP + 1;
- if (e < HIBITS) { /* Integer part is in manh. */
- if (e < 0) { /* |x|<1 */
- *iptr = zero[ux.bits.sign];
- return x;
- }
- if ((GETFRAC(ux.bits.manh, HIBITS - 1 - e)|ux.bits.manl) == 0) {
- /* x is an integer. */
- *iptr = x;
- return zero[ux.bits.sign];
- }
- /* Clear all but the top e+1 bits. */
- ux.bits.manh >>= HIBITS - 1 - e;
- ux.bits.manh <<= HIBITS - 1 - e;
- ux.bits.manl = 0;
- *iptr = ux.e;
- return x - ux.e;
- } else if (e >= LDBL_MANT_DIG - 1) { /* x has no fraction part. */
+ /* no fractional part */
+ if (e >= LDBL_MANT_DIG-1) {
*iptr = x;
- if (e == LDBL_MAX_EXP && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)) /* nan */
+ if (isnan(x))
return x;
- return zero[ux.bits.sign];
- } else { /* Fraction part is in manl. */
- if (GETFRAC(ux.bits.manl, LDBL_MANT_DIG - 1 - e) == 0) {
- /* x is integral. */
- *iptr = x;
- return zero[ux.bits.sign];
- }
- /* Clear all but the top e+1 bits. */
- ux.bits.manl >>= LDBL_MANT_DIG - 1 - e;
- ux.bits.manl <<= LDBL_MANT_DIG - 1 - e;
- *iptr = ux.e;
- return x - ux.e;
+ return s ? -0.0 : 0.0;
+ }
+
+ /* no integral part*/
+ if (e < 0) {
+ *iptr = s ? -0.0 : 0.0;
+ return x;
+ }
+
+ /* raises spurious inexact */
+ absx = s ? -x : x;
+ y = absx + TOINT - TOINT - absx;
+ if (y == 0) {
+ *iptr = x;
+ return s ? -0.0 : 0.0;
}
+ if (y > 0)
+ y -= 1;
+ if (s)
+ y = -y;
+ *iptr = x + y;
+ return -y;
}
#endif
-/* origin: FreeBSD /usr/src/lib/msun/src/e_remainder.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.
- * ====================================================
- */
-/* remainder(x,p)
- * Return :
- * returns x REM p = x - [x/p]*p as if in infinite
- * precise arithmetic, where [x/p] is the (infinite bit)
- * integer nearest x/p (in half way case choose the even one).
- * Method :
- * Based on fmod() return x-[x/p]chopped*p exactlp.
- */
+#include <math.h>
-#include "libm.h"
-
-double remainder(double x, double p)
+double remainder(double x, double y)
{
- int32_t hx,hp;
- uint32_t sx,lx,lp;
- double p_half;
-
- EXTRACT_WORDS(hx, lx, x);
- EXTRACT_WORDS(hp, lp, p);
- sx = hx & 0x80000000;
- hp &= 0x7fffffff;
- hx &= 0x7fffffff;
-
- /* purge off exception values */
- if ((hp|lp) == 0 || /* p = 0 */
- hx >= 0x7ff00000 || /* x not finite */
- (hp >= 0x7ff00000 && (hp-0x7ff00000 | lp) != 0)) /* p is NaN */
- return (x*p)/(x*p);
-
- if (hp <= 0x7fdfffff)
- x = fmod(x, p+p); /* now x < 2p */
- if (((hx-hp)|(lx-lp)) == 0)
- return 0.0*x;
- x = fabs(x);
- p = fabs(p);
- if (hp < 0x00200000) {
- if (x + x > p) {
- x -= p;
- if (x + x >= p)
- x -= p;
- }
- } else {
- p_half = 0.5*p;
- if (x > p_half) {
- x -= p;
- if (x >= p_half)
- x -= p;
- }
- }
- GET_HIGH_WORD(hx, x);
- if ((hx&0x7fffffff) == 0)
- hx = 0;
- SET_HIGH_WORD(x, hx^sx);
- return x;
+ int q;
+ return remquo(x, y, &q);
}
-/* origin: FreeBSD /usr/src/lib/msun/src/e_remainderf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+#include <math.h>
-#include "libm.h"
-
-float remainderf(float x, float p)
+float remainderf(float x, float y)
{
- int32_t hx,hp;
- uint32_t sx;
- float p_half;
-
- GET_FLOAT_WORD(hx, x);
- GET_FLOAT_WORD(hp, p);
- sx = hx & 0x80000000;
- hp &= 0x7fffffff;
- hx &= 0x7fffffff;
-
- /* purge off exception values */
- if (hp == 0 || hx >= 0x7f800000 || hp > 0x7f800000) /* p = 0, x not finite, p is NaN */
- return (x*p)/(x*p);
-
- if (hp <= 0x7effffff)
- x = fmodf(x, p + p); /* now x < 2p */
- if (hx - hp == 0)
- return 0.0f*x;
- x = fabsf(x);
- p = fabsf(p);
- if (hp < 0x01000000) {
- if (x + x > p) {
- x -= p;
- if (x + x >= p)
- x -= p;
- }
- } else {
- p_half = 0.5f*p;
- if (x > p_half) {
- x -= p;
- if (x >= p_half)
- x -= p;
- }
- }
- GET_FLOAT_WORD(hx, x);
- if ((hx & 0x7fffffff) == 0)
- hx = 0;
- SET_FLOAT_WORD(x, hx ^ sx);
- return x;
+ int q;
+ return remquof(x, y, &q);
}
-/* origin: FreeBSD /usr/src/lib/msun/src/s_remquo.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.
- * ====================================================
- */
-/*
- * 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.
- */
-
-#include "libm.h"
-
-static const double Zero[] = {0.0, -0.0,};
+#include <math.h>
+#include <stdint.h>
double remquo(double x, double y, int *quo)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- uint32_t lx,ly,lz,q,sxy;
-
- EXTRACT_WORDS(hx, lx, x);
- EXTRACT_WORDS(hy, ly, y);
- sxy = (hx ^ hy) & 0x80000000;
- sx = hx & 0x80000000; /* sign of x */
- hx ^= sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
+ union {double f; uint64_t i;} ux = {x}, uy = {y};
+ int ex = ux.i>>52 & 0x7ff;
+ int ey = uy.i>>52 & 0x7ff;
+ int sx = ux.i>>63;
+ int sy = uy.i>>63;
+ uint32_t q;
+ uint64_t i;
+ uint64_t uxi = ux.i;
- /* purge off exception values */
- if ((hy|ly) == 0 || hx >= 0x7ff00000 || /* y = 0, or x not finite */
- (hy|((ly|-ly)>>31)) > 0x7ff00000) /* or y is NaN */
+ *quo = 0;
+ if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
return (x*y)/(x*y);
- if (hx <= hy) {
- if (hx < hy || lx < ly) { /* |x| < |y| return x or x-y */
- q = 0;
- goto fixup;
- }
- if (lx == ly) { /* |x| = |y| return x*0 */
- *quo = sxy ? -1 : 1;
- return Zero[(uint32_t)sx>>31];
- }
- }
+ if (ux.i<<1 == 0)
+ return x;
- // FIXME: use ilogb?
-
- /* determine ix = ilogb(x) */
- if (hx < 0x00100000) { /* subnormal x */
- if (hx == 0) {
- for (ix = -1043, i=lx; i>0; i<<=1) ix--;
- } else {
- 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--;
- } else {
- for (iy = -1022, i=hy<<11; i>0; i<<=1) iy--;
- }
- } else
- iy = (hy>>20) - 1023;
-
- /* set up {hx,lx}, {hy,ly} and align y to x */
- if (ix >= -1022)
- hx = 0x00100000|(0x000fffff&hx);
- else { /* subnormal x, shift x to normal */
- n = -1022 - ix;
- if (n <= 31) {
- hx = (hx<<n)|(lx>>(32-n));
- lx <<= n;
- } else {
- hx = lx<<(n-32);
- lx = 0;
- }
+ /* normalize x and y */
+ if (!ex) {
+ for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
+ uxi <<= -ex + 1;
+ } else {
+ uxi &= -1ULL >> 12;
+ uxi |= 1ULL << 52;
}
- if (iy >= -1022)
- hy = 0x00100000|(0x000fffff&hy);
- else { /* subnormal y, shift y to normal */
- n = -1022 - iy;
- if (n <= 31) {
- hy = (hy<<n)|(ly>>(32-n));
- ly <<= n;
- } else {
- hy = ly<<(n-32);
- ly = 0;
- }
+ if (!ey) {
+ for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
+ uy.i <<= -ey + 1;
+ } else {
+ uy.i &= -1ULL >> 12;
+ uy.i |= 1ULL << 52;
}
- /* fix point fmod */
- n = ix - iy;
q = 0;
- while (n--) {
- hz = hx - hy;
- lz = lx - ly;
- if (lx < ly)
- hz--;
- if (hz < 0) {
- hx = hx + hx + (lx>>31);
- lx = lx + lx;
- } else {
- hx = hz + hz + (lz>>31);
- lx = lz + lz;
+ if (ex < ey) {
+ if (ex+1 == ey)
+ goto end;
+ return x;
+ }
+
+ /* x mod y */
+ for (; ex > ey; ex--) {
+ i = uxi - uy.i;
+ if (i >> 63 == 0) {
+ uxi = i;
q++;
}
+ uxi <<= 1;
q <<= 1;
}
- hz = hx - hy;
- lz = lx - ly;
- if (lx < ly)
- hz--;
- if (hz >= 0) {
- hx = hz;
- lx = lz;
+ i = uxi - uy.i;
+ if (i >> 63 == 0) {
+ uxi = i;
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[(uint32_t)sx>>31];
- }
- while (hx < 0x00100000) { /* normalize x */
- hx = hx + hx + (lx>>31);
- lx = lx + lx;
- iy--;
+ if (uxi == 0)
+ ex = -60;
+ else
+ for (; uxi>>52 == 0; uxi <<= 1, ex--);
+end:
+ /* scale result and decide between |x| and |x|-|y| */
+ if (ex > 0) {
+ uxi -= 1ULL << 52;
+ uxi |= (uint64_t)ex << 52;
+ } else {
+ uxi >>= -ex + 1;
}
- if (iy >= -1022) { /* normalize output */
- hx = (hx-0x00100000)|((iy+1023)<<20);
- } else { /* subnormal output */
- n = -1022 - iy;
- if (n <= 20) {
- lx = (lx>>n)|((uint32_t)hx<<(32-n));
- hx >>= n;
- } else if (n <= 31) {
- lx = (hx<<(32-n))|(lx>>n);
- hx = 0;
- } else {
- lx = hx>>(n-32);
- hx = 0;
- }
- }
-fixup:
- INSERT_WORDS(x, hx, lx);
- y = fabs(y);
- if (y < 0x1p-1021) {
- if (x + x > y || (x + x == y && (q & 1))) {
- q++;
- x -= y;
- }
- } else if (x > 0.5*y || (x == 0.5*y && (q & 1))) {
- q++;
+ ux.i = uxi;
+ x = ux.f;
+ if (sy)
+ y = -y;
+ if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
x -= y;
+ q++;
}
- GET_HIGH_WORD(hx, x);
- SET_HIGH_WORD(x, hx ^ sx);
q &= 0x7fffffff;
- *quo = sxy ? -q : q;
- return x;
+ *quo = sx^sy ? -(int)q : (int)q;
+ return sx ? -x : x;
}
-/* origin: FreeBSD /usr/src/lib/msun/src/s_remquof.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.
- * ====================================================
- */
-/*
- * 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.
- */
-
-#include "libm.h"
-
-static const float Zero[] = {0.0, -0.0,};
+#include <math.h>
+#include <stdint.h>
float remquof(float x, float y, int *quo)
{
- int32_t n,hx,hy,hz,ix,iy,sx,i;
- uint32_t q,sxy;
+ union {float f; uint32_t i;} ux = {x}, uy = {y};
+ int ex = ux.i>>23 & 0xff;
+ int ey = uy.i>>23 & 0xff;
+ int sx = ux.i>>31;
+ int sy = uy.i>>31;
+ uint32_t q;
+ uint32_t i;
+ uint32_t uxi = ux.i;
- GET_FLOAT_WORD(hx, x);
- GET_FLOAT_WORD(hy, y);
- sxy = (hx ^ hy) & 0x80000000;
- sx = hx & 0x80000000; /* sign of x */
- hx ^= sx; /* |x| */
- hy &= 0x7fffffff; /* |y| */
-
- /* purge off exception values */
- if (hy == 0 || hx >= 0x7f800000 || hy > 0x7f800000) /* y=0,NaN;or x not finite */
+ *quo = 0;
+ if (uy.i<<1 == 0 || isnan(y) || ex == 0xff)
return (x*y)/(x*y);
- if (hx < hy) { /* |x| < |y| return x or x-y */
- q = 0;
- goto fixup;
- } else if(hx==hy) { /* |x| = |y| return x*0*/
- *quo = sxy ? -1 : 1;
- return Zero[(uint32_t)sx>>31];
- }
+ if (ux.i<<1 == 0)
+ return x;
- /* determine ix = ilogb(x) */
- if (hx < 0x00800000) { /* subnormal x */
- 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--;
- } else
- iy = (hy>>23) - 127;
-
- /* set up {hx,lx}, {hy,ly} and align y to x */
- if (ix >= -126)
- hx = 0x00800000|(0x007fffff&hx);
- else { /* subnormal x, shift x to normal */
- n = -126 - ix;
- hx <<= n;
+ /* normalize x and y */
+ if (!ex) {
+ for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1);
+ uxi <<= -ex + 1;
+ } else {
+ uxi &= -1U >> 9;
+ uxi |= 1U << 23;
}
- if (iy >= -126)
- hy = 0x00800000|(0x007fffff&hy);
- else { /* subnormal y, shift y to normal */
- n = -126 - iy;
- hy <<= n;
+ if (!ey) {
+ for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1);
+ uy.i <<= -ey + 1;
+ } else {
+ uy.i &= -1U >> 9;
+ uy.i |= 1U << 23;
}
- /* fix point fmod */
- n = ix - iy;
q = 0;
- while (n--) {
- hz = hx - hy;
- if (hz < 0)
- hx = hx << 1;
- else {
- hx = hz << 1;
+ if (ex < ey) {
+ if (ex+1 == ey)
+ goto end;
+ return x;
+ }
+
+ /* x mod y */
+ for (; ex > ey; ex--) {
+ i = uxi - uy.i;
+ if (i >> 31 == 0) {
+ uxi = i;
q++;
}
+ uxi <<= 1;
q <<= 1;
}
- hz = hx - hy;
- if (hz >= 0) {
- hx = hz;
+ i = uxi - uy.i;
+ if (i >> 31 == 0) {
+ uxi = i;
q++;
}
-
- /* convert back to floating value and restore the sign */
- if (hx == 0) { /* return sign(x)*0 */
- q &= 0x7fffffff;
- *quo = sxy ? -q : q;
- return Zero[(uint32_t)sx>>31];
- }
- while (hx < 0x00800000) { /* normalize x */
- hx <<= 1;
- iy--;
+ if (uxi == 0)
+ ex = -30;
+ else
+ for (; uxi>>23 == 0; uxi <<= 1, ex--);
+end:
+ /* scale result and decide between |x| and |x|-|y| */
+ if (ex > 0) {
+ uxi -= 1U << 23;
+ uxi |= (uint32_t)ex << 23;
+ } else {
+ uxi >>= -ex + 1;
}
- if (iy >= -126) { /* normalize output */
- hx = (hx-0x00800000)|((iy+127)<<23);
- } else { /* subnormal output */
- n = -126 - iy;
- hx >>= n;
- }
-fixup:
- SET_FLOAT_WORD(x,hx);
- y = fabsf(y);
- if (y < 0x1p-125f) {
- if (x + x > y || (x + x == y && (q & 1))) {
- q++;
- x -= y;
- }
- } else if (x > 0.5f*y || (x == 0.5f*y && (q & 1))) {
- q++;
+ ux.i = uxi;
+ x = ux.f;
+ if (sy)
+ y = -y;
+ if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
x -= y;
+ q++;
}
- GET_FLOAT_WORD(hx, x);
- SET_FLOAT_WORD(x, hx ^ sx);
q &= 0x7fffffff;
- *quo = sxy ? -q : q;
- return x;
+ *quo = sx^sy ? -(int)q : (int)q;
+ return sx ? -x : x;
}
-/* 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
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
projects / musl / commitdiff