1#include <math.h>
2#include <stdint.h>
3
4double remquo(double x, double y, int *quo)
5{
6	union {double f; uint64_t i;} ux = {x}, uy = {y};
7	int ex = ux.i>>52 & 0x7ff;
8	int ey = uy.i>>52 & 0x7ff;
9	int sx = ux.i>>63;
10	int sy = uy.i>>63;
11	uint32_t q;
12	uint64_t i;
13	uint64_t uxi = ux.i;
14
15	*quo = 0;
16	if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff)
17		return (x*y)/(x*y);
18	if (ux.i<<1 == 0)
19		return x;
20
21	/* normalize x and y */
22	if (!ex) {
23		for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1);
24		uxi <<= -ex + 1;
25	} else {
26		uxi &= -1ULL >> 12;
27		uxi |= 1ULL << 52;
28	}
29	if (!ey) {
30		for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1);
31		uy.i <<= -ey + 1;
32	} else {
33		uy.i &= -1ULL >> 12;
34		uy.i |= 1ULL << 52;
35	}
36
37	q = 0;
38	if (ex < ey) {
39		if (ex+1 == ey)
40			goto end;
41		return x;
42	}
43
44	/* x mod y */
45	for (; ex > ey; ex--) {
46		i = uxi - uy.i;
47		if (i >> 63 == 0) {
48			uxi = i;
49			q++;
50		}
51		uxi <<= 1;
52		q <<= 1;
53	}
54	i = uxi - uy.i;
55	if (i >> 63 == 0) {
56		uxi = i;
57		q++;
58	}
59	if (uxi == 0)
60		ex = -60;
61	else
62		for (; uxi>>52 == 0; uxi <<= 1, ex--);
63end:
64	/* scale result and decide between |x| and |x|-|y| */
65	if (ex > 0) {
66		uxi -= 1ULL << 52;
67		uxi |= (uint64_t)ex << 52;
68	} else {
69		uxi >>= -ex + 1;
70	}
71	ux.i = uxi;
72	x = ux.f;
73	if (sy)
74		y = -y;
75	if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
76		x -= y;
77		q++;
78	}
79	q &= 0x7fffffff;
80	*quo = sx^sy ? -(int)q : (int)q;
81	return sx ? -x : x;
82}
83