1 #include <math.h>
2 #include <stdint.h>
3 
remquo(double x,double y,int * quo)4 double 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--);
63 end:
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