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