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