summaryrefslogtreecommitdiff
path: root/src/math/remainderf.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/remainderf.c
parentd1a2ead878c27ac4ec600740320f8b76e1f961e9 (diff)
downloadmusl-ee2ee92d62c43f6658d37ddea4c316d2089d0fe9.tar.gz
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/remainderf.c')
-rw-r--r--src/math/remainderf.c60
1 files changed, 4 insertions, 56 deletions
diff --git a/src/math/remainderf.c b/src/math/remainderf.c
index 77671039..b418bbff 100644
--- a/src/math/remainderf.c
+++ b/src/math/remainderf.c
@@ -1,59 +1,7 @@
-/* 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);
}