1 |
272 |
jeremybenn |
/* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
2 |
|
|
|
3 |
|
|
This file is part of GCC.
|
4 |
|
|
|
5 |
|
|
GCC is free software; you can redistribute it and/or modify it under
|
6 |
|
|
the terms of the GNU General Public License as published by the Free
|
7 |
|
|
Software Foundation; either version 3, or (at your option) any later
|
8 |
|
|
version.
|
9 |
|
|
|
10 |
|
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
11 |
|
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
12 |
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
13 |
|
|
for more details.
|
14 |
|
|
|
15 |
|
|
Under Section 7 of GPL version 3, you are granted additional
|
16 |
|
|
permissions described in the GCC Runtime Library Exception, version
|
17 |
|
|
3.1, as published by the Free Software Foundation.
|
18 |
|
|
|
19 |
|
|
You should have received a copy of the GNU General Public License and
|
20 |
|
|
a copy of the GCC Runtime Library Exception along with this program;
|
21 |
|
|
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
22 |
|
|
<http://www.gnu.org/licenses/>. */
|
23 |
|
|
|
24 |
|
|
#define BID_128RES
|
25 |
|
|
#include "bid_internal.h"
|
26 |
|
|
|
27 |
|
|
/*****************************************************************************
|
28 |
|
|
* BID128 nextup
|
29 |
|
|
****************************************************************************/
|
30 |
|
|
|
31 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
32 |
|
|
void
|
33 |
|
|
bid128_nextup (UINT128 * pres,
|
34 |
|
|
UINT128 *
|
35 |
|
|
px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
36 |
|
|
UINT128 x = *px;
|
37 |
|
|
#else
|
38 |
|
|
UINT128
|
39 |
|
|
bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
40 |
|
|
_EXC_INFO_PARAM) {
|
41 |
|
|
#endif
|
42 |
|
|
|
43 |
|
|
UINT128 res;
|
44 |
|
|
UINT64 x_sign;
|
45 |
|
|
UINT64 x_exp;
|
46 |
|
|
int exp;
|
47 |
|
|
BID_UI64DOUBLE tmp1;
|
48 |
|
|
int x_nr_bits;
|
49 |
|
|
int q1, ind;
|
50 |
|
|
UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
|
51 |
|
|
|
52 |
|
|
BID_SWAP128 (x);
|
53 |
|
|
// unpack the argument
|
54 |
|
|
x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
55 |
|
|
C1.w[1] = x.w[1] & MASK_COEFF;
|
56 |
|
|
C1.w[0] = x.w[0];
|
57 |
|
|
|
58 |
|
|
// check for NaN or Infinity
|
59 |
|
|
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
|
60 |
|
|
// x is special
|
61 |
|
|
if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
62 |
|
|
// if x = NaN, then res = Q (x)
|
63 |
|
|
// check first for non-canonical NaN payload
|
64 |
|
|
if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
65 |
|
|
(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
66 |
|
|
&& (x.w[0] > 0x38c15b09ffffffffull))) {
|
67 |
|
|
x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
68 |
|
|
x.w[0] = 0x0ull;
|
69 |
|
|
}
|
70 |
|
|
if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
71 |
|
|
// set invalid flag
|
72 |
|
|
*pfpsf |= INVALID_EXCEPTION;
|
73 |
|
|
// return quiet (x)
|
74 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
75 |
|
|
res.w[0] = x.w[0];
|
76 |
|
|
} else { // x is QNaN
|
77 |
|
|
// return x
|
78 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
79 |
|
|
res.w[0] = x.w[0];
|
80 |
|
|
}
|
81 |
|
|
} else { // x is not NaN, so it must be infinity
|
82 |
|
|
if (!x_sign) { // x is +inf
|
83 |
|
|
res.w[1] = 0x7800000000000000ull; // +inf
|
84 |
|
|
res.w[0] = 0x0000000000000000ull;
|
85 |
|
|
} else { // x is -inf
|
86 |
|
|
res.w[1] = 0xdfffed09bead87c0ull; // -MAXFP = -999...99 * 10^emax
|
87 |
|
|
res.w[0] = 0x378d8e63ffffffffull;
|
88 |
|
|
}
|
89 |
|
|
}
|
90 |
|
|
BID_RETURN (res);
|
91 |
|
|
}
|
92 |
|
|
// check for non-canonical values (treated as zero)
|
93 |
|
|
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
94 |
|
|
// non-canonical
|
95 |
|
|
x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
96 |
|
|
C1.w[1] = 0; // significand high
|
97 |
|
|
C1.w[0] = 0; // significand low
|
98 |
|
|
} else { // G0_G1 != 11
|
99 |
|
|
x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
100 |
|
|
if (C1.w[1] > 0x0001ed09bead87c0ull ||
|
101 |
|
|
(C1.w[1] == 0x0001ed09bead87c0ull
|
102 |
|
|
&& C1.w[0] > 0x378d8e63ffffffffull)) {
|
103 |
|
|
// x is non-canonical if coefficient is larger than 10^34 -1
|
104 |
|
|
C1.w[1] = 0;
|
105 |
|
|
C1.w[0] = 0;
|
106 |
|
|
} else { // canonical
|
107 |
|
|
;
|
108 |
|
|
}
|
109 |
|
|
}
|
110 |
|
|
|
111 |
|
|
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
|
112 |
|
|
// x is +/-0
|
113 |
|
|
res.w[1] = 0x0000000000000000ull; // +1 * 10^emin
|
114 |
|
|
res.w[0] = 0x0000000000000001ull;
|
115 |
|
|
} else { // x is not special and is not zero
|
116 |
|
|
if (x.w[1] == 0x5fffed09bead87c0ull
|
117 |
|
|
&& x.w[0] == 0x378d8e63ffffffffull) {
|
118 |
|
|
// x = +MAXFP = 999...99 * 10^emax
|
119 |
|
|
res.w[1] = 0x7800000000000000ull; // +inf
|
120 |
|
|
res.w[0] = 0x0000000000000000ull;
|
121 |
|
|
} else if (x.w[1] == 0x8000000000000000ull
|
122 |
|
|
&& x.w[0] == 0x0000000000000001ull) {
|
123 |
|
|
// x = -MINFP = 1...99 * 10^emin
|
124 |
|
|
res.w[1] = 0x8000000000000000ull; // -0
|
125 |
|
|
res.w[0] = 0x0000000000000000ull;
|
126 |
|
|
} else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
127 |
|
|
// can add/subtract 1 ulp to the significand
|
128 |
|
|
|
129 |
|
|
// Note: we could check here if x >= 10^34 to speed up the case q1 = 34
|
130 |
|
|
// q1 = nr. of decimal digits in x
|
131 |
|
|
// determine first the nr. of bits in x
|
132 |
|
|
if (C1.w[1] == 0) {
|
133 |
|
|
if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
|
134 |
|
|
// split the 64-bit value in two 32-bit halves to avoid rnd errors
|
135 |
|
|
if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
|
136 |
|
|
tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
|
137 |
|
|
x_nr_bits =
|
138 |
|
|
33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
139 |
|
|
0x3ff);
|
140 |
|
|
} else { // x < 2^32
|
141 |
|
|
tmp1.d = (double) (C1.w[0]); // exact conversion
|
142 |
|
|
x_nr_bits =
|
143 |
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
144 |
|
|
0x3ff);
|
145 |
|
|
}
|
146 |
|
|
} else { // if x < 2^53
|
147 |
|
|
tmp1.d = (double) C1.w[0]; // exact conversion
|
148 |
|
|
x_nr_bits =
|
149 |
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
150 |
|
|
}
|
151 |
|
|
} else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
|
152 |
|
|
tmp1.d = (double) C1.w[1]; // exact conversion
|
153 |
|
|
x_nr_bits =
|
154 |
|
|
65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
155 |
|
|
}
|
156 |
|
|
q1 = nr_digits[x_nr_bits - 1].digits;
|
157 |
|
|
if (q1 == 0) {
|
158 |
|
|
q1 = nr_digits[x_nr_bits - 1].digits1;
|
159 |
|
|
if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
|
160 |
|
|
|| (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
|
161 |
|
|
&& C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
|
162 |
|
|
q1++;
|
163 |
|
|
}
|
164 |
|
|
// if q1 < P34 then pad the significand with zeros
|
165 |
|
|
if (q1 < P34) {
|
166 |
|
|
exp = (x_exp >> 49) - 6176;
|
167 |
|
|
if (exp + 6176 > P34 - q1) {
|
168 |
|
|
ind = P34 - q1; // 1 <= ind <= P34 - 1
|
169 |
|
|
// pad with P34 - q1 zeros, until exponent = emin
|
170 |
|
|
// C1 = C1 * 10^ind
|
171 |
|
|
if (q1 <= 19) { // 64-bit C1
|
172 |
|
|
if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
|
173 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
174 |
|
|
} else { // 128-bit 10^ind and 64-bit C1
|
175 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
176 |
|
|
}
|
177 |
|
|
} else { // C1 is (most likely) 128-bit
|
178 |
|
|
if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
|
179 |
|
|
__mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
180 |
|
|
} else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
|
181 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
182 |
|
|
} else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
|
183 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
184 |
|
|
}
|
185 |
|
|
}
|
186 |
|
|
x_exp = x_exp - ((UINT64) ind << 49);
|
187 |
|
|
} else { // pad with zeros until the exponent reaches emin
|
188 |
|
|
ind = exp + 6176;
|
189 |
|
|
// C1 = C1 * 10^ind
|
190 |
|
|
if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
|
191 |
|
|
if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
|
192 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
193 |
|
|
} else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
|
194 |
|
|
__mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
195 |
|
|
}
|
196 |
|
|
} else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
|
197 |
|
|
// 64-bit C1, 128-bit 10^ind
|
198 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
199 |
|
|
}
|
200 |
|
|
x_exp = EXP_MIN;
|
201 |
|
|
}
|
202 |
|
|
}
|
203 |
|
|
if (!x_sign) { // x > 0
|
204 |
|
|
// add 1 ulp (add 1 to the significand)
|
205 |
|
|
C1.w[0]++;
|
206 |
|
|
if (C1.w[0] == 0)
|
207 |
|
|
C1.w[1]++;
|
208 |
|
|
if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
|
209 |
|
|
C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
|
210 |
|
|
C1.w[0] = 0x38c15b0a00000000ull;
|
211 |
|
|
x_exp = x_exp + EXP_P1;
|
212 |
|
|
}
|
213 |
|
|
} else { // x < 0
|
214 |
|
|
// subtract 1 ulp (subtract 1 from the significand)
|
215 |
|
|
C1.w[0]--;
|
216 |
|
|
if (C1.w[0] == 0xffffffffffffffffull)
|
217 |
|
|
C1.w[1]--;
|
218 |
|
|
if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
|
219 |
|
|
C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
|
220 |
|
|
C1.w[0] = 0x378d8e63ffffffffull;
|
221 |
|
|
x_exp = x_exp - EXP_P1;
|
222 |
|
|
}
|
223 |
|
|
}
|
224 |
|
|
// assemble the result
|
225 |
|
|
res.w[1] = x_sign | x_exp | C1.w[1];
|
226 |
|
|
res.w[0] = C1.w[0];
|
227 |
|
|
} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
228 |
|
|
} // end x is not special and is not zero
|
229 |
|
|
BID_RETURN (res);
|
230 |
|
|
}
|
231 |
|
|
|
232 |
|
|
/*****************************************************************************
|
233 |
|
|
* BID128 nextdown
|
234 |
|
|
****************************************************************************/
|
235 |
|
|
|
236 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
237 |
|
|
void
|
238 |
|
|
bid128_nextdown (UINT128 * pres,
|
239 |
|
|
UINT128 *
|
240 |
|
|
px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
|
241 |
|
|
UINT128 x = *px;
|
242 |
|
|
#else
|
243 |
|
|
UINT128
|
244 |
|
|
bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
245 |
|
|
_EXC_INFO_PARAM) {
|
246 |
|
|
#endif
|
247 |
|
|
|
248 |
|
|
UINT128 res;
|
249 |
|
|
UINT64 x_sign;
|
250 |
|
|
UINT64 x_exp;
|
251 |
|
|
int exp;
|
252 |
|
|
BID_UI64DOUBLE tmp1;
|
253 |
|
|
int x_nr_bits;
|
254 |
|
|
int q1, ind;
|
255 |
|
|
UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
|
256 |
|
|
|
257 |
|
|
BID_SWAP128 (x);
|
258 |
|
|
// unpack the argument
|
259 |
|
|
x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
|
260 |
|
|
C1.w[1] = x.w[1] & MASK_COEFF;
|
261 |
|
|
C1.w[0] = x.w[0];
|
262 |
|
|
|
263 |
|
|
// check for NaN or Infinity
|
264 |
|
|
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
|
265 |
|
|
// x is special
|
266 |
|
|
if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
267 |
|
|
// if x = NaN, then res = Q (x)
|
268 |
|
|
// check first for non-canonical NaN payload
|
269 |
|
|
if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
270 |
|
|
(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
271 |
|
|
&& (x.w[0] > 0x38c15b09ffffffffull))) {
|
272 |
|
|
x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
273 |
|
|
x.w[0] = 0x0ull;
|
274 |
|
|
}
|
275 |
|
|
if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
276 |
|
|
// set invalid flag
|
277 |
|
|
*pfpsf |= INVALID_EXCEPTION;
|
278 |
|
|
// return quiet (x)
|
279 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
280 |
|
|
res.w[0] = x.w[0];
|
281 |
|
|
} else { // x is QNaN
|
282 |
|
|
// return x
|
283 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
284 |
|
|
res.w[0] = x.w[0];
|
285 |
|
|
}
|
286 |
|
|
} else { // x is not NaN, so it must be infinity
|
287 |
|
|
if (!x_sign) { // x is +inf
|
288 |
|
|
res.w[1] = 0x5fffed09bead87c0ull; // +MAXFP = +999...99 * 10^emax
|
289 |
|
|
res.w[0] = 0x378d8e63ffffffffull;
|
290 |
|
|
} else { // x is -inf
|
291 |
|
|
res.w[1] = 0xf800000000000000ull; // -inf
|
292 |
|
|
res.w[0] = 0x0000000000000000ull;
|
293 |
|
|
}
|
294 |
|
|
}
|
295 |
|
|
BID_RETURN (res);
|
296 |
|
|
}
|
297 |
|
|
// check for non-canonical values (treated as zero)
|
298 |
|
|
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
299 |
|
|
// non-canonical
|
300 |
|
|
x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
301 |
|
|
C1.w[1] = 0; // significand high
|
302 |
|
|
C1.w[0] = 0; // significand low
|
303 |
|
|
} else { // G0_G1 != 11
|
304 |
|
|
x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
305 |
|
|
if (C1.w[1] > 0x0001ed09bead87c0ull ||
|
306 |
|
|
(C1.w[1] == 0x0001ed09bead87c0ull
|
307 |
|
|
&& C1.w[0] > 0x378d8e63ffffffffull)) {
|
308 |
|
|
// x is non-canonical if coefficient is larger than 10^34 -1
|
309 |
|
|
C1.w[1] = 0;
|
310 |
|
|
C1.w[0] = 0;
|
311 |
|
|
} else { // canonical
|
312 |
|
|
;
|
313 |
|
|
}
|
314 |
|
|
}
|
315 |
|
|
|
316 |
|
|
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
|
317 |
|
|
// x is +/-0
|
318 |
|
|
res.w[1] = 0x8000000000000000ull; // -1 * 10^emin
|
319 |
|
|
res.w[0] = 0x0000000000000001ull;
|
320 |
|
|
} else { // x is not special and is not zero
|
321 |
|
|
if (x.w[1] == 0xdfffed09bead87c0ull
|
322 |
|
|
&& x.w[0] == 0x378d8e63ffffffffull) {
|
323 |
|
|
// x = -MAXFP = -999...99 * 10^emax
|
324 |
|
|
res.w[1] = 0xf800000000000000ull; // -inf
|
325 |
|
|
res.w[0] = 0x0000000000000000ull;
|
326 |
|
|
} else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) { // +MINFP
|
327 |
|
|
res.w[1] = 0x0000000000000000ull; // +0
|
328 |
|
|
res.w[0] = 0x0000000000000000ull;
|
329 |
|
|
} else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
330 |
|
|
// can add/subtract 1 ulp to the significand
|
331 |
|
|
|
332 |
|
|
// Note: we could check here if x >= 10^34 to speed up the case q1 = 34
|
333 |
|
|
// q1 = nr. of decimal digits in x
|
334 |
|
|
// determine first the nr. of bits in x
|
335 |
|
|
if (C1.w[1] == 0) {
|
336 |
|
|
if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
|
337 |
|
|
// split the 64-bit value in two 32-bit halves to avoid rnd errors
|
338 |
|
|
if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
|
339 |
|
|
tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
|
340 |
|
|
x_nr_bits =
|
341 |
|
|
33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
342 |
|
|
0x3ff);
|
343 |
|
|
} else { // x < 2^32
|
344 |
|
|
tmp1.d = (double) (C1.w[0]); // exact conversion
|
345 |
|
|
x_nr_bits =
|
346 |
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
|
347 |
|
|
0x3ff);
|
348 |
|
|
}
|
349 |
|
|
} else { // if x < 2^53
|
350 |
|
|
tmp1.d = (double) C1.w[0]; // exact conversion
|
351 |
|
|
x_nr_bits =
|
352 |
|
|
1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
353 |
|
|
}
|
354 |
|
|
} else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
|
355 |
|
|
tmp1.d = (double) C1.w[1]; // exact conversion
|
356 |
|
|
x_nr_bits =
|
357 |
|
|
65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
|
358 |
|
|
}
|
359 |
|
|
q1 = nr_digits[x_nr_bits - 1].digits;
|
360 |
|
|
if (q1 == 0) {
|
361 |
|
|
q1 = nr_digits[x_nr_bits - 1].digits1;
|
362 |
|
|
if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
|
363 |
|
|
|| (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
|
364 |
|
|
&& C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
|
365 |
|
|
q1++;
|
366 |
|
|
}
|
367 |
|
|
// if q1 < P then pad the significand with zeros
|
368 |
|
|
if (q1 < P34) {
|
369 |
|
|
exp = (x_exp >> 49) - 6176;
|
370 |
|
|
if (exp + 6176 > P34 - q1) {
|
371 |
|
|
ind = P34 - q1; // 1 <= ind <= P34 - 1
|
372 |
|
|
// pad with P34 - q1 zeros, until exponent = emin
|
373 |
|
|
// C1 = C1 * 10^ind
|
374 |
|
|
if (q1 <= 19) { // 64-bit C1
|
375 |
|
|
if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
|
376 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
377 |
|
|
} else { // 128-bit 10^ind and 64-bit C1
|
378 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
379 |
|
|
}
|
380 |
|
|
} else { // C1 is (most likely) 128-bit
|
381 |
|
|
if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
|
382 |
|
|
__mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
383 |
|
|
} else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
|
384 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
385 |
|
|
} else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
|
386 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
387 |
|
|
}
|
388 |
|
|
}
|
389 |
|
|
x_exp = x_exp - ((UINT64) ind << 49);
|
390 |
|
|
} else { // pad with zeros until the exponent reaches emin
|
391 |
|
|
ind = exp + 6176;
|
392 |
|
|
// C1 = C1 * 10^ind
|
393 |
|
|
if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
|
394 |
|
|
if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
|
395 |
|
|
__mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
|
396 |
|
|
} else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
|
397 |
|
|
__mul_128x64_to_128 (C1, ten2k64[ind], C1);
|
398 |
|
|
}
|
399 |
|
|
} else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
|
400 |
|
|
// 64-bit C1, 128-bit 10^ind
|
401 |
|
|
__mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
|
402 |
|
|
}
|
403 |
|
|
x_exp = EXP_MIN;
|
404 |
|
|
}
|
405 |
|
|
}
|
406 |
|
|
if (x_sign) { // x < 0
|
407 |
|
|
// add 1 ulp (add 1 to the significand)
|
408 |
|
|
C1.w[0]++;
|
409 |
|
|
if (C1.w[0] == 0)
|
410 |
|
|
C1.w[1]++;
|
411 |
|
|
if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
|
412 |
|
|
C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
|
413 |
|
|
C1.w[0] = 0x38c15b0a00000000ull;
|
414 |
|
|
x_exp = x_exp + EXP_P1;
|
415 |
|
|
}
|
416 |
|
|
} else { // x > 0
|
417 |
|
|
// subtract 1 ulp (subtract 1 from the significand)
|
418 |
|
|
C1.w[0]--;
|
419 |
|
|
if (C1.w[0] == 0xffffffffffffffffull)
|
420 |
|
|
C1.w[1]--;
|
421 |
|
|
if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
|
422 |
|
|
C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
|
423 |
|
|
C1.w[0] = 0x378d8e63ffffffffull;
|
424 |
|
|
x_exp = x_exp - EXP_P1;
|
425 |
|
|
}
|
426 |
|
|
}
|
427 |
|
|
// assemble the result
|
428 |
|
|
res.w[1] = x_sign | x_exp | C1.w[1];
|
429 |
|
|
res.w[0] = C1.w[0];
|
430 |
|
|
} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
|
431 |
|
|
} // end x is not special and is not zero
|
432 |
|
|
BID_RETURN (res);
|
433 |
|
|
}
|
434 |
|
|
|
435 |
|
|
/*****************************************************************************
|
436 |
|
|
* BID128 nextafter
|
437 |
|
|
****************************************************************************/
|
438 |
|
|
|
439 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
440 |
|
|
void
|
441 |
|
|
bid128_nextafter (UINT128 * pres, UINT128 * px,
|
442 |
|
|
UINT128 *
|
443 |
|
|
py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
|
444 |
|
|
{
|
445 |
|
|
UINT128 x = *px;
|
446 |
|
|
UINT128 y = *py;
|
447 |
|
|
UINT128 xnswp = *px;
|
448 |
|
|
UINT128 ynswp = *py;
|
449 |
|
|
#else
|
450 |
|
|
UINT128
|
451 |
|
|
bid128_nextafter (UINT128 x,
|
452 |
|
|
UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
|
453 |
|
|
_EXC_INFO_PARAM) {
|
454 |
|
|
UINT128 xnswp = x;
|
455 |
|
|
UINT128 ynswp = y;
|
456 |
|
|
#endif
|
457 |
|
|
|
458 |
|
|
UINT128 res;
|
459 |
|
|
UINT128 tmp1, tmp2, tmp3;
|
460 |
|
|
FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
|
461 |
|
|
int res1, res2;
|
462 |
|
|
UINT64 x_exp;
|
463 |
|
|
|
464 |
|
|
|
465 |
|
|
BID_SWAP128 (x);
|
466 |
|
|
BID_SWAP128 (y);
|
467 |
|
|
// check for NaNs
|
468 |
|
|
if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
|
469 |
|
|
|| ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
|
470 |
|
|
// x is special or y is special
|
471 |
|
|
if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
|
472 |
|
|
// if x = NaN, then res = Q (x)
|
473 |
|
|
// check first for non-canonical NaN payload
|
474 |
|
|
if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
475 |
|
|
(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
476 |
|
|
&& (x.w[0] > 0x38c15b09ffffffffull))) {
|
477 |
|
|
x.w[1] = x.w[1] & 0xffffc00000000000ull;
|
478 |
|
|
x.w[0] = 0x0ull;
|
479 |
|
|
}
|
480 |
|
|
if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
|
481 |
|
|
// set invalid flag
|
482 |
|
|
*pfpsf |= INVALID_EXCEPTION;
|
483 |
|
|
// return quiet (x)
|
484 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
485 |
|
|
res.w[0] = x.w[0];
|
486 |
|
|
} else { // x is QNaN
|
487 |
|
|
// return x
|
488 |
|
|
res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
489 |
|
|
res.w[0] = x.w[0];
|
490 |
|
|
if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
491 |
|
|
// set invalid flag
|
492 |
|
|
*pfpsf |= INVALID_EXCEPTION;
|
493 |
|
|
}
|
494 |
|
|
}
|
495 |
|
|
BID_RETURN (res)
|
496 |
|
|
} else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
|
497 |
|
|
// if x = NaN, then res = Q (x)
|
498 |
|
|
// check first for non-canonical NaN payload
|
499 |
|
|
if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
|
500 |
|
|
(((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
|
501 |
|
|
&& (y.w[0] > 0x38c15b09ffffffffull))) {
|
502 |
|
|
y.w[1] = y.w[1] & 0xffffc00000000000ull;
|
503 |
|
|
y.w[0] = 0x0ull;
|
504 |
|
|
}
|
505 |
|
|
if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
|
506 |
|
|
// set invalid flag
|
507 |
|
|
*pfpsf |= INVALID_EXCEPTION;
|
508 |
|
|
// return quiet (x)
|
509 |
|
|
res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
|
510 |
|
|
res.w[0] = y.w[0];
|
511 |
|
|
} else { // x is QNaN
|
512 |
|
|
// return x
|
513 |
|
|
res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
|
514 |
|
|
res.w[0] = y.w[0];
|
515 |
|
|
}
|
516 |
|
|
BID_RETURN (res)
|
517 |
|
|
} else { // at least one is infinity
|
518 |
|
|
if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
|
519 |
|
|
x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
|
520 |
|
|
x.w[0] = 0x0ull;
|
521 |
|
|
}
|
522 |
|
|
if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
|
523 |
|
|
y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
|
524 |
|
|
y.w[0] = 0x0ull;
|
525 |
|
|
}
|
526 |
|
|
}
|
527 |
|
|
}
|
528 |
|
|
// neither x nor y is NaN
|
529 |
|
|
|
530 |
|
|
// if not infinity, check for non-canonical values x (treated as zero)
|
531 |
|
|
if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
|
532 |
|
|
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
|
533 |
|
|
// non-canonical
|
534 |
|
|
x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
|
535 |
|
|
x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
|
536 |
|
|
x.w[0] = 0x0ull;
|
537 |
|
|
} else { // G0_G1 != 11
|
538 |
|
|
x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
|
539 |
|
|
if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
|
540 |
|
|
((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
|
541 |
|
|
&& x.w[0] > 0x378d8e63ffffffffull)) {
|
542 |
|
|
// x is non-canonical if coefficient is larger than 10^34 -1
|
543 |
|
|
x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
|
544 |
|
|
x.w[0] = 0x0ull;
|
545 |
|
|
} else { // canonical
|
546 |
|
|
;
|
547 |
|
|
}
|
548 |
|
|
}
|
549 |
|
|
}
|
550 |
|
|
// no need to check for non-canonical y
|
551 |
|
|
|
552 |
|
|
// neither x nor y is NaN
|
553 |
|
|
tmp_fpsf = *pfpsf; // save fpsf
|
554 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
555 |
|
|
bid128_quiet_equal (&res1, &xnswp,
|
556 |
|
|
&ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
557 |
|
|
_EXC_INFO_ARG);
|
558 |
|
|
bid128_quiet_greater (&res2, &xnswp,
|
559 |
|
|
&ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
560 |
|
|
_EXC_INFO_ARG);
|
561 |
|
|
#else
|
562 |
|
|
res1 =
|
563 |
|
|
bid128_quiet_equal (xnswp,
|
564 |
|
|
ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
565 |
|
|
_EXC_INFO_ARG);
|
566 |
|
|
res2 =
|
567 |
|
|
bid128_quiet_greater (xnswp,
|
568 |
|
|
ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
569 |
|
|
_EXC_INFO_ARG);
|
570 |
|
|
#endif
|
571 |
|
|
*pfpsf = tmp_fpsf; // restore fpsf
|
572 |
|
|
|
573 |
|
|
if (res1) { // x = y
|
574 |
|
|
// return x with the sign of y
|
575 |
|
|
res.w[1] =
|
576 |
|
|
(x.w[1] & 0x7fffffffffffffffull) | (y.
|
577 |
|
|
w[1] & 0x8000000000000000ull);
|
578 |
|
|
res.w[0] = x.w[0];
|
579 |
|
|
} else if (res2) { // x > y
|
580 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
581 |
|
|
bid128_nextdown (&res,
|
582 |
|
|
&xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
583 |
|
|
_EXC_INFO_ARG);
|
584 |
|
|
#else
|
585 |
|
|
res =
|
586 |
|
|
bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
587 |
|
|
_EXC_INFO_ARG);
|
588 |
|
|
#endif
|
589 |
|
|
BID_SWAP128 (res);
|
590 |
|
|
} else { // x < y
|
591 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
592 |
|
|
bid128_nextup (&res,
|
593 |
|
|
&xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
594 |
|
|
#else
|
595 |
|
|
res =
|
596 |
|
|
bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
|
597 |
|
|
#endif
|
598 |
|
|
BID_SWAP128 (res);
|
599 |
|
|
}
|
600 |
|
|
// if the operand x is finite but the result is infinite, signal
|
601 |
|
|
// overflow and inexact
|
602 |
|
|
if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
|
603 |
|
|
&& ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
|
604 |
|
|
// set the inexact flag
|
605 |
|
|
*pfpsf |= INEXACT_EXCEPTION;
|
606 |
|
|
// set the overflow flag
|
607 |
|
|
*pfpsf |= OVERFLOW_EXCEPTION;
|
608 |
|
|
}
|
609 |
|
|
// if the result is in (-10^emin, 10^emin), and is different from the
|
610 |
|
|
// operand x, signal underflow and inexact
|
611 |
|
|
tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull;
|
612 |
|
|
tmp1.w[LOW_128W] = 0x38c15b0a00000000ull; // +100...0[34] * 10^emin
|
613 |
|
|
tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
|
614 |
|
|
tmp2.w[LOW_128W] = res.w[0];
|
615 |
|
|
tmp3.w[HIGH_128W] = res.w[1];
|
616 |
|
|
tmp3.w[LOW_128W] = res.w[0];
|
617 |
|
|
tmp_fpsf = *pfpsf; // save fpsf
|
618 |
|
|
#if DECIMAL_CALL_BY_REFERENCE
|
619 |
|
|
bid128_quiet_greater (&res1, &tmp1,
|
620 |
|
|
&tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
621 |
|
|
_EXC_INFO_ARG);
|
622 |
|
|
bid128_quiet_not_equal (&res2, &xnswp,
|
623 |
|
|
&tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
624 |
|
|
_EXC_INFO_ARG);
|
625 |
|
|
#else
|
626 |
|
|
res1 =
|
627 |
|
|
bid128_quiet_greater (tmp1,
|
628 |
|
|
tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
629 |
|
|
_EXC_INFO_ARG);
|
630 |
|
|
res2 =
|
631 |
|
|
bid128_quiet_not_equal (xnswp,
|
632 |
|
|
tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
|
633 |
|
|
_EXC_INFO_ARG);
|
634 |
|
|
#endif
|
635 |
|
|
*pfpsf = tmp_fpsf; // restore fpsf
|
636 |
|
|
if (res1 && res2) {
|
637 |
|
|
// set the inexact flag
|
638 |
|
|
*pfpsf |= INEXACT_EXCEPTION;
|
639 |
|
|
// set the underflow flag
|
640 |
|
|
*pfpsf |= UNDERFLOW_EXCEPTION;
|
641 |
|
|
}
|
642 |
|
|
BID_RETURN (res);
|
643 |
|
|
}
|