math: fix pow signed shift ub
[oweals/musl.git] / src / math / fmod.c
1 #include <math.h>
2 #include <stdint.h>
3
4 double 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 }