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