Add a "fork" of musl as //fusl.

fusl/src/math/fmod.c

+#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;
+}