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