summaryrefslogtreecommitdiff
path: root/src/math/fmodf.c
diff options
context:
space:
mode:
authorSzabolcs Nagy <nsz@port70.net>2013-09-03 04:09:12 +0000
committerSzabolcs Nagy <nsz@port70.net>2013-09-05 11:30:07 +0000
commitee2ee92d62c43f6658d37ddea4c316d2089d0fe9 (patch)
tree1a3e5a63fd40fa763f87d5dd9efd021117ea1d11 /src/math/fmodf.c
parentd1a2ead878c27ac4ec600740320f8b76e1f961e9 (diff)
downloadmusl-ee2ee92d62c43f6658d37ddea4c316d2089d0fe9.tar.gz
musl-ee2ee92d62c43f6658d37ddea4c316d2089d0fe9.tar.bz2
musl-ee2ee92d62c43f6658d37ddea4c316d2089d0fe9.tar.xz
musl-ee2ee92d62c43f6658d37ddea4c316d2089d0fe9.zip
math: rewrite remainder functions (remainder, remquo, fmod, modf)
* results are exact * modfl follows truncl (raises inexact flag spuriously now) * modf and modff only had cosmetic cleanup * remainder is just a wrapper around remquo now * using iterative shift+subtract for remquo and fmod * ld80 and ld128 are supported as well
Diffstat (limited to 'src/math/fmodf.c')
-rw-r--r--src/math/fmodf.c139
1 files changed, 50 insertions, 89 deletions
diff --git a/src/math/fmodf.c b/src/math/fmodf.c
index 837e9371..ff58f933 100644
--- a/src/math/fmodf.c
+++ b/src/math/fmodf.c
@@ -1,104 +1,65 @@
-/* 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;
}