aboutsummaryrefslogtreecommitdiff
path: root/lib/libm_dbl/fmod.c
diff options
context:
space:
mode:
authorScott Shawcroft <scott@tannewt.org>2017-08-28 13:59:26 -0700
committerGitHub <noreply@github.com>2017-08-28 13:59:26 -0700
commit447a4b1ecd5655357c9daaace03b91355d2e2c6d (patch)
tree46ca08a93e9870aab4fbb533843f3c99796cbdf4 /lib/libm_dbl/fmod.c
parent4e63d55d8d0df6d7600bb684fa90a8cb54708ba8 (diff)
parentc679c80c710fd4a0e37ad87200d8a721e19727aa (diff)
Merge pull request #208 from adafruit/merge-v1.9.2
Merge MicroPython v1.9.2
Diffstat (limited to 'lib/libm_dbl/fmod.c')
-rw-r--r--lib/libm_dbl/fmod.c68
1 files changed, 68 insertions, 0 deletions
diff --git a/lib/libm_dbl/fmod.c b/lib/libm_dbl/fmod.c
new file mode 100644
index 000000000..6849722ba
--- /dev/null
+++ b/lib/libm_dbl/fmod.c
@@ -0,0 +1,68 @@
+#include <math.h>
+#include <stdint.h>
+
+double fmod(double x, double 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;
+ uint64_t i;
+
+ /* 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;
+
+ if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
+ return (x*y)/(x*y);
+ if (uxi<<1 <= uy.i<<1) {
+ if (uxi<<1 == uy.i<<1)
+ return 0*x;
+ return x;
+ }
+
+ /* 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 (!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;
+ }
+
+ /* 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;
+ }
+ i = uxi - uy.i;
+ if (i >> 63 == 0) {
+ if (i == 0)
+ return 0*x;
+ uxi = i;
+ }
+ for (; uxi>>52 == 0; uxi <<= 1, ex--);
+
+ /* scale result */
+ if (ex > 0) {
+ uxi -= 1ULL << 52;
+ uxi |= (uint64_t)ex << 52;
+ } else {
+ uxi >>= -ex + 1;
+ }
+ uxi |= (uint64_t)sx << 63;
+ ux.i = uxi;
+ return ux.f;
+}