1#include <math.h>
2#include <stdint.h>
3
4double fmod(double x, double y)
5{
6	union {double f; uint64_t i;} ux = {x}, uy = {y};
7	int ex = ux.i>>52 & 0x7ff;
8	int ey = uy.i>>52 & 0x7ff;
9	int sx = ux.i>>63;
10	uint64_t i;
11
12	/* in the followings uxi should be ux.i, but then gcc wrongly adds */
13	/* float load/store to inner loops ruining performance and code size */
14	uint64_t uxi = ux.i;
15
16	if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
17		return (x*y)/(x*y);
18	if (uxi<<1 <= uy.i<<1) {
19		if (uxi<<1 == uy.i<<1)
20			return 0*x;
21		return x;
22	}
23
24	/* normalize x and y */
25	if (!ex) {
26		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
27		uxi <<= -ex + 1;
28	} else {
29		uxi &= -1ULL >> 12;
30		uxi |= 1ULL << 52;
31	}
32	if (!ey) {
33		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
34		uy.i <<= -ey + 1;
35	} else {
36		uy.i &= -1ULL >> 12;
37		uy.i |= 1ULL << 52;
38	}
39
40	/* x mod y */
41	for (; ex > ey; ex--) {
42		i = uxi - uy.i;
43		if (i >> 63 == 0) {
44			if (i == 0)
45				return 0*x;
46			uxi = i;
47		}
48		uxi <<= 1;
49	}
50	i = uxi - uy.i;
51	if (i >> 63 == 0) {
52		if (i == 0)
53			return 0*x;
54		uxi = i;
55	}
56	for (; uxi>>52 == 0; uxi <<= 1, ex--);
57
58	/* scale result */
59	if (ex > 0) {
60		uxi -= 1ULL << 52;
61		uxi |= (uint64_t)ex << 52;
62	} else {
63		uxi >>= -ex + 1;
64	}
65	uxi |= (uint64_t)sx << 63;
66	ux.i = uxi;
67	return ux.f;
68}
69