summaryrefslogtreecommitdiff
path: root/src/math
diff options
context:
space:
mode:
authornsz <nsz@port70.net>2012-03-29 14:03:18 +0200
committernsz <nsz@port70.net>2012-03-29 14:03:18 +0200
commitf6ceccd92247575e4a35bc94f581a570b8052d43 (patch)
tree817b7531ea3b62c845f7e3731e48439d4be9745a /src/math
parent9f58d06007818c43b9504e959ef81ab5f113b374 (diff)
downloadmusl-f6ceccd92247575e4a35bc94f581a570b8052d43.tar.gz
musl-f6ceccd92247575e4a35bc94f581a570b8052d43.tar.bz2
musl-f6ceccd92247575e4a35bc94f581a570b8052d43.tar.xz
musl-f6ceccd92247575e4a35bc94f581a570b8052d43.zip
math: rewrite modf.c and clean up modff.c
cleaner implementation with unions and unsigned arithmetic
Diffstat (limited to 'src/math')
-rw-r--r--src/math/modf.c89
-rw-r--r--src/math/modff.c28
2 files changed, 45 insertions, 72 deletions
diff --git a/src/math/modf.c b/src/math/modf.c
index cca3b652..de45069f 100644
--- a/src/math/modf.c
+++ b/src/math/modf.c
@@ -1,68 +1,37 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_modf.c */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-/*
- * modf(double x, double *iptr)
- * return fraction part of x, and return x's integral part in *iptr.
- * Method:
- * Bit twiddling.
- *
- * Exception:
- * No exception.
- */
-
-#include "libm.h"
+#include <math.h>
+#include <stdint.h>
double modf(double x, double *iptr)
{
- int32_t i0,i1,j0;
- uint32_t i;
+ union {double x; uint64_t n;} u = {x};
+ uint64_t mask;
+ int e;
- EXTRACT_WORDS(i0, i1, x);
- j0 = ((i0>>20) & 0x7ff) - 0x3ff; /* exponent of x */
- if (j0 < 20) { /* integer part in high x */
- if (j0 < 0) { /* |x| < 1 */
- INSERT_WORDS(*iptr, i0 & 0x80000000, 0); /* *iptr = +-0 */
- return x;
- }
- i = 0x000fffff >> j0;
- if (((i0&i)|i1) == 0) { /* x is integral */
- uint32_t high;
- *iptr = x;
- GET_HIGH_WORD(high, x);
- INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
- return x;
- }
- INSERT_WORDS(*iptr, i0&~i, 0);
- return x - *iptr;
- } else if (j0 > 51) { /* no fraction part */
- uint32_t high;
- if (j0 == 0x400) { /* inf/NaN */
- *iptr = x;
- return 0.0 / x;
- }
+ e = (int)(u.n>>52 & 0x7ff) - 0x3ff;
+
+ /* no fractional part */
+ if (e >= 52) {
*iptr = x;
- GET_HIGH_WORD(high, x);
- INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
- return x;
- } else { /* fraction part in low x */
- i = (uint32_t)0xffffffff >> (j0 - 20);
- if ((i1&i) == 0) { /* x is integral */
- uint32_t high;
- *iptr = x;
- GET_HIGH_WORD(high, x);
- INSERT_WORDS(x, high & 0x80000000, 0); /* return +-0 */
+ if (e == 0x400 && u.n<<12 != 0) /* nan */
return x;
- }
- INSERT_WORDS(*iptr, i0, i1&~i);
- return x - *iptr;
+ u.n &= (uint64_t)1<<63;
+ return u.x;
+ }
+
+ /* no integral part*/
+ if (e < 0) {
+ u.n &= (uint64_t)1<<63;
+ *iptr = u.x;
+ return x;
+ }
+
+ mask = (uint64_t)-1>>12 >> e;
+ if ((u.n & mask) == 0) {
+ *iptr = x;
+ u.n &= (uint64_t)1<<63;
+ return u.x;
}
+ u.n &= ~mask;
+ *iptr = u.x;
+ return x - *iptr;
}
diff --git a/src/math/modff.c b/src/math/modff.c
index bf6e4ced..84d0b82a 100644
--- a/src/math/modff.c
+++ b/src/math/modff.c
@@ -1,33 +1,37 @@
-#include "libm.h"
+#include <math.h>
+#include <stdint.h>
float modff(float x, float *iptr)
{
- uint32_t u, mask;
+ union {float x; uint32_t n;} u = {x};
+ uint32_t mask;
int e;
- GET_FLOAT_WORD(u, x);
- e = (int)(u>>23 & 0xff) - 0x7f;
+ e = (int)(u.n>>23 & 0xff) - 0x7f;
/* no fractional part */
if (e >= 23) {
*iptr = x;
- if (e == 0x80 && u<<9 != 0) /* nan */
+ if (e == 0x80 && u.n<<9 != 0) { /* nan */
return x;
- SET_FLOAT_WORD(x, u & 0x80000000);
- return x;
+ }
+ u.n &= 0x80000000;
+ return u.x;
}
/* no integral part */
if (e < 0) {
- SET_FLOAT_WORD(*iptr, u & 0x80000000);
+ u.n &= 0x80000000;
+ *iptr = u.x;
return x;
}
mask = 0x007fffff>>e;
- if ((u & mask) == 0) {
+ if ((u.n & mask) == 0) {
*iptr = x;
- SET_FLOAT_WORD(x, u & 0x80000000);
- return x;
+ u.n &= 0x80000000;
+ return u.x;
}
- SET_FLOAT_WORD(*iptr, u & ~mask);
+ u.n &= ~mask;
+ *iptr = u.x;
return x - *iptr;
}