Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 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 |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 9 | long double remquol(long double x, long double y, int *quo) |
| 10 | { |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 17 | |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 18 | *quo = 0; |
| 19 | if (y == 0 || isnan(y) || ex == 0x7fff) |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 20 | return (x*y)/(x*y); |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 34 | } |
| 35 | |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 36 | q = 0; |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 61 | q++; |
| 62 | } |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 88 | } |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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 |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 105 | } |
| 106 | |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 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 | } |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 120 | q &= 0x7fffffff; |
Szabolcs Nagy | ee2ee92 | 2013-09-03 04:09:12 +0000 | [diff] [blame] | 121 | *quo = sx^sy ? -(int)q : (int)q; |
| 122 | return sx ? -x : x; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 123 | } |
| 124 | #endif |