1 /*****************************************************************************/
2 /*****************************************************************************/
3 // fmodf from musl-0.9.15
4 /*****************************************************************************/
5 /*****************************************************************************/
6 
7 #include "libm.h"
8 
fmodf(float x,float y)9 float fmodf(float x, float y)
10 {
11 	union {float f; uint32_t i;} ux = {x}, uy = {y};
12 	int ex = ux.i>>23 & 0xff;
13 	int ey = uy.i>>23 & 0xff;
14 	uint32_t sx = ux.i & 0x80000000;
15 	uint32_t i;
16 	uint32_t uxi = ux.i;
17 
18 	if (uy.i<<1 == 0 || isnan(y) || ex == 0xff)
19 		return (x*y)/(x*y);
20 	if (uxi<<1 <= uy.i<<1) {
21 		if (uxi<<1 == uy.i<<1)
22 			return 0*x;
23 		return x;
24 	}
25 
26 	/* normalize x and y */
27 	if (!ex) {
28 		for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1);
29 		uxi <<= -ex + 1;
30 	} else {
31 		uxi &= -1U >> 9;
32 		uxi |= 1U << 23;
33 	}
34 	if (!ey) {
35 		for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1);
36 		uy.i <<= -ey + 1;
37 	} else {
38 		uy.i &= -1U >> 9;
39 		uy.i |= 1U << 23;
40 	}
41 
42 	/* x mod y */
43 	for (; ex > ey; ex--) {
44 		i = uxi - uy.i;
45 		if (i >> 31 == 0) {
46 			if (i == 0)
47 				return 0*x;
48 			uxi = i;
49 		}
50 		uxi <<= 1;
51 	}
52 	i = uxi - uy.i;
53 	if (i >> 31 == 0) {
54 		if (i == 0)
55 			return 0*x;
56 		uxi = i;
57 	}
58 	for (; uxi>>23 == 0; uxi <<= 1, ex--);
59 
60 	/* scale result up */
61 	if (ex > 0) {
62 		uxi -= 1U << 23;
63 		uxi |= (uint32_t)ex << 23;
64 	} else {
65 		uxi >>= -ex + 1;
66 	}
67 	uxi |= sx;
68 	ux.i = uxi;
69 	return ux.f;
70 }
71