1#include "libm.h" 2 3#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 4long double remquol(long double x, long double y, int* quo) { 5 return remquo(x, y, quo); 6} 7#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 8long double remquol(long double x, long double y, int* quo) { 9 union ldshape ux = {x}, uy = {y}; 10 int ex = ux.i.se & 0x7fff; 11 int ey = uy.i.se & 0x7fff; 12 int sx = ux.i.se >> 15; 13 int sy = uy.i.se >> 15; 14 uint32_t q; 15 16 *quo = 0; 17 if (y == 0 || isnan(y) || ex == 0x7fff) 18 return (x * y) / (x * y); 19 if (x == 0) 20 return x; 21 22 /* normalize x and y */ 23 if (!ex) { 24 ux.i.se = ex; 25 ux.f *= 0x1p120f; 26 ex = ux.i.se - 120; 27 } 28 if (!ey) { 29 uy.i.se = ey; 30 uy.f *= 0x1p120f; 31 ey = uy.i.se - 120; 32 } 33 34 q = 0; 35 if (ex >= ey) { 36/* x mod y */ 37#if LDBL_MANT_DIG == 64 38 uint64_t i, zx, my; 39 zx = ux.i.m; 40 my = uy.i.m; 41 for (; ex > ey; ex--) { 42 i = zx - my; 43 if (zx >= my) { 44 zx = 2 * i; 45 q++; 46 q <<= 1; 47 } else if (2 * zx < zx) { 48 zx = 2 * zx - my; 49 q <<= 1; 50 q++; 51 } else { 52 zx = 2 * zx; 53 q <<= 1; 54 } 55 } 56 i = zx - my; 57 if (zx >= my) { 58 zx = i; 59 q++; 60 } 61 if (zx == 0) 62 ex = -120; 63 else 64 for (; zx >> 63 == 0; zx *= 2, ex--) 65 ; 66 ux.i.m = zx; 67#elif LDBL_MANT_DIG == 113 68 uint64_t hi, lo, xhi, xlo, yhi, ylo; 69 xhi = (ux.i2.hi & -1ULL >> 16) | 1ULL << 48; 70 yhi = (uy.i2.hi & -1ULL >> 16) | 1ULL << 48; 71 xlo = ux.i2.lo; 72 ylo = ux.i2.lo; 73 for (; ex > ey; ex--) { 74 hi = xhi - yhi; 75 lo = xlo - ylo; 76 if (xlo < ylo) 77 hi -= 1; 78 if (hi >> 63 == 0) { 79 xhi = 2 * hi + (lo >> 63); 80 xlo = 2 * lo; 81 q++; 82 } else { 83 xhi = 2 * xhi + (xlo >> 63); 84 xlo = 2 * xlo; 85 } 86 q <<= 1; 87 } 88 hi = xhi - yhi; 89 lo = xlo - ylo; 90 if (xlo < ylo) 91 hi -= 1; 92 if (hi >> 63 == 0) { 93 xhi = hi; 94 xlo = lo; 95 q++; 96 } 97 if ((xhi | xlo) == 0) 98 ex = -120; 99 else 100 for (; xhi >> 48 == 0; xhi = 2 * xhi + (xlo >> 63), xlo = 2 * xlo, ex--) 101 ; 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