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