1 #include "libm.h" 2 3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 4 long 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 9 long 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 125