1 #include "libm.h"
2 
3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
fmodl(long double x,long double y)4 long double fmodl(long double x, long double y) {
5     return fmod(x, y);
6 }
7 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
fmodl(long double x,long double y)8 long double fmodl(long double x, long double y) {
9     union ldshape ux = {x}, uy = {y};
10     int ex = ux.i.se & 0x7fff;
11     int ey = uy.i.se & 0x7fff;
12     int sx = ux.i.se & 0x8000;
13 
14     if (y == 0 || isnan(y) || ex == 0x7fff)
15         return (x * y) / (x * y);
16     ux.i.se = ex;
17     uy.i.se = ey;
18     if (ux.f <= uy.f) {
19         if (ux.f == uy.f)
20             return 0 * x;
21         return x;
22     }
23 
24     /* normalize x and y */
25     if (!ex) {
26         ux.f *= 0x1p120f;
27         ex = ux.i.se - 120;
28     }
29     if (!ey) {
30         uy.f *= 0x1p120f;
31         ey = uy.i.se - 120;
32     }
33 
34 /* x mod y */
35 #if LDBL_MANT_DIG == 64
36     uint64_t i, zx, my;
37     zx = ux.i.m;
38     my = uy.i.m;
39     for (; ex > ey; ex--) {
40         i = zx - my;
41         if (zx >= my) {
42             if (i == 0)
43                 return 0 * x;
44             zx = 2 * i;
45         } else if (2 * zx < zx) {
46             zx = 2 * zx - my;
47         } else {
48             zx = 2 * zx;
49         }
50     }
51     i = zx - my;
52     if (zx >= my) {
53         if (i == 0)
54             return 0 * x;
55         zx = i;
56     }
57     for (; zx >> 63 == 0; zx *= 2, ex--)
58         ;
59     ux.i.m = zx;
60 #elif LDBL_MANT_DIG == 113
61     uint64_t hi, lo, xhi, xlo, yhi, ylo;
62     xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48;
63     yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48;
64     xlo = ux.i2.lo;
65     ylo = uy.i2.lo;
66     for (; ex > ey; ex--) {
67         hi = xhi - yhi;
68         lo = xlo - ylo;
69         if (xlo < ylo)
70             hi -= 1;
71         if (hi >> 63 == 0) {
72             if ((hi | lo) == 0)
73                 return 0 * x;
74             xhi = 2 * hi + (lo >> 63);
75             xlo = 2 * lo;
76         } else {
77             xhi = 2 * xhi + (xlo >> 63);
78             xlo = 2 * xlo;
79         }
80     }
81     hi = xhi - yhi;
82     lo = xlo - ylo;
83     if (xlo < ylo)
84         hi -= 1;
85     if (hi >> 63 == 0) {
86         if ((hi | lo) == 0)
87             return 0 * x;
88         xhi = hi;
89         xlo = lo;
90     }
91     for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--)
92         ;
93     ux.i2.hi = xhi;
94     ux.i2.lo = xlo;
95 #endif
96 
97     /* scale result */
98     if (ex <= 0) {
99         ux.i.se = (ex + 120) | sx;
100         ux.f *= 0x1p-120f;
101     } else
102         ux.i.se = ex | sx;
103     return ux.f;
104 }
105 #endif
106