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