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 |