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