]>
Commit | Line | Data |
---|---|---|
a945c346 | 1 | /* Copyright (C) 2007-2024 Free Software Foundation, Inc. |
200359e8 L |
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 | |
748086b7 | 7 | Software Foundation; either version 3, or (at your option) any later |
200359e8 L |
8 | version. |
9 | ||
200359e8 L |
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 | ||
748086b7 JJ |
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/>. */ | |
200359e8 L |
23 | |
24 | #include "bid_internal.h" | |
25 | ||
26 | /***************************************************************************** | |
27 | * BID128_to_int32_rnint | |
28 | ****************************************************************************/ | |
29 | ||
b2a00c89 | 30 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rnint, x) |
200359e8 | 31 | |
b2a00c89 L |
32 | int res; |
33 | UINT64 x_sign; | |
34 | UINT64 x_exp; | |
35 | int exp; // unbiased exponent | |
200359e8 | 36 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
37 | UINT64 tmp64; |
38 | BID_UI64DOUBLE tmp1; | |
39 | unsigned int x_nr_bits; | |
40 | int q, ind, shift; | |
41 | UINT128 C1, C; | |
42 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
43 | UINT256 fstar; | |
44 | UINT256 P256; | |
200359e8 L |
45 | |
46 | // unpack x | |
b2a00c89 L |
47 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
48 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
49 | C1.w[1] = x.w[1] & MASK_COEFF; | |
50 | C1.w[0] = x.w[0]; | |
200359e8 L |
51 | |
52 | // check for NaN or Infinity | |
b2a00c89 | 53 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 54 | // x is special |
b2a00c89 L |
55 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
56 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
57 | // set invalid flag | |
58 | *pfpsf |= INVALID_EXCEPTION; | |
59 | // return Integer Indefinite | |
60 | res = 0x80000000; | |
61 | } else { // x is QNaN | |
62 | // set invalid flag | |
63 | *pfpsf |= INVALID_EXCEPTION; | |
64 | // return Integer Indefinite | |
65 | res = 0x80000000; | |
66 | } | |
67 | BID_RETURN (res); | |
68 | } else { // x is not a NaN, so it must be infinity | |
69 | if (!x_sign) { // x is +inf | |
70 | // set invalid flag | |
71 | *pfpsf |= INVALID_EXCEPTION; | |
72 | // return Integer Indefinite | |
73 | res = 0x80000000; | |
74 | } else { // x is -inf | |
75 | // set invalid flag | |
76 | *pfpsf |= INVALID_EXCEPTION; | |
77 | // return Integer Indefinite | |
78 | res = 0x80000000; | |
200359e8 | 79 | } |
b2a00c89 L |
80 | BID_RETURN (res); |
81 | } | |
82 | } | |
200359e8 | 83 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
84 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
85 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
86 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
87 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
88 | res = 0x00000000; | |
89 | BID_RETURN (res); | |
90 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
91 | // x is 0 | |
92 | res = 0x00000000; | |
93 | BID_RETURN (res); | |
94 | } else { // x is not special and is not zero | |
95 | ||
96 | // q = nr. of decimal digits in x | |
97 | // determine first the nr. of bits in x | |
98 | if (C1.w[1] == 0) { | |
99 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
100 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
101 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
102 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
103 | x_nr_bits = | |
104 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
105 | } else { // x < 2^32 | |
106 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
107 | x_nr_bits = |
108 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
109 | } | |
b2a00c89 L |
110 | } else { // if x < 2^53 |
111 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 112 | x_nr_bits = |
b2a00c89 | 113 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 114 | } |
b2a00c89 L |
115 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
116 | tmp1.d = (double) C1.w[1]; // exact conversion | |
117 | x_nr_bits = | |
118 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
119 | } | |
120 | q = nr_digits[x_nr_bits - 1].digits; | |
121 | if (q == 0) { | |
122 | q = nr_digits[x_nr_bits - 1].digits1; | |
123 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
124 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
125 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
126 | q++; | |
127 | } | |
128 | exp = (x_exp >> 49) - 6176; | |
129 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
130 | // set invalid flag | |
131 | *pfpsf |= INVALID_EXCEPTION; | |
132 | // return Integer Indefinite | |
133 | res = 0x80000000; | |
134 | BID_RETURN (res); | |
135 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
136 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
137 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
138 | // the cases that do not fit are identified here; the ones that fit | |
139 | // fall through and will be handled with other cases further, | |
140 | // under '1 <= q + exp <= 10' | |
141 | if (x_sign) { // if n < 0 and q + exp = 10 | |
142 | // if n < -2^31 - 1/2 then n is too large | |
143 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2 | |
144 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=34 | |
145 | if (q <= 11) { | |
146 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
147 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
148 | if (tmp64 > 0x500000005ull) { | |
149 | // set invalid flag | |
150 | *pfpsf |= INVALID_EXCEPTION; | |
151 | // return Integer Indefinite | |
152 | res = 0x80000000; | |
153 | BID_RETURN (res); | |
154 | } | |
155 | // else cases that can be rounded to a 32-bit int fall through | |
156 | // to '1 <= q + exp <= 10' | |
157 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
158 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005 <=> | |
159 | // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23 | |
160 | // (scale 2^31+1/2 up) | |
161 | tmp64 = 0x500000005ull; | |
162 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
163 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
164 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
165 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
166 | } | |
167 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
168 | // set invalid flag | |
169 | *pfpsf |= INVALID_EXCEPTION; | |
170 | // return Integer Indefinite | |
171 | res = 0x80000000; | |
172 | BID_RETURN (res); | |
173 | } | |
174 | // else cases that can be rounded to a 32-bit int fall through | |
175 | // to '1 <= q + exp <= 10' | |
176 | } | |
177 | } else { // if n > 0 and q + exp = 10 | |
178 | // if n >= 2^31 - 1/2 then n is too large | |
179 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2 | |
180 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34 | |
181 | if (q <= 11) { | |
182 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
183 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
184 | if (tmp64 >= 0x4fffffffbull) { | |
185 | // set invalid flag | |
186 | *pfpsf |= INVALID_EXCEPTION; | |
187 | // return Integer Indefinite | |
188 | res = 0x80000000; | |
189 | BID_RETURN (res); | |
190 | } | |
191 | // else cases that can be rounded to a 32-bit int fall through | |
192 | // to '1 <= q + exp <= 10' | |
193 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
194 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=> | |
195 | // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23 | |
196 | // (scale 2^31-1/2 up) | |
197 | tmp64 = 0x4fffffffbull; | |
198 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
199 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
200 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
201 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
202 | } | |
203 | if (C1.w[1] > C.w[1] | |
204 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
205 | // set invalid flag | |
206 | *pfpsf |= INVALID_EXCEPTION; | |
207 | // return Integer Indefinite | |
208 | res = 0x80000000; | |
209 | BID_RETURN (res); | |
210 | } | |
211 | // else cases that can be rounded to a 32-bit int fall through | |
212 | // to '1 <= q + exp <= 10' | |
213 | } | |
200359e8 | 214 | } |
b2a00c89 L |
215 | } |
216 | // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2 | |
217 | // Note: some of the cases tested for above fall through to this point | |
218 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) | |
219 | // return 0 | |
220 | res = 0x00000000; | |
221 | BID_RETURN (res); | |
222 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) | |
223 | // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1) | |
224 | // res = 0 | |
225 | // else | |
226 | // res = +/-1 | |
227 | ind = q - 1; | |
228 | if (ind <= 18) { // 0 <= ind <= 18 | |
229 | if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) { | |
230 | res = 0x00000000; // return 0 | |
231 | } else if (x_sign) { // n < 0 | |
232 | res = 0xffffffff; // return -1 | |
233 | } else { // n > 0 | |
234 | res = 0x00000001; // return +1 | |
235 | } | |
236 | } else { // 19 <= ind <= 33 | |
237 | if ((C1.w[1] < midpoint128[ind - 19].w[1]) | |
238 | || ((C1.w[1] == midpoint128[ind - 19].w[1]) | |
239 | && (C1.w[0] <= midpoint128[ind - 19].w[0]))) { | |
240 | res = 0x00000000; // return 0 | |
241 | } else if (x_sign) { // n < 0 | |
242 | res = 0xffffffff; // return -1 | |
243 | } else { // n > 0 | |
244 | res = 0x00000001; // return +1 | |
200359e8 L |
245 | } |
246 | } | |
b2a00c89 L |
247 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) |
248 | // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded | |
249 | // to nearest to a 32-bit signed integer | |
250 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
251 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
252 | // chop off ind digits from the lower part of C1 | |
253 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
254 | tmp64 = C1.w[0]; | |
255 | if (ind <= 19) { | |
256 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
257 | } else { | |
258 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
259 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
260 | } | |
261 | if (C1.w[0] < tmp64) | |
262 | C1.w[1]++; | |
263 | // calculate C* and f* | |
264 | // C* is actually floor(C*) in this case | |
265 | // C* and f* need shifting and masking, as shown by | |
266 | // shiftright128[] and maskhigh128[] | |
267 | // 1 <= x <= 33 | |
268 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
269 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
270 | // the approximation of 10^(-x) was rounded up to 118 bits | |
271 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
272 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
273 | Cstar.w[1] = P256.w[3]; | |
274 | Cstar.w[0] = P256.w[2]; | |
275 | fstar.w[3] = 0; | |
276 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
277 | fstar.w[1] = P256.w[1]; | |
278 | fstar.w[0] = P256.w[0]; | |
279 | } else { // 22 <= ind - 1 <= 33 | |
280 | Cstar.w[1] = 0; | |
281 | Cstar.w[0] = P256.w[3]; | |
282 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
283 | fstar.w[2] = P256.w[2]; | |
284 | fstar.w[1] = P256.w[1]; | |
285 | fstar.w[0] = P256.w[0]; | |
286 | } | |
287 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
288 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
289 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
290 | // if floor(C*) is even then C* = floor(C*) - logical right | |
291 | // shift; C* has p decimal digits, correct by Prop. 1) | |
292 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
293 | // shift; C* has p decimal digits, correct by Pr. 1) | |
200359e8 | 294 | // else |
b2a00c89 L |
295 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
296 | // correct by Property 1) | |
297 | // n = C* * 10^(e+x) | |
298 | ||
299 | // shift right C* by Ex-128 = shiftright128[ind] | |
300 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
301 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
302 | Cstar.w[0] = | |
303 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
304 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
305 | } else { // 22 <= ind - 1 <= 33 | |
306 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
307 | } | |
308 | // if the result was a midpoint it was rounded away from zero, so | |
309 | // it will need a correction | |
310 | // check for midpoints | |
311 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
312 | && (fstar.w[1] || fstar.w[0]) | |
313 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
314 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
315 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
316 | // the result is a midpoint; round to nearest | |
317 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
318 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
319 | Cstar.w[0]--; // Cstar.w[0] is now even | |
320 | } // else MP in [ODD, EVEN] | |
200359e8 | 321 | } |
b2a00c89 L |
322 | if (x_sign) |
323 | res = -Cstar.w[0]; | |
324 | else | |
325 | res = Cstar.w[0]; | |
326 | } else if (exp == 0) { | |
327 | // 1 <= q <= 10 | |
328 | // res = +/-C (exact) | |
329 | if (x_sign) | |
330 | res = -C1.w[0]; | |
331 | else | |
332 | res = C1.w[0]; | |
333 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
334 | // res = +/-C * 10^exp (exact) | |
335 | if (x_sign) | |
336 | res = -C1.w[0] * ten2k64[exp]; | |
337 | else | |
338 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
339 | } |
340 | } | |
b2a00c89 L |
341 | } |
342 | ||
343 | BID_RETURN (res); | |
200359e8 L |
344 | } |
345 | ||
346 | /***************************************************************************** | |
347 | * BID128_to_int32_xrnint | |
348 | ****************************************************************************/ | |
349 | ||
b2a00c89 L |
350 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrnint, |
351 | x) | |
200359e8 | 352 | |
b2a00c89 L |
353 | int res; |
354 | UINT64 x_sign; | |
355 | UINT64 x_exp; | |
356 | int exp; // unbiased exponent | |
200359e8 | 357 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
358 | UINT64 tmp64, tmp64A; |
359 | BID_UI64DOUBLE tmp1; | |
360 | unsigned int x_nr_bits; | |
361 | int q, ind, shift; | |
362 | UINT128 C1, C; | |
363 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
364 | UINT256 fstar; | |
365 | UINT256 P256; | |
200359e8 L |
366 | |
367 | // unpack x | |
b2a00c89 L |
368 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
369 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
370 | C1.w[1] = x.w[1] & MASK_COEFF; | |
371 | C1.w[0] = x.w[0]; | |
200359e8 L |
372 | |
373 | // check for NaN or Infinity | |
b2a00c89 | 374 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 375 | // x is special |
b2a00c89 L |
376 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
377 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
378 | // set invalid flag | |
379 | *pfpsf |= INVALID_EXCEPTION; | |
380 | // return Integer Indefinite | |
381 | res = 0x80000000; | |
382 | } else { // x is QNaN | |
383 | // set invalid flag | |
384 | *pfpsf |= INVALID_EXCEPTION; | |
385 | // return Integer Indefinite | |
386 | res = 0x80000000; | |
200359e8 | 387 | } |
b2a00c89 L |
388 | BID_RETURN (res); |
389 | } else { // x is not a NaN, so it must be infinity | |
390 | if (!x_sign) { // x is +inf | |
391 | // set invalid flag | |
392 | *pfpsf |= INVALID_EXCEPTION; | |
393 | // return Integer Indefinite | |
394 | res = 0x80000000; | |
395 | } else { // x is -inf | |
396 | // set invalid flag | |
397 | *pfpsf |= INVALID_EXCEPTION; | |
398 | // return Integer Indefinite | |
399 | res = 0x80000000; | |
400 | } | |
401 | BID_RETURN (res); | |
402 | } | |
403 | } | |
200359e8 | 404 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
405 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
406 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
407 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
408 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
409 | res = 0x00000000; | |
410 | BID_RETURN (res); | |
411 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
412 | // x is 0 | |
413 | res = 0x00000000; | |
414 | BID_RETURN (res); | |
415 | } else { // x is not special and is not zero | |
416 | ||
417 | // q = nr. of decimal digits in x | |
418 | // determine first the nr. of bits in x | |
419 | if (C1.w[1] == 0) { | |
420 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
421 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
422 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
423 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
424 | x_nr_bits = | |
425 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
426 | } else { // x < 2^32 | |
427 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
428 | x_nr_bits = |
429 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
430 | } | |
b2a00c89 L |
431 | } else { // if x < 2^53 |
432 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 433 | x_nr_bits = |
b2a00c89 | 434 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 435 | } |
b2a00c89 L |
436 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
437 | tmp1.d = (double) C1.w[1]; // exact conversion | |
438 | x_nr_bits = | |
439 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
440 | } | |
441 | q = nr_digits[x_nr_bits - 1].digits; | |
442 | if (q == 0) { | |
443 | q = nr_digits[x_nr_bits - 1].digits1; | |
444 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
445 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
446 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
447 | q++; | |
448 | } | |
449 | exp = (x_exp >> 49) - 6176; | |
450 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
451 | // set invalid flag | |
452 | *pfpsf |= INVALID_EXCEPTION; | |
453 | // return Integer Indefinite | |
454 | res = 0x80000000; | |
455 | BID_RETURN (res); | |
456 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
457 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
458 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
459 | // the cases that do not fit are identified here; the ones that fit | |
460 | // fall through and will be handled with other cases further, | |
461 | // under '1 <= q + exp <= 10' | |
462 | if (x_sign) { // if n < 0 and q + exp = 10 | |
463 | // if n < -2^31 - 1/2 then n is too large | |
464 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2 | |
465 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=34 | |
466 | if (q <= 11) { | |
467 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
468 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
469 | if (tmp64 > 0x500000005ull) { | |
470 | // set invalid flag | |
471 | *pfpsf |= INVALID_EXCEPTION; | |
472 | // return Integer Indefinite | |
473 | res = 0x80000000; | |
474 | BID_RETURN (res); | |
475 | } | |
476 | // else cases that can be rounded to a 32-bit int fall through | |
477 | // to '1 <= q + exp <= 10' | |
478 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
479 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005 <=> | |
480 | // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23 | |
481 | // (scale 2^31+1/2 up) | |
482 | tmp64 = 0x500000005ull; | |
483 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
484 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
485 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
486 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
487 | } | |
488 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
489 | // set invalid flag | |
490 | *pfpsf |= INVALID_EXCEPTION; | |
491 | // return Integer Indefinite | |
492 | res = 0x80000000; | |
493 | BID_RETURN (res); | |
494 | } | |
495 | // else cases that can be rounded to a 32-bit int fall through | |
496 | // to '1 <= q + exp <= 10' | |
497 | } | |
498 | } else { // if n > 0 and q + exp = 10 | |
499 | // if n >= 2^31 - 1/2 then n is too large | |
500 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2 | |
501 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34 | |
502 | if (q <= 11) { | |
503 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
504 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
505 | if (tmp64 >= 0x4fffffffbull) { | |
506 | // set invalid flag | |
507 | *pfpsf |= INVALID_EXCEPTION; | |
508 | // return Integer Indefinite | |
509 | res = 0x80000000; | |
510 | BID_RETURN (res); | |
511 | } | |
512 | // else cases that can be rounded to a 32-bit int fall through | |
513 | // to '1 <= q + exp <= 10' | |
514 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
515 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=> | |
516 | // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23 | |
517 | // (scale 2^31-1/2 up) | |
518 | tmp64 = 0x4fffffffbull; | |
519 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
520 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
521 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
522 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
523 | } | |
524 | if (C1.w[1] > C.w[1] | |
525 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
526 | // set invalid flag | |
527 | *pfpsf |= INVALID_EXCEPTION; | |
528 | // return Integer Indefinite | |
529 | res = 0x80000000; | |
530 | BID_RETURN (res); | |
531 | } | |
532 | // else cases that can be rounded to a 32-bit int fall through | |
533 | // to '1 <= q + exp <= 10' | |
534 | } | |
200359e8 | 535 | } |
b2a00c89 L |
536 | } |
537 | // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2 | |
538 | // Note: some of the cases tested for above fall through to this point | |
539 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) | |
540 | // set inexact flag | |
541 | *pfpsf |= INEXACT_EXCEPTION; | |
542 | // return 0 | |
543 | res = 0x00000000; | |
544 | BID_RETURN (res); | |
545 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) | |
546 | // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1) | |
547 | // res = 0 | |
548 | // else | |
549 | // res = +/-1 | |
550 | ind = q - 1; | |
551 | if (ind <= 18) { // 0 <= ind <= 18 | |
552 | if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) { | |
553 | res = 0x00000000; // return 0 | |
554 | } else if (x_sign) { // n < 0 | |
555 | res = 0xffffffff; // return -1 | |
556 | } else { // n > 0 | |
557 | res = 0x00000001; // return +1 | |
558 | } | |
559 | } else { // 19 <= ind <= 33 | |
560 | if ((C1.w[1] < midpoint128[ind - 19].w[1]) | |
561 | || ((C1.w[1] == midpoint128[ind - 19].w[1]) | |
562 | && (C1.w[0] <= midpoint128[ind - 19].w[0]))) { | |
563 | res = 0x00000000; // return 0 | |
564 | } else if (x_sign) { // n < 0 | |
565 | res = 0xffffffff; // return -1 | |
566 | } else { // n > 0 | |
567 | res = 0x00000001; // return +1 | |
200359e8 L |
568 | } |
569 | } | |
b2a00c89 L |
570 | // set inexact flag |
571 | *pfpsf |= INEXACT_EXCEPTION; | |
572 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
573 | // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded | |
574 | // to nearest to a 32-bit signed integer | |
575 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
576 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
577 | // chop off ind digits from the lower part of C1 | |
578 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
579 | tmp64 = C1.w[0]; | |
580 | if (ind <= 19) { | |
581 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
582 | } else { | |
583 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
584 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
585 | } | |
586 | if (C1.w[0] < tmp64) | |
587 | C1.w[1]++; | |
588 | // calculate C* and f* | |
589 | // C* is actually floor(C*) in this case | |
590 | // C* and f* need shifting and masking, as shown by | |
591 | // shiftright128[] and maskhigh128[] | |
592 | // 1 <= x <= 33 | |
593 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
594 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
595 | // the approximation of 10^(-x) was rounded up to 118 bits | |
596 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
597 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
598 | Cstar.w[1] = P256.w[3]; | |
599 | Cstar.w[0] = P256.w[2]; | |
600 | fstar.w[3] = 0; | |
601 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
602 | fstar.w[1] = P256.w[1]; | |
603 | fstar.w[0] = P256.w[0]; | |
604 | } else { // 22 <= ind - 1 <= 33 | |
605 | Cstar.w[1] = 0; | |
606 | Cstar.w[0] = P256.w[3]; | |
607 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
608 | fstar.w[2] = P256.w[2]; | |
609 | fstar.w[1] = P256.w[1]; | |
610 | fstar.w[0] = P256.w[0]; | |
611 | } | |
612 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
613 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
614 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
615 | // if floor(C*) is even then C* = floor(C*) - logical right | |
616 | // shift; C* has p decimal digits, correct by Prop. 1) | |
617 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
618 | // shift; C* has p decimal digits, correct by Pr. 1) | |
200359e8 | 619 | // else |
b2a00c89 L |
620 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
621 | // correct by Property 1) | |
622 | // n = C* * 10^(e+x) | |
623 | ||
624 | // shift right C* by Ex-128 = shiftright128[ind] | |
625 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
626 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
627 | Cstar.w[0] = | |
628 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
629 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
630 | } else { // 22 <= ind - 1 <= 33 | |
631 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
632 | } | |
633 | // determine inexactness of the rounding of C* | |
634 | // if (0 < f* - 1/2 < 10^(-x)) then | |
635 | // the result is exact | |
636 | // else // if (f* - 1/2 > T*) then | |
637 | // the result is inexact | |
638 | if (ind - 1 <= 2) { | |
639 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
640 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
641 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
642 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
643 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
644 | // set the inexact flag |
645 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
646 | } // else the result is exact |
647 | } else { // the result is inexact; f2* <= 1/2 | |
648 | // set the inexact flag | |
649 | *pfpsf |= INEXACT_EXCEPTION; | |
650 | } | |
651 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
652 | if (fstar.w[3] > 0x0 || | |
653 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
654 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
655 | (fstar.w[1] || fstar.w[0]))) { | |
656 | // f2* > 1/2 and the result may be exact | |
657 | // Calculate f2* - 1/2 | |
658 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
659 | tmp64A = fstar.w[3]; | |
660 | if (tmp64 > fstar.w[2]) | |
661 | tmp64A--; | |
662 | if (tmp64A || tmp64 | |
663 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
664 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
665 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
666 | // set the inexact flag |
667 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
668 | } // else the result is exact |
669 | } else { // the result is inexact; f2* <= 1/2 | |
670 | // set the inexact flag | |
671 | *pfpsf |= INEXACT_EXCEPTION; | |
672 | } | |
673 | } else { // if 22 <= ind <= 33 | |
674 | if (fstar.w[3] > onehalf128[ind - 1] || | |
675 | (fstar.w[3] == onehalf128[ind - 1] && | |
676 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
677 | // f2* > 1/2 and the result may be exact | |
678 | // Calculate f2* - 1/2 | |
679 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
680 | if (tmp64 || fstar.w[2] | |
681 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
682 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
683 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
684 | // set the inexact flag |
685 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
686 | } // else the result is exact |
687 | } else { // the result is inexact; f2* <= 1/2 | |
688 | // set the inexact flag | |
689 | *pfpsf |= INEXACT_EXCEPTION; | |
690 | } | |
200359e8 | 691 | } |
b2a00c89 L |
692 | // if the result was a midpoint it was rounded away from zero, so |
693 | // it will need a correction | |
694 | // check for midpoints | |
695 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
696 | && (fstar.w[1] || fstar.w[0]) | |
697 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
698 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
699 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
700 | // the result is a midpoint; round to nearest | |
701 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
702 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
703 | Cstar.w[0]--; // Cstar.w[0] is now even | |
704 | } // else MP in [ODD, EVEN] | |
705 | } | |
706 | if (x_sign) | |
707 | res = -Cstar.w[0]; | |
708 | else | |
709 | res = Cstar.w[0]; | |
710 | } else if (exp == 0) { | |
711 | // 1 <= q <= 10 | |
712 | // res = +/-C (exact) | |
713 | if (x_sign) | |
714 | res = -C1.w[0]; | |
715 | else | |
716 | res = C1.w[0]; | |
717 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
718 | // res = +/-C * 10^exp (exact) | |
719 | if (x_sign) | |
720 | res = -C1.w[0] * ten2k64[exp]; | |
721 | else | |
722 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
723 | } |
724 | } | |
b2a00c89 L |
725 | } |
726 | ||
727 | BID_RETURN (res); | |
200359e8 L |
728 | } |
729 | ||
730 | /***************************************************************************** | |
731 | * BID128_to_int32_floor | |
732 | ****************************************************************************/ | |
733 | ||
b2a00c89 | 734 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_floor, x) |
200359e8 | 735 | |
b2a00c89 L |
736 | int res; |
737 | UINT64 x_sign; | |
738 | UINT64 x_exp; | |
739 | int exp; // unbiased exponent | |
200359e8 | 740 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
741 | UINT64 tmp64, tmp64A; |
742 | BID_UI64DOUBLE tmp1; | |
743 | unsigned int x_nr_bits; | |
744 | int q, ind, shift; | |
745 | UINT128 C1, C; | |
746 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
747 | UINT256 fstar; | |
748 | UINT256 P256; | |
749 | int is_inexact_lt_midpoint = 0; | |
750 | int is_inexact_gt_midpoint = 0; | |
751 | int is_midpoint_lt_even = 0; | |
752 | int is_midpoint_gt_even = 0; | |
200359e8 L |
753 | |
754 | // unpack x | |
b2a00c89 L |
755 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
756 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
757 | C1.w[1] = x.w[1] & MASK_COEFF; | |
758 | C1.w[0] = x.w[0]; | |
200359e8 L |
759 | |
760 | // check for NaN or Infinity | |
b2a00c89 | 761 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 762 | // x is special |
b2a00c89 L |
763 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
764 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
765 | // set invalid flag | |
766 | *pfpsf |= INVALID_EXCEPTION; | |
767 | // return Integer Indefinite | |
768 | res = 0x80000000; | |
769 | } else { // x is QNaN | |
770 | // set invalid flag | |
771 | *pfpsf |= INVALID_EXCEPTION; | |
772 | // return Integer Indefinite | |
773 | res = 0x80000000; | |
774 | } | |
775 | BID_RETURN (res); | |
776 | } else { // x is not a NaN, so it must be infinity | |
777 | if (!x_sign) { // x is +inf | |
778 | // set invalid flag | |
779 | *pfpsf |= INVALID_EXCEPTION; | |
780 | // return Integer Indefinite | |
781 | res = 0x80000000; | |
782 | } else { // x is -inf | |
783 | // set invalid flag | |
784 | *pfpsf |= INVALID_EXCEPTION; | |
785 | // return Integer Indefinite | |
786 | res = 0x80000000; | |
200359e8 | 787 | } |
b2a00c89 L |
788 | BID_RETURN (res); |
789 | } | |
790 | } | |
200359e8 | 791 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
792 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
793 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
794 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
795 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
796 | res = 0x00000000; | |
797 | BID_RETURN (res); | |
798 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
799 | // x is 0 | |
800 | res = 0x00000000; | |
801 | BID_RETURN (res); | |
802 | } else { // x is not special and is not zero | |
803 | ||
804 | // q = nr. of decimal digits in x | |
805 | // determine first the nr. of bits in x | |
806 | if (C1.w[1] == 0) { | |
807 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
808 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
809 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
810 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
811 | x_nr_bits = | |
812 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
813 | } else { // x < 2^32 | |
814 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
815 | x_nr_bits = |
816 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
817 | } | |
b2a00c89 L |
818 | } else { // if x < 2^53 |
819 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 820 | x_nr_bits = |
b2a00c89 | 821 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 822 | } |
b2a00c89 L |
823 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
824 | tmp1.d = (double) C1.w[1]; // exact conversion | |
825 | x_nr_bits = | |
826 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
827 | } | |
828 | q = nr_digits[x_nr_bits - 1].digits; | |
829 | if (q == 0) { | |
830 | q = nr_digits[x_nr_bits - 1].digits1; | |
831 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
832 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
833 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
834 | q++; | |
835 | } | |
836 | exp = (x_exp >> 49) - 6176; | |
837 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
838 | // set invalid flag | |
839 | *pfpsf |= INVALID_EXCEPTION; | |
840 | // return Integer Indefinite | |
841 | res = 0x80000000; | |
842 | BID_RETURN (res); | |
843 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
844 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
845 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
846 | // the cases that do not fit are identified here; the ones that fit | |
847 | // fall through and will be handled with other cases further, | |
848 | // under '1 <= q + exp <= 10' | |
849 | if (x_sign) { // if n < 0 and q + exp = 10 | |
850 | // if n < -2^31 then n is too large | |
851 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 | |
852 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=34 | |
853 | if (q <= 11) { | |
854 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
855 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
856 | if (tmp64 > 0x500000000ull) { | |
857 | // set invalid flag | |
858 | *pfpsf |= INVALID_EXCEPTION; | |
859 | // return Integer Indefinite | |
860 | res = 0x80000000; | |
861 | BID_RETURN (res); | |
862 | } | |
863 | // else cases that can be rounded to a 32-bit int fall through | |
864 | // to '1 <= q + exp <= 10' | |
865 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
866 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000 <=> | |
867 | // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
868 | // (scale 2^31 up) | |
869 | tmp64 = 0x500000000ull; | |
870 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
871 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
872 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
873 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
874 | } | |
875 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
876 | // set invalid flag | |
877 | *pfpsf |= INVALID_EXCEPTION; | |
878 | // return Integer Indefinite | |
879 | res = 0x80000000; | |
880 | BID_RETURN (res); | |
881 | } | |
882 | // else cases that can be rounded to a 32-bit int fall through | |
883 | // to '1 <= q + exp <= 10' | |
884 | } | |
885 | } else { // if n > 0 and q + exp = 10 | |
886 | // if n >= 2^31 then n is too large | |
887 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31 | |
888 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34 | |
889 | if (q <= 11) { | |
890 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
891 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
892 | if (tmp64 >= 0x500000000ull) { | |
893 | // set invalid flag | |
894 | *pfpsf |= INVALID_EXCEPTION; | |
895 | // return Integer Indefinite | |
896 | res = 0x80000000; | |
897 | BID_RETURN (res); | |
898 | } | |
899 | // else cases that can be rounded to a 32-bit int fall through | |
900 | // to '1 <= q + exp <= 10' | |
901 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
902 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=> | |
903 | // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
904 | // (scale 2^31 up) | |
905 | tmp64 = 0x500000000ull; | |
906 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
907 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
908 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
909 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
910 | } | |
911 | if (C1.w[1] > C.w[1] | |
912 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
913 | // set invalid flag | |
914 | *pfpsf |= INVALID_EXCEPTION; | |
915 | // return Integer Indefinite | |
916 | res = 0x80000000; | |
917 | BID_RETURN (res); | |
918 | } | |
919 | // else cases that can be rounded to a 32-bit int fall through | |
920 | // to '1 <= q + exp <= 10' | |
921 | } | |
200359e8 | 922 | } |
b2a00c89 L |
923 | } |
924 | // n is not too large to be converted to int32: -2^31 <= n < 2^31 | |
925 | // Note: some of the cases tested for above fall through to this point | |
926 | if ((q + exp) <= 0) { | |
927 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
928 | // return 0 | |
929 | if (x_sign) | |
930 | res = 0xffffffff; | |
931 | else | |
932 | res = 0x00000000; | |
933 | BID_RETURN (res); | |
934 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
935 | // -2^31 <= x <= -1 or 1 <= x < 2^31 so x can be rounded | |
936 | // toward negative infinity to a 32-bit signed integer | |
937 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
938 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
939 | // chop off ind digits from the lower part of C1 | |
940 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
941 | tmp64 = C1.w[0]; | |
942 | if (ind <= 19) { | |
943 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
944 | } else { | |
945 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
946 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
947 | } | |
948 | if (C1.w[0] < tmp64) | |
949 | C1.w[1]++; | |
950 | // calculate C* and f* | |
951 | // C* is actually floor(C*) in this case | |
952 | // C* and f* need shifting and masking, as shown by | |
953 | // shiftright128[] and maskhigh128[] | |
954 | // 1 <= x <= 33 | |
955 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
956 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
957 | // the approximation of 10^(-x) was rounded up to 118 bits | |
958 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
959 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
960 | Cstar.w[1] = P256.w[3]; | |
961 | Cstar.w[0] = P256.w[2]; | |
962 | fstar.w[3] = 0; | |
963 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
964 | fstar.w[1] = P256.w[1]; | |
965 | fstar.w[0] = P256.w[0]; | |
966 | } else { // 22 <= ind - 1 <= 33 | |
967 | Cstar.w[1] = 0; | |
968 | Cstar.w[0] = P256.w[3]; | |
969 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
970 | fstar.w[2] = P256.w[2]; | |
971 | fstar.w[1] = P256.w[1]; | |
972 | fstar.w[0] = P256.w[0]; | |
973 | } | |
974 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
975 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
976 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
977 | // if floor(C*) is even then C* = floor(C*) - logical right | |
978 | // shift; C* has p decimal digits, correct by Prop. 1) | |
979 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
980 | // shift; C* has p decimal digits, correct by Pr. 1) | |
981 | // else | |
982 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
983 | // correct by Property 1) | |
984 | // n = C* * 10^(e+x) | |
985 | ||
986 | // shift right C* by Ex-128 = shiftright128[ind] | |
987 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
988 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
989 | Cstar.w[0] = | |
990 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
991 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
992 | } else { // 22 <= ind - 1 <= 33 | |
993 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
994 | } | |
995 | // determine inexactness of the rounding of C* | |
996 | // if (0 < f* - 1/2 < 10^(-x)) then | |
997 | // the result is exact | |
998 | // else // if (f* - 1/2 > T*) then | |
999 | // the result is inexact | |
1000 | if (ind - 1 <= 2) { | |
1001 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
1002 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
1003 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
1004 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
1005 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
1006 | is_inexact_lt_midpoint = 1; | |
1007 | } // else the result is exact | |
1008 | } else { // the result is inexact; f2* <= 1/2 | |
1009 | is_inexact_gt_midpoint = 1; | |
1010 | } | |
1011 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
1012 | if (fstar.w[3] > 0x0 || | |
1013 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
1014 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
1015 | (fstar.w[1] || fstar.w[0]))) { | |
1016 | // f2* > 1/2 and the result may be exact | |
1017 | // Calculate f2* - 1/2 | |
1018 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
1019 | tmp64A = fstar.w[3]; | |
1020 | if (tmp64 > fstar.w[2]) | |
1021 | tmp64A--; | |
1022 | if (tmp64A || tmp64 | |
1023 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1024 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1025 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
1026 | is_inexact_lt_midpoint = 1; | |
1027 | } // else the result is exact | |
1028 | } else { // the result is inexact; f2* <= 1/2 | |
1029 | is_inexact_gt_midpoint = 1; | |
1030 | } | |
1031 | } else { // if 22 <= ind <= 33 | |
1032 | if (fstar.w[3] > onehalf128[ind - 1] || | |
1033 | (fstar.w[3] == onehalf128[ind - 1] && | |
1034 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
1035 | // f2* > 1/2 and the result may be exact | |
1036 | // Calculate f2* - 1/2 | |
1037 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
1038 | if (tmp64 || fstar.w[2] | |
1039 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1040 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1041 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
1042 | is_inexact_lt_midpoint = 1; | |
1043 | } // else the result is exact | |
1044 | } else { // the result is inexact; f2* <= 1/2 | |
1045 | is_inexact_gt_midpoint = 1; | |
200359e8 L |
1046 | } |
1047 | } | |
b2a00c89 L |
1048 | |
1049 | // if the result was a midpoint it was rounded away from zero, so | |
1050 | // it will need a correction | |
1051 | // check for midpoints | |
1052 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
1053 | && (fstar.w[1] || fstar.w[0]) | |
1054 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
1055 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1056 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
1057 | // the result is a midpoint; round to nearest | |
1058 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
1059 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
1060 | Cstar.w[0]--; // Cstar.w[0] is now even | |
1061 | is_midpoint_gt_even = 1; | |
1062 | is_inexact_lt_midpoint = 0; | |
1063 | is_inexact_gt_midpoint = 0; | |
1064 | } else { // else MP in [ODD, EVEN] | |
1065 | is_midpoint_lt_even = 1; | |
1066 | is_inexact_lt_midpoint = 0; | |
1067 | is_inexact_gt_midpoint = 0; | |
1068 | } | |
1069 | } | |
1070 | // general correction for RM | |
1071 | if (x_sign && (is_midpoint_gt_even || is_inexact_lt_midpoint)) { | |
1072 | Cstar.w[0] = Cstar.w[0] + 1; | |
1073 | } else if (!x_sign | |
1074 | && (is_midpoint_lt_even || is_inexact_gt_midpoint)) { | |
1075 | Cstar.w[0] = Cstar.w[0] - 1; | |
1076 | } else { | |
1077 | ; // the result is already correct | |
1078 | } | |
200359e8 | 1079 | if (x_sign) |
b2a00c89 | 1080 | res = -Cstar.w[0]; |
200359e8 | 1081 | else |
b2a00c89 L |
1082 | res = Cstar.w[0]; |
1083 | } else if (exp == 0) { | |
1084 | // 1 <= q <= 10 | |
1085 | // res = +/-C (exact) | |
1086 | if (x_sign) | |
1087 | res = -C1.w[0]; | |
1088 | else | |
1089 | res = C1.w[0]; | |
1090 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
1091 | // res = +/-C * 10^exp (exact) | |
1092 | if (x_sign) | |
1093 | res = -C1.w[0] * ten2k64[exp]; | |
1094 | else | |
1095 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
1096 | } |
1097 | } | |
b2a00c89 L |
1098 | } |
1099 | ||
1100 | BID_RETURN (res); | |
200359e8 L |
1101 | } |
1102 | ||
1103 | ||
1104 | /***************************************************************************** | |
1105 | * BID128_to_int32_xfloor | |
1106 | ****************************************************************************/ | |
1107 | ||
b2a00c89 L |
1108 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xfloor, |
1109 | x) | |
200359e8 | 1110 | |
b2a00c89 L |
1111 | int res; |
1112 | UINT64 x_sign; | |
1113 | UINT64 x_exp; | |
1114 | int exp; // unbiased exponent | |
200359e8 | 1115 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
1116 | UINT64 tmp64, tmp64A; |
1117 | BID_UI64DOUBLE tmp1; | |
1118 | unsigned int x_nr_bits; | |
1119 | int q, ind, shift; | |
1120 | UINT128 C1, C; | |
1121 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
1122 | UINT256 fstar; | |
1123 | UINT256 P256; | |
1124 | int is_inexact_lt_midpoint = 0; | |
1125 | int is_inexact_gt_midpoint = 0; | |
1126 | int is_midpoint_lt_even = 0; | |
1127 | int is_midpoint_gt_even = 0; | |
200359e8 L |
1128 | |
1129 | // unpack x | |
b2a00c89 L |
1130 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
1131 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
1132 | C1.w[1] = x.w[1] & MASK_COEFF; | |
1133 | C1.w[0] = x.w[0]; | |
200359e8 L |
1134 | |
1135 | // check for NaN or Infinity | |
b2a00c89 | 1136 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 1137 | // x is special |
b2a00c89 L |
1138 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
1139 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
1140 | // set invalid flag | |
1141 | *pfpsf |= INVALID_EXCEPTION; | |
1142 | // return Integer Indefinite | |
1143 | res = 0x80000000; | |
1144 | } else { // x is QNaN | |
1145 | // set invalid flag | |
1146 | *pfpsf |= INVALID_EXCEPTION; | |
1147 | // return Integer Indefinite | |
1148 | res = 0x80000000; | |
1149 | } | |
1150 | BID_RETURN (res); | |
1151 | } else { // x is not a NaN, so it must be infinity | |
1152 | if (!x_sign) { // x is +inf | |
1153 | // set invalid flag | |
1154 | *pfpsf |= INVALID_EXCEPTION; | |
1155 | // return Integer Indefinite | |
1156 | res = 0x80000000; | |
1157 | } else { // x is -inf | |
1158 | // set invalid flag | |
1159 | *pfpsf |= INVALID_EXCEPTION; | |
1160 | // return Integer Indefinite | |
1161 | res = 0x80000000; | |
200359e8 | 1162 | } |
b2a00c89 L |
1163 | BID_RETURN (res); |
1164 | } | |
1165 | } | |
200359e8 | 1166 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
1167 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
1168 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
1169 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
1170 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
1171 | res = 0x00000000; | |
1172 | BID_RETURN (res); | |
1173 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
1174 | // x is 0 | |
1175 | res = 0x00000000; | |
1176 | BID_RETURN (res); | |
1177 | } else { // x is not special and is not zero | |
1178 | ||
1179 | // q = nr. of decimal digits in x | |
1180 | // determine first the nr. of bits in x | |
1181 | if (C1.w[1] == 0) { | |
1182 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
1183 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
1184 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
1185 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
1186 | x_nr_bits = | |
1187 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1188 | } else { // x < 2^32 | |
1189 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
1190 | x_nr_bits = |
1191 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1192 | } | |
b2a00c89 L |
1193 | } else { // if x < 2^53 |
1194 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 1195 | x_nr_bits = |
b2a00c89 | 1196 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 1197 | } |
b2a00c89 L |
1198 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
1199 | tmp1.d = (double) C1.w[1]; // exact conversion | |
1200 | x_nr_bits = | |
1201 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1202 | } | |
1203 | q = nr_digits[x_nr_bits - 1].digits; | |
1204 | if (q == 0) { | |
1205 | q = nr_digits[x_nr_bits - 1].digits1; | |
1206 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
1207 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
1208 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
1209 | q++; | |
1210 | } | |
1211 | exp = (x_exp >> 49) - 6176; | |
1212 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
1213 | // set invalid flag | |
1214 | *pfpsf |= INVALID_EXCEPTION; | |
1215 | // return Integer Indefinite | |
1216 | res = 0x80000000; | |
1217 | BID_RETURN (res); | |
1218 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
1219 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
1220 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
1221 | // the cases that do not fit are identified here; the ones that fit | |
1222 | // fall through and will be handled with other cases further, | |
1223 | // under '1 <= q + exp <= 10' | |
1224 | if (x_sign) { // if n < 0 and q + exp = 10 | |
1225 | // if n < -2^31 then n is too large | |
1226 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 | |
1227 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=34 | |
1228 | if (q <= 11) { | |
1229 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
1230 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
1231 | if (tmp64 > 0x500000000ull) { | |
1232 | // set invalid flag | |
1233 | *pfpsf |= INVALID_EXCEPTION; | |
1234 | // return Integer Indefinite | |
1235 | res = 0x80000000; | |
1236 | BID_RETURN (res); | |
1237 | } | |
1238 | // else cases that can be rounded to a 32-bit int fall through | |
1239 | // to '1 <= q + exp <= 10' | |
1240 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
1241 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000 <=> | |
1242 | // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
1243 | // (scale 2^31 up) | |
1244 | tmp64 = 0x500000000ull; | |
1245 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
1246 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
1247 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
1248 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
1249 | } | |
1250 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
1251 | // set invalid flag | |
1252 | *pfpsf |= INVALID_EXCEPTION; | |
1253 | // return Integer Indefinite | |
1254 | res = 0x80000000; | |
1255 | BID_RETURN (res); | |
1256 | } | |
1257 | // else cases that can be rounded to a 32-bit int fall through | |
1258 | // to '1 <= q + exp <= 10' | |
1259 | } | |
1260 | } else { // if n > 0 and q + exp = 10 | |
1261 | // if n >= 2^31 then n is too large | |
1262 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31 | |
1263 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34 | |
1264 | if (q <= 11) { | |
1265 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
1266 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
1267 | if (tmp64 >= 0x500000000ull) { | |
1268 | // set invalid flag | |
1269 | *pfpsf |= INVALID_EXCEPTION; | |
1270 | // return Integer Indefinite | |
1271 | res = 0x80000000; | |
1272 | BID_RETURN (res); | |
1273 | } | |
1274 | // else cases that can be rounded to a 32-bit int fall through | |
1275 | // to '1 <= q + exp <= 10' | |
1276 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
1277 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=> | |
1278 | // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
1279 | // (scale 2^31 up) | |
1280 | tmp64 = 0x500000000ull; | |
1281 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
1282 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
1283 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
1284 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
1285 | } | |
1286 | if (C1.w[1] > C.w[1] | |
1287 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
1288 | // set invalid flag | |
1289 | *pfpsf |= INVALID_EXCEPTION; | |
1290 | // return Integer Indefinite | |
1291 | res = 0x80000000; | |
1292 | BID_RETURN (res); | |
1293 | } | |
1294 | // else cases that can be rounded to a 32-bit int fall through | |
1295 | // to '1 <= q + exp <= 10' | |
200359e8 L |
1296 | } |
1297 | } | |
b2a00c89 L |
1298 | } |
1299 | // n is not too large to be converted to int32: -2^31 <= n < 2^31 | |
1300 | // Note: some of the cases tested for above fall through to this point | |
1301 | if ((q + exp) <= 0) { | |
1302 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
1303 | // set inexact flag | |
1304 | *pfpsf |= INEXACT_EXCEPTION; | |
1305 | // return 0 | |
1306 | if (x_sign) | |
1307 | res = 0xffffffff; | |
1308 | else | |
1309 | res = 0x00000000; | |
1310 | BID_RETURN (res); | |
1311 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
1312 | // -2^31 <= x <= -1 or 1 <= x < 2^31 so x can be rounded | |
1313 | // toward negative infinity to a 32-bit signed integer | |
1314 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
1315 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
1316 | // chop off ind digits from the lower part of C1 | |
1317 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
1318 | tmp64 = C1.w[0]; | |
1319 | if (ind <= 19) { | |
1320 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
1321 | } else { | |
1322 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
1323 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
1324 | } | |
1325 | if (C1.w[0] < tmp64) | |
1326 | C1.w[1]++; | |
1327 | // calculate C* and f* | |
1328 | // C* is actually floor(C*) in this case | |
1329 | // C* and f* need shifting and masking, as shown by | |
1330 | // shiftright128[] and maskhigh128[] | |
1331 | // 1 <= x <= 33 | |
1332 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
1333 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
1334 | // the approximation of 10^(-x) was rounded up to 118 bits | |
1335 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
1336 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
1337 | Cstar.w[1] = P256.w[3]; | |
1338 | Cstar.w[0] = P256.w[2]; | |
1339 | fstar.w[3] = 0; | |
1340 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
1341 | fstar.w[1] = P256.w[1]; | |
1342 | fstar.w[0] = P256.w[0]; | |
1343 | } else { // 22 <= ind - 1 <= 33 | |
1344 | Cstar.w[1] = 0; | |
1345 | Cstar.w[0] = P256.w[3]; | |
1346 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
1347 | fstar.w[2] = P256.w[2]; | |
1348 | fstar.w[1] = P256.w[1]; | |
1349 | fstar.w[0] = P256.w[0]; | |
1350 | } | |
1351 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
1352 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
1353 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
1354 | // if floor(C*) is even then C* = floor(C*) - logical right | |
1355 | // shift; C* has p decimal digits, correct by Prop. 1) | |
1356 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
1357 | // shift; C* has p decimal digits, correct by Pr. 1) | |
1358 | // else | |
1359 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
1360 | // correct by Property 1) | |
1361 | // n = C* * 10^(e+x) | |
1362 | ||
1363 | // shift right C* by Ex-128 = shiftright128[ind] | |
1364 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
1365 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
1366 | Cstar.w[0] = | |
1367 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
1368 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
1369 | } else { // 22 <= ind - 1 <= 33 | |
1370 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
1371 | } | |
1372 | // determine inexactness of the rounding of C* | |
1373 | // if (0 < f* - 1/2 < 10^(-x)) then | |
1374 | // the result is exact | |
1375 | // else // if (f* - 1/2 > T*) then | |
1376 | // the result is inexact | |
1377 | if (ind - 1 <= 2) { | |
1378 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
1379 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
1380 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
1381 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
1382 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
1383 | // set the inexact flag |
1384 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
1385 | is_inexact_lt_midpoint = 1; |
1386 | } // else the result is exact | |
1387 | } else { // the result is inexact; f2* <= 1/2 | |
1388 | // set the inexact flag | |
1389 | *pfpsf |= INEXACT_EXCEPTION; | |
1390 | is_inexact_gt_midpoint = 1; | |
1391 | } | |
1392 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
1393 | if (fstar.w[3] > 0x0 || | |
1394 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
1395 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
1396 | (fstar.w[1] || fstar.w[0]))) { | |
1397 | // f2* > 1/2 and the result may be exact | |
1398 | // Calculate f2* - 1/2 | |
1399 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
1400 | tmp64A = fstar.w[3]; | |
1401 | if (tmp64 > fstar.w[2]) | |
1402 | tmp64A--; | |
1403 | if (tmp64A || tmp64 | |
1404 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1405 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1406 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
1407 | // set the inexact flag |
1408 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
1409 | is_inexact_lt_midpoint = 1; |
1410 | } // else the result is exact | |
1411 | } else { // the result is inexact; f2* <= 1/2 | |
1412 | // set the inexact flag | |
1413 | *pfpsf |= INEXACT_EXCEPTION; | |
1414 | is_inexact_gt_midpoint = 1; | |
1415 | } | |
1416 | } else { // if 22 <= ind <= 33 | |
1417 | if (fstar.w[3] > onehalf128[ind - 1] || | |
1418 | (fstar.w[3] == onehalf128[ind - 1] && | |
1419 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
1420 | // f2* > 1/2 and the result may be exact | |
1421 | // Calculate f2* - 1/2 | |
1422 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
1423 | if (tmp64 || fstar.w[2] | |
1424 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1425 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1426 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
1427 | // set the inexact flag |
1428 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
1429 | is_inexact_lt_midpoint = 1; |
1430 | } // else the result is exact | |
1431 | } else { // the result is inexact; f2* <= 1/2 | |
1432 | // set the inexact flag | |
1433 | *pfpsf |= INEXACT_EXCEPTION; | |
1434 | is_inexact_gt_midpoint = 1; | |
1435 | } | |
1436 | } | |
1437 | ||
1438 | // if the result was a midpoint it was rounded away from zero, so | |
1439 | // it will need a correction | |
1440 | // check for midpoints | |
1441 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
1442 | && (fstar.w[1] || fstar.w[0]) | |
1443 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
1444 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1445 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
1446 | // the result is a midpoint; round to nearest | |
1447 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
1448 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
1449 | Cstar.w[0]--; // Cstar.w[0] is now even | |
1450 | is_midpoint_gt_even = 1; | |
1451 | is_inexact_lt_midpoint = 0; | |
1452 | is_inexact_gt_midpoint = 0; | |
1453 | } else { // else MP in [ODD, EVEN] | |
1454 | is_midpoint_lt_even = 1; | |
1455 | is_inexact_lt_midpoint = 0; | |
1456 | is_inexact_gt_midpoint = 0; | |
1457 | } | |
1458 | } | |
1459 | // general correction for RM | |
1460 | if (x_sign && (is_midpoint_gt_even || is_inexact_lt_midpoint)) { | |
1461 | Cstar.w[0] = Cstar.w[0] + 1; | |
1462 | } else if (!x_sign | |
1463 | && (is_midpoint_lt_even || is_inexact_gt_midpoint)) { | |
1464 | Cstar.w[0] = Cstar.w[0] - 1; | |
1465 | } else { | |
1466 | ; // the result is already correct | |
200359e8 | 1467 | } |
b2a00c89 L |
1468 | if (x_sign) |
1469 | res = -Cstar.w[0]; | |
1470 | else | |
1471 | res = Cstar.w[0]; | |
1472 | } else if (exp == 0) { | |
1473 | // 1 <= q <= 10 | |
1474 | // res = +/-C (exact) | |
1475 | if (x_sign) | |
1476 | res = -C1.w[0]; | |
1477 | else | |
1478 | res = C1.w[0]; | |
1479 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
1480 | // res = +/-C * 10^exp (exact) | |
1481 | if (x_sign) | |
1482 | res = -C1.w[0] * ten2k64[exp]; | |
1483 | else | |
1484 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
1485 | } |
1486 | } | |
b2a00c89 L |
1487 | } |
1488 | ||
1489 | BID_RETURN (res); | |
200359e8 L |
1490 | } |
1491 | ||
1492 | /***************************************************************************** | |
1493 | * BID128_to_int32_ceil | |
1494 | ****************************************************************************/ | |
1495 | ||
b2a00c89 | 1496 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_ceil, x) |
200359e8 | 1497 | |
b2a00c89 L |
1498 | int res; |
1499 | UINT64 x_sign; | |
1500 | UINT64 x_exp; | |
1501 | int exp; // unbiased exponent | |
200359e8 | 1502 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
1503 | UINT64 tmp64, tmp64A; |
1504 | BID_UI64DOUBLE tmp1; | |
1505 | unsigned int x_nr_bits; | |
1506 | int q, ind, shift; | |
1507 | UINT128 C1, C; | |
1508 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
1509 | UINT256 fstar; | |
1510 | UINT256 P256; | |
1511 | int is_inexact_lt_midpoint = 0; | |
1512 | int is_inexact_gt_midpoint = 0; | |
1513 | int is_midpoint_lt_even = 0; | |
1514 | int is_midpoint_gt_even = 0; | |
200359e8 L |
1515 | |
1516 | // unpack x | |
b2a00c89 L |
1517 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
1518 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
1519 | C1.w[1] = x.w[1] & MASK_COEFF; | |
1520 | C1.w[0] = x.w[0]; | |
200359e8 L |
1521 | |
1522 | // check for NaN or Infinity | |
b2a00c89 | 1523 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 1524 | // x is special |
b2a00c89 L |
1525 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
1526 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
1527 | // set invalid flag | |
1528 | *pfpsf |= INVALID_EXCEPTION; | |
1529 | // return Integer Indefinite | |
1530 | res = 0x80000000; | |
1531 | } else { // x is QNaN | |
1532 | // set invalid flag | |
1533 | *pfpsf |= INVALID_EXCEPTION; | |
1534 | // return Integer Indefinite | |
1535 | res = 0x80000000; | |
1536 | } | |
1537 | BID_RETURN (res); | |
1538 | } else { // x is not a NaN, so it must be infinity | |
1539 | if (!x_sign) { // x is +inf | |
1540 | // set invalid flag | |
1541 | *pfpsf |= INVALID_EXCEPTION; | |
1542 | // return Integer Indefinite | |
1543 | res = 0x80000000; | |
1544 | } else { // x is -inf | |
1545 | // set invalid flag | |
1546 | *pfpsf |= INVALID_EXCEPTION; | |
1547 | // return Integer Indefinite | |
1548 | res = 0x80000000; | |
200359e8 | 1549 | } |
b2a00c89 L |
1550 | BID_RETURN (res); |
1551 | } | |
1552 | } | |
200359e8 | 1553 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
1554 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
1555 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
1556 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
1557 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
1558 | res = 0x00000000; | |
1559 | BID_RETURN (res); | |
1560 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
1561 | // x is 0 | |
1562 | res = 0x00000000; | |
1563 | BID_RETURN (res); | |
1564 | } else { // x is not special and is not zero | |
1565 | ||
1566 | // q = nr. of decimal digits in x | |
1567 | // determine first the nr. of bits in x | |
1568 | if (C1.w[1] == 0) { | |
1569 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
1570 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
1571 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
1572 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
1573 | x_nr_bits = | |
1574 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1575 | } else { // x < 2^32 | |
1576 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
1577 | x_nr_bits = |
1578 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1579 | } | |
b2a00c89 L |
1580 | } else { // if x < 2^53 |
1581 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 1582 | x_nr_bits = |
b2a00c89 | 1583 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 1584 | } |
b2a00c89 L |
1585 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
1586 | tmp1.d = (double) C1.w[1]; // exact conversion | |
1587 | x_nr_bits = | |
1588 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1589 | } | |
1590 | q = nr_digits[x_nr_bits - 1].digits; | |
1591 | if (q == 0) { | |
1592 | q = nr_digits[x_nr_bits - 1].digits1; | |
1593 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
1594 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
1595 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
1596 | q++; | |
1597 | } | |
1598 | exp = (x_exp >> 49) - 6176; | |
1599 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
1600 | // set invalid flag | |
1601 | *pfpsf |= INVALID_EXCEPTION; | |
1602 | // return Integer Indefinite | |
1603 | res = 0x80000000; | |
1604 | BID_RETURN (res); | |
1605 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
1606 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
1607 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
1608 | // the cases that do not fit are identified here; the ones that fit | |
1609 | // fall through and will be handled with other cases further, | |
1610 | // under '1 <= q + exp <= 10' | |
1611 | if (x_sign) { // if n < 0 and q + exp = 10 | |
1612 | // if n <= -2^31-1 then n is too large | |
1613 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1 | |
1614 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34 | |
1615 | if (q <= 11) { | |
1616 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
1617 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
1618 | if (tmp64 >= 0x50000000aull) { | |
1619 | // set invalid flag | |
1620 | *pfpsf |= INVALID_EXCEPTION; | |
1621 | // return Integer Indefinite | |
1622 | res = 0x80000000; | |
1623 | BID_RETURN (res); | |
1624 | } | |
1625 | // else cases that can be rounded to a 32-bit int fall through | |
1626 | // to '1 <= q + exp <= 10' | |
1627 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
1628 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=> | |
1629 | // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23 | |
1630 | // (scale 2^31+1 up) | |
1631 | tmp64 = 0x50000000aull; | |
1632 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
1633 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
1634 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
1635 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
1636 | } | |
1637 | if (C1.w[1] > C.w[1] | |
1638 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
1639 | // set invalid flag | |
1640 | *pfpsf |= INVALID_EXCEPTION; | |
1641 | // return Integer Indefinite | |
1642 | res = 0x80000000; | |
1643 | BID_RETURN (res); | |
1644 | } | |
1645 | // else cases that can be rounded to a 32-bit int fall through | |
1646 | // to '1 <= q + exp <= 10' | |
1647 | } | |
1648 | } else { // if n > 0 and q + exp = 10 | |
1649 | // if n > 2^31 - 1 then n is too large | |
1650 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1 | |
1651 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=34 | |
1652 | if (q <= 11) { | |
1653 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
1654 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
1655 | if (tmp64 > 0x4fffffff6ull) { | |
1656 | // set invalid flag | |
1657 | *pfpsf |= INVALID_EXCEPTION; | |
1658 | // return Integer Indefinite | |
1659 | res = 0x80000000; | |
1660 | BID_RETURN (res); | |
1661 | } | |
1662 | // else cases that can be rounded to a 32-bit int fall through | |
1663 | // to '1 <= q + exp <= 10' | |
1664 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
1665 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6 <=> | |
1666 | // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23 | |
1667 | // (scale 2^31 up) | |
1668 | tmp64 = 0x4fffffff6ull; | |
1669 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
1670 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
1671 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
1672 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
1673 | } | |
1674 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
1675 | // set invalid flag | |
1676 | *pfpsf |= INVALID_EXCEPTION; | |
1677 | // return Integer Indefinite | |
1678 | res = 0x80000000; | |
1679 | BID_RETURN (res); | |
1680 | } | |
1681 | // else cases that can be rounded to a 32-bit int fall through | |
1682 | // to '1 <= q + exp <= 10' | |
1683 | } | |
200359e8 | 1684 | } |
b2a00c89 L |
1685 | } |
1686 | // n is not too large to be converted to int32: -2^31-1 < n <= 2^31-1 | |
1687 | // Note: some of the cases tested for above fall through to this point | |
1688 | if ((q + exp) <= 0) { | |
1689 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
1690 | // return 0 | |
1691 | if (x_sign) | |
1692 | res = 0x00000000; | |
1693 | else | |
1694 | res = 0x00000001; | |
1695 | BID_RETURN (res); | |
1696 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
1697 | // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded | |
1698 | // toward positive infinity to a 32-bit signed integer | |
1699 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
1700 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
1701 | // chop off ind digits from the lower part of C1 | |
1702 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
1703 | tmp64 = C1.w[0]; | |
1704 | if (ind <= 19) { | |
1705 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
1706 | } else { | |
1707 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
1708 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
1709 | } | |
1710 | if (C1.w[0] < tmp64) | |
1711 | C1.w[1]++; | |
1712 | // calculate C* and f* | |
1713 | // C* is actually floor(C*) in this case | |
1714 | // C* and f* need shifting and masking, as shown by | |
1715 | // shiftright128[] and maskhigh128[] | |
1716 | // 1 <= x <= 33 | |
1717 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
1718 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
1719 | // the approximation of 10^(-x) was rounded up to 118 bits | |
1720 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
1721 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
1722 | Cstar.w[1] = P256.w[3]; | |
1723 | Cstar.w[0] = P256.w[2]; | |
1724 | fstar.w[3] = 0; | |
1725 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
1726 | fstar.w[1] = P256.w[1]; | |
1727 | fstar.w[0] = P256.w[0]; | |
1728 | } else { // 22 <= ind - 1 <= 33 | |
1729 | Cstar.w[1] = 0; | |
1730 | Cstar.w[0] = P256.w[3]; | |
1731 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
1732 | fstar.w[2] = P256.w[2]; | |
1733 | fstar.w[1] = P256.w[1]; | |
1734 | fstar.w[0] = P256.w[0]; | |
1735 | } | |
1736 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
1737 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
1738 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
1739 | // if floor(C*) is even then C* = floor(C*) - logical right | |
1740 | // shift; C* has p decimal digits, correct by Prop. 1) | |
1741 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
1742 | // shift; C* has p decimal digits, correct by Pr. 1) | |
1743 | // else | |
1744 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
1745 | // correct by Property 1) | |
1746 | // n = C* * 10^(e+x) | |
1747 | ||
1748 | // shift right C* by Ex-128 = shiftright128[ind] | |
1749 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
1750 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
1751 | Cstar.w[0] = | |
1752 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
1753 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
1754 | } else { // 22 <= ind - 1 <= 33 | |
1755 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
1756 | } | |
1757 | // determine inexactness of the rounding of C* | |
1758 | // if (0 < f* - 1/2 < 10^(-x)) then | |
1759 | // the result is exact | |
1760 | // else // if (f* - 1/2 > T*) then | |
1761 | // the result is inexact | |
1762 | if (ind - 1 <= 2) { | |
1763 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
1764 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
1765 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
1766 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
1767 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
1768 | is_inexact_lt_midpoint = 1; | |
1769 | } // else the result is exact | |
1770 | } else { // the result is inexact; f2* <= 1/2 | |
1771 | is_inexact_gt_midpoint = 1; | |
1772 | } | |
1773 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
1774 | if (fstar.w[3] > 0x0 || | |
1775 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
1776 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
1777 | (fstar.w[1] || fstar.w[0]))) { | |
1778 | // f2* > 1/2 and the result may be exact | |
1779 | // Calculate f2* - 1/2 | |
1780 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
1781 | tmp64A = fstar.w[3]; | |
1782 | if (tmp64 > fstar.w[2]) | |
1783 | tmp64A--; | |
1784 | if (tmp64A || tmp64 | |
1785 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1786 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1787 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
1788 | is_inexact_lt_midpoint = 1; | |
1789 | } // else the result is exact | |
1790 | } else { // the result is inexact; f2* <= 1/2 | |
1791 | is_inexact_gt_midpoint = 1; | |
1792 | } | |
1793 | } else { // if 22 <= ind <= 33 | |
1794 | if (fstar.w[3] > onehalf128[ind - 1] || | |
1795 | (fstar.w[3] == onehalf128[ind - 1] && | |
1796 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
1797 | // f2* > 1/2 and the result may be exact | |
1798 | // Calculate f2* - 1/2 | |
1799 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
1800 | if (tmp64 || fstar.w[2] | |
1801 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
1802 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1803 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
1804 | is_inexact_lt_midpoint = 1; | |
1805 | } // else the result is exact | |
1806 | } else { // the result is inexact; f2* <= 1/2 | |
1807 | is_inexact_gt_midpoint = 1; | |
200359e8 L |
1808 | } |
1809 | } | |
b2a00c89 L |
1810 | |
1811 | // if the result was a midpoint it was rounded away from zero, so | |
1812 | // it will need a correction | |
1813 | // check for midpoints | |
1814 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
1815 | && (fstar.w[1] || fstar.w[0]) | |
1816 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
1817 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
1818 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
1819 | // the result is a midpoint; round to nearest | |
1820 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
1821 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
1822 | Cstar.w[0]--; // Cstar.w[0] is now even | |
1823 | is_midpoint_gt_even = 1; | |
1824 | is_inexact_lt_midpoint = 0; | |
1825 | is_inexact_gt_midpoint = 0; | |
1826 | } else { // else MP in [ODD, EVEN] | |
1827 | is_midpoint_lt_even = 1; | |
1828 | is_inexact_lt_midpoint = 0; | |
1829 | is_inexact_gt_midpoint = 0; | |
1830 | } | |
1831 | } | |
1832 | // general correction for RM | |
1833 | if (x_sign && (is_midpoint_lt_even || is_inexact_gt_midpoint)) { | |
1834 | Cstar.w[0] = Cstar.w[0] - 1; | |
1835 | } else if (!x_sign | |
1836 | && (is_midpoint_gt_even || is_inexact_lt_midpoint)) { | |
1837 | Cstar.w[0] = Cstar.w[0] + 1; | |
1838 | } else { | |
1839 | ; // the result is already correct | |
1840 | } | |
200359e8 | 1841 | if (x_sign) |
b2a00c89 | 1842 | res = -Cstar.w[0]; |
200359e8 | 1843 | else |
b2a00c89 L |
1844 | res = Cstar.w[0]; |
1845 | } else if (exp == 0) { | |
1846 | // 1 <= q <= 10 | |
1847 | // res = +/-C (exact) | |
1848 | if (x_sign) | |
1849 | res = -C1.w[0]; | |
1850 | else | |
1851 | res = C1.w[0]; | |
1852 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
1853 | // res = +/-C * 10^exp (exact) | |
1854 | if (x_sign) | |
1855 | res = -C1.w[0] * ten2k64[exp]; | |
1856 | else | |
1857 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
1858 | } |
1859 | } | |
b2a00c89 L |
1860 | } |
1861 | ||
1862 | BID_RETURN (res); | |
200359e8 L |
1863 | } |
1864 | ||
1865 | /***************************************************************************** | |
1866 | * BID128_to_int32_xceil | |
1867 | ****************************************************************************/ | |
1868 | ||
b2a00c89 | 1869 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xceil, x) |
200359e8 | 1870 | |
b2a00c89 L |
1871 | int res; |
1872 | UINT64 x_sign; | |
1873 | UINT64 x_exp; | |
1874 | int exp; // unbiased exponent | |
200359e8 | 1875 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
1876 | UINT64 tmp64, tmp64A; |
1877 | BID_UI64DOUBLE tmp1; | |
1878 | unsigned int x_nr_bits; | |
1879 | int q, ind, shift; | |
1880 | UINT128 C1, C; | |
1881 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
1882 | UINT256 fstar; | |
1883 | UINT256 P256; | |
1884 | int is_inexact_lt_midpoint = 0; | |
1885 | int is_inexact_gt_midpoint = 0; | |
1886 | int is_midpoint_lt_even = 0; | |
1887 | int is_midpoint_gt_even = 0; | |
200359e8 L |
1888 | |
1889 | // unpack x | |
b2a00c89 L |
1890 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
1891 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
1892 | C1.w[1] = x.w[1] & MASK_COEFF; | |
1893 | C1.w[0] = x.w[0]; | |
200359e8 L |
1894 | |
1895 | // check for NaN or Infinity | |
b2a00c89 | 1896 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 1897 | // x is special |
b2a00c89 L |
1898 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
1899 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
1900 | // set invalid flag | |
1901 | *pfpsf |= INVALID_EXCEPTION; | |
1902 | // return Integer Indefinite | |
1903 | res = 0x80000000; | |
1904 | } else { // x is QNaN | |
1905 | // set invalid flag | |
1906 | *pfpsf |= INVALID_EXCEPTION; | |
1907 | // return Integer Indefinite | |
1908 | res = 0x80000000; | |
200359e8 | 1909 | } |
b2a00c89 L |
1910 | BID_RETURN (res); |
1911 | } else { // x is not a NaN, so it must be infinity | |
1912 | if (!x_sign) { // x is +inf | |
1913 | // set invalid flag | |
1914 | *pfpsf |= INVALID_EXCEPTION; | |
1915 | // return Integer Indefinite | |
1916 | res = 0x80000000; | |
1917 | } else { // x is -inf | |
1918 | // set invalid flag | |
1919 | *pfpsf |= INVALID_EXCEPTION; | |
1920 | // return Integer Indefinite | |
1921 | res = 0x80000000; | |
1922 | } | |
1923 | BID_RETURN (res); | |
1924 | } | |
1925 | } | |
200359e8 | 1926 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
1927 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
1928 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
1929 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
1930 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
1931 | res = 0x00000000; | |
1932 | BID_RETURN (res); | |
1933 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
1934 | // x is 0 | |
1935 | res = 0x00000000; | |
1936 | BID_RETURN (res); | |
1937 | } else { // x is not special and is not zero | |
1938 | ||
1939 | // q = nr. of decimal digits in x | |
1940 | // determine first the nr. of bits in x | |
1941 | if (C1.w[1] == 0) { | |
1942 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
1943 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
1944 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
1945 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
1946 | x_nr_bits = | |
1947 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1948 | } else { // x < 2^32 | |
1949 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
1950 | x_nr_bits = |
1951 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1952 | } | |
b2a00c89 L |
1953 | } else { // if x < 2^53 |
1954 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 1955 | x_nr_bits = |
b2a00c89 | 1956 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 1957 | } |
b2a00c89 L |
1958 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
1959 | tmp1.d = (double) C1.w[1]; // exact conversion | |
1960 | x_nr_bits = | |
1961 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1962 | } | |
1963 | q = nr_digits[x_nr_bits - 1].digits; | |
1964 | if (q == 0) { | |
1965 | q = nr_digits[x_nr_bits - 1].digits1; | |
1966 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
1967 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
1968 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
1969 | q++; | |
1970 | } | |
1971 | exp = (x_exp >> 49) - 6176; | |
1972 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
1973 | // set invalid flag | |
1974 | *pfpsf |= INVALID_EXCEPTION; | |
1975 | // return Integer Indefinite | |
1976 | res = 0x80000000; | |
1977 | BID_RETURN (res); | |
1978 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
1979 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
1980 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
1981 | // the cases that do not fit are identified here; the ones that fit | |
1982 | // fall through and will be handled with other cases further, | |
1983 | // under '1 <= q + exp <= 10' | |
1984 | if (x_sign) { // if n < 0 and q + exp = 10 | |
1985 | // if n <= -2^31-1 then n is too large | |
1986 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1 | |
1987 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34 | |
1988 | if (q <= 11) { | |
1989 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
1990 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
1991 | if (tmp64 >= 0x50000000aull) { | |
1992 | // set invalid flag | |
1993 | *pfpsf |= INVALID_EXCEPTION; | |
1994 | // return Integer Indefinite | |
1995 | res = 0x80000000; | |
1996 | BID_RETURN (res); | |
1997 | } | |
1998 | // else cases that can be rounded to a 32-bit int fall through | |
1999 | // to '1 <= q + exp <= 10' | |
2000 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2001 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=> | |
2002 | // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23 | |
2003 | // (scale 2^31+1 up) | |
2004 | tmp64 = 0x50000000aull; | |
2005 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2006 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2007 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2008 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2009 | } | |
2010 | if (C1.w[1] > C.w[1] | |
2011 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
2012 | // set invalid flag | |
2013 | *pfpsf |= INVALID_EXCEPTION; | |
2014 | // return Integer Indefinite | |
2015 | res = 0x80000000; | |
2016 | BID_RETURN (res); | |
2017 | } | |
2018 | // else cases that can be rounded to a 32-bit int fall through | |
2019 | // to '1 <= q + exp <= 10' | |
2020 | } | |
2021 | } else { // if n > 0 and q + exp = 10 | |
2022 | // if n > 2^31 - 1 then n is too large | |
2023 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1 | |
2024 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=34 | |
2025 | if (q <= 11) { | |
2026 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
2027 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
2028 | if (tmp64 > 0x4fffffff6ull) { | |
2029 | // set invalid flag | |
2030 | *pfpsf |= INVALID_EXCEPTION; | |
2031 | // return Integer Indefinite | |
2032 | res = 0x80000000; | |
2033 | BID_RETURN (res); | |
2034 | } | |
2035 | // else cases that can be rounded to a 32-bit int fall through | |
2036 | // to '1 <= q + exp <= 10' | |
2037 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2038 | // 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6 <=> | |
2039 | // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23 | |
2040 | // (scale 2^31 up) | |
2041 | tmp64 = 0x4fffffff6ull; | |
2042 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2043 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2044 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2045 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2046 | } | |
2047 | if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) { | |
2048 | // set invalid flag | |
2049 | *pfpsf |= INVALID_EXCEPTION; | |
2050 | // return Integer Indefinite | |
2051 | res = 0x80000000; | |
2052 | BID_RETURN (res); | |
2053 | } | |
2054 | // else cases that can be rounded to a 32-bit int fall through | |
2055 | // to '1 <= q + exp <= 10' | |
200359e8 L |
2056 | } |
2057 | } | |
b2a00c89 L |
2058 | } |
2059 | // n is not too large to be converted to int32: -2^31-1 < n <= 2^31-1 | |
2060 | // Note: some of the cases tested for above fall through to this point | |
2061 | if ((q + exp) <= 0) { | |
2062 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
2063 | // set inexact flag | |
2064 | *pfpsf |= INEXACT_EXCEPTION; | |
2065 | // return 0 | |
2066 | if (x_sign) | |
2067 | res = 0x00000000; | |
2068 | else | |
2069 | res = 0x00000001; | |
2070 | BID_RETURN (res); | |
2071 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
2072 | // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded | |
2073 | // toward positive infinity to a 32-bit signed integer | |
2074 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
2075 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
2076 | // chop off ind digits from the lower part of C1 | |
2077 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
2078 | tmp64 = C1.w[0]; | |
2079 | if (ind <= 19) { | |
2080 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
2081 | } else { | |
2082 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
2083 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
2084 | } | |
2085 | if (C1.w[0] < tmp64) | |
2086 | C1.w[1]++; | |
2087 | // calculate C* and f* | |
2088 | // C* is actually floor(C*) in this case | |
2089 | // C* and f* need shifting and masking, as shown by | |
2090 | // shiftright128[] and maskhigh128[] | |
2091 | // 1 <= x <= 33 | |
2092 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
2093 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
2094 | // the approximation of 10^(-x) was rounded up to 118 bits | |
2095 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
2096 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2097 | Cstar.w[1] = P256.w[3]; | |
2098 | Cstar.w[0] = P256.w[2]; | |
2099 | fstar.w[3] = 0; | |
2100 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
2101 | fstar.w[1] = P256.w[1]; | |
2102 | fstar.w[0] = P256.w[0]; | |
2103 | } else { // 22 <= ind - 1 <= 33 | |
2104 | Cstar.w[1] = 0; | |
2105 | Cstar.w[0] = P256.w[3]; | |
2106 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
2107 | fstar.w[2] = P256.w[2]; | |
2108 | fstar.w[1] = P256.w[1]; | |
2109 | fstar.w[0] = P256.w[0]; | |
2110 | } | |
2111 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
2112 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
2113 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
2114 | // if floor(C*) is even then C* = floor(C*) - logical right | |
2115 | // shift; C* has p decimal digits, correct by Prop. 1) | |
2116 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
2117 | // shift; C* has p decimal digits, correct by Pr. 1) | |
2118 | // else | |
2119 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
2120 | // correct by Property 1) | |
2121 | // n = C* * 10^(e+x) | |
2122 | ||
2123 | // shift right C* by Ex-128 = shiftright128[ind] | |
2124 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
2125 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2126 | Cstar.w[0] = | |
2127 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
2128 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
2129 | } else { // 22 <= ind - 1 <= 33 | |
2130 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
2131 | } | |
2132 | // determine inexactness of the rounding of C* | |
2133 | // if (0 < f* - 1/2 < 10^(-x)) then | |
2134 | // the result is exact | |
2135 | // else // if (f* - 1/2 > T*) then | |
2136 | // the result is inexact | |
2137 | if (ind - 1 <= 2) { | |
2138 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
2139 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
2140 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
2141 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
2142 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2143 | // set the inexact flag |
2144 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2145 | is_inexact_lt_midpoint = 1; |
2146 | } // else the result is exact | |
2147 | } else { // the result is inexact; f2* <= 1/2 | |
2148 | // set the inexact flag | |
2149 | *pfpsf |= INEXACT_EXCEPTION; | |
2150 | is_inexact_gt_midpoint = 1; | |
2151 | } | |
2152 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
2153 | if (fstar.w[3] > 0x0 || | |
2154 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
2155 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
2156 | (fstar.w[1] || fstar.w[0]))) { | |
2157 | // f2* > 1/2 and the result may be exact | |
2158 | // Calculate f2* - 1/2 | |
2159 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
2160 | tmp64A = fstar.w[3]; | |
2161 | if (tmp64 > fstar.w[2]) | |
2162 | tmp64A--; | |
2163 | if (tmp64A || tmp64 | |
2164 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2165 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2166 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2167 | // set the inexact flag |
2168 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2169 | is_inexact_lt_midpoint = 1; |
2170 | } // else the result is exact | |
2171 | } else { // the result is inexact; f2* <= 1/2 | |
2172 | // set the inexact flag | |
2173 | *pfpsf |= INEXACT_EXCEPTION; | |
2174 | is_inexact_gt_midpoint = 1; | |
2175 | } | |
2176 | } else { // if 22 <= ind <= 33 | |
2177 | if (fstar.w[3] > onehalf128[ind - 1] || | |
2178 | (fstar.w[3] == onehalf128[ind - 1] && | |
2179 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
2180 | // f2* > 1/2 and the result may be exact | |
2181 | // Calculate f2* - 1/2 | |
2182 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
2183 | if (tmp64 || fstar.w[2] | |
2184 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2185 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2186 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2187 | // set the inexact flag |
2188 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2189 | is_inexact_lt_midpoint = 1; |
2190 | } // else the result is exact | |
2191 | } else { // the result is inexact; f2* <= 1/2 | |
2192 | // set the inexact flag | |
2193 | *pfpsf |= INEXACT_EXCEPTION; | |
2194 | is_inexact_gt_midpoint = 1; | |
2195 | } | |
2196 | } | |
2197 | ||
2198 | // if the result was a midpoint it was rounded away from zero, so | |
2199 | // it will need a correction | |
2200 | // check for midpoints | |
2201 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
2202 | && (fstar.w[1] || fstar.w[0]) | |
2203 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
2204 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2205 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
2206 | // the result is a midpoint; round to nearest | |
2207 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
2208 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
2209 | Cstar.w[0]--; // Cstar.w[0] is now even | |
2210 | is_midpoint_gt_even = 1; | |
2211 | is_inexact_lt_midpoint = 0; | |
2212 | is_inexact_gt_midpoint = 0; | |
2213 | } else { // else MP in [ODD, EVEN] | |
2214 | is_midpoint_lt_even = 1; | |
2215 | is_inexact_lt_midpoint = 0; | |
2216 | is_inexact_gt_midpoint = 0; | |
2217 | } | |
200359e8 | 2218 | } |
b2a00c89 L |
2219 | // general correction for RM |
2220 | if (x_sign && (is_midpoint_lt_even || is_inexact_gt_midpoint)) { | |
2221 | Cstar.w[0] = Cstar.w[0] - 1; | |
2222 | } else if (!x_sign | |
2223 | && (is_midpoint_gt_even || is_inexact_lt_midpoint)) { | |
2224 | Cstar.w[0] = Cstar.w[0] + 1; | |
2225 | } else { | |
2226 | ; // the result is already correct | |
2227 | } | |
2228 | if (x_sign) | |
2229 | res = -Cstar.w[0]; | |
2230 | else | |
2231 | res = Cstar.w[0]; | |
2232 | } else if (exp == 0) { | |
2233 | // 1 <= q <= 10 | |
2234 | // res = +/-C (exact) | |
2235 | if (x_sign) | |
2236 | res = -C1.w[0]; | |
2237 | else | |
2238 | res = C1.w[0]; | |
2239 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
2240 | // res = +/-C * 10^exp (exact) | |
2241 | if (x_sign) | |
2242 | res = -C1.w[0] * ten2k64[exp]; | |
2243 | else | |
2244 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
2245 | } |
2246 | } | |
b2a00c89 L |
2247 | } |
2248 | ||
2249 | BID_RETURN (res); | |
200359e8 L |
2250 | } |
2251 | ||
2252 | /***************************************************************************** | |
2253 | * BID128_to_int32_int | |
2254 | ****************************************************************************/ | |
2255 | ||
b2a00c89 | 2256 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_int, x) |
200359e8 | 2257 | |
b2a00c89 L |
2258 | int res; |
2259 | UINT64 x_sign; | |
2260 | UINT64 x_exp; | |
2261 | int exp; // unbiased exponent | |
200359e8 | 2262 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
2263 | UINT64 tmp64, tmp64A; |
2264 | BID_UI64DOUBLE tmp1; | |
2265 | unsigned int x_nr_bits; | |
2266 | int q, ind, shift; | |
2267 | UINT128 C1, C; | |
2268 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
2269 | UINT256 fstar; | |
2270 | UINT256 P256; | |
2271 | int is_inexact_gt_midpoint = 0; | |
2272 | int is_midpoint_lt_even = 0; | |
200359e8 L |
2273 | |
2274 | // unpack x | |
b2a00c89 L |
2275 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
2276 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
2277 | C1.w[1] = x.w[1] & MASK_COEFF; | |
2278 | C1.w[0] = x.w[0]; | |
200359e8 L |
2279 | |
2280 | // check for NaN or Infinity | |
b2a00c89 | 2281 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 2282 | // x is special |
b2a00c89 L |
2283 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
2284 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
2285 | // set invalid flag | |
2286 | *pfpsf |= INVALID_EXCEPTION; | |
2287 | // return Integer Indefinite | |
2288 | res = 0x80000000; | |
2289 | } else { // x is QNaN | |
2290 | // set invalid flag | |
2291 | *pfpsf |= INVALID_EXCEPTION; | |
2292 | // return Integer Indefinite | |
2293 | res = 0x80000000; | |
2294 | } | |
2295 | BID_RETURN (res); | |
2296 | } else { // x is not a NaN, so it must be infinity | |
2297 | if (!x_sign) { // x is +inf | |
2298 | // set invalid flag | |
2299 | *pfpsf |= INVALID_EXCEPTION; | |
2300 | // return Integer Indefinite | |
2301 | res = 0x80000000; | |
2302 | } else { // x is -inf | |
2303 | // set invalid flag | |
2304 | *pfpsf |= INVALID_EXCEPTION; | |
2305 | // return Integer Indefinite | |
2306 | res = 0x80000000; | |
200359e8 | 2307 | } |
b2a00c89 L |
2308 | BID_RETURN (res); |
2309 | } | |
2310 | } | |
200359e8 | 2311 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
2312 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
2313 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
2314 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
2315 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
2316 | res = 0x00000000; | |
2317 | BID_RETURN (res); | |
2318 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
2319 | // x is 0 | |
2320 | res = 0x00000000; | |
2321 | BID_RETURN (res); | |
2322 | } else { // x is not special and is not zero | |
2323 | ||
2324 | // q = nr. of decimal digits in x | |
2325 | // determine first the nr. of bits in x | |
2326 | if (C1.w[1] == 0) { | |
2327 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
2328 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
2329 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
2330 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
2331 | x_nr_bits = | |
2332 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2333 | } else { // x < 2^32 | |
2334 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
2335 | x_nr_bits = |
2336 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2337 | } | |
b2a00c89 L |
2338 | } else { // if x < 2^53 |
2339 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 2340 | x_nr_bits = |
b2a00c89 | 2341 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 2342 | } |
b2a00c89 L |
2343 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
2344 | tmp1.d = (double) C1.w[1]; // exact conversion | |
2345 | x_nr_bits = | |
2346 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2347 | } | |
2348 | q = nr_digits[x_nr_bits - 1].digits; | |
2349 | if (q == 0) { | |
2350 | q = nr_digits[x_nr_bits - 1].digits1; | |
2351 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
2352 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
2353 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
2354 | q++; | |
2355 | } | |
2356 | exp = (x_exp >> 49) - 6176; | |
2357 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
2358 | // set invalid flag | |
2359 | *pfpsf |= INVALID_EXCEPTION; | |
2360 | // return Integer Indefinite | |
2361 | res = 0x80000000; | |
2362 | BID_RETURN (res); | |
2363 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
2364 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
2365 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
2366 | // the cases that do not fit are identified here; the ones that fit | |
2367 | // fall through and will be handled with other cases further, | |
2368 | // under '1 <= q + exp <= 10' | |
2369 | if (x_sign) { // if n < 0 and q + exp = 10 | |
2370 | // if n <= -2^31 - 1 then n is too large | |
2371 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1 | |
2372 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34 | |
2373 | if (q <= 11) { | |
2374 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
2375 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
2376 | if (tmp64 >= 0x50000000aull) { | |
2377 | // set invalid flag | |
2378 | *pfpsf |= INVALID_EXCEPTION; | |
2379 | // return Integer Indefinite | |
2380 | res = 0x80000000; | |
2381 | BID_RETURN (res); | |
2382 | } | |
2383 | // else cases that can be rounded to a 32-bit int fall through | |
2384 | // to '1 <= q + exp <= 10' | |
2385 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2386 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=> | |
2387 | // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23 | |
2388 | // (scale 2^31+1 up) | |
2389 | tmp64 = 0x50000000aull; | |
2390 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2391 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2392 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2393 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2394 | } | |
2395 | if (C1.w[1] > C.w[1] | |
2396 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
2397 | // set invalid flag | |
2398 | *pfpsf |= INVALID_EXCEPTION; | |
2399 | // return Integer Indefinite | |
2400 | res = 0x80000000; | |
2401 | BID_RETURN (res); | |
2402 | } | |
2403 | // else cases that can be rounded to a 32-bit int fall through | |
2404 | // to '1 <= q + exp <= 10' | |
2405 | } | |
2406 | } else { // if n > 0 and q + exp = 10 | |
2407 | // if n >= 2^31 then n is too large | |
2408 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31 | |
2409 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34 | |
2410 | if (q <= 11) { | |
2411 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
2412 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
2413 | if (tmp64 >= 0x500000000ull) { | |
2414 | // set invalid flag | |
2415 | *pfpsf |= INVALID_EXCEPTION; | |
2416 | // return Integer Indefinite | |
2417 | res = 0x80000000; | |
2418 | BID_RETURN (res); | |
2419 | } | |
2420 | // else cases that can be rounded to a 32-bit int fall through | |
2421 | // to '1 <= q + exp <= 10' | |
2422 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2423 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=> | |
2424 | // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
2425 | // (scale 2^31-1/2 up) | |
2426 | tmp64 = 0x500000000ull; | |
2427 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2428 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2429 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2430 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2431 | } | |
2432 | if (C1.w[1] > C.w[1] | |
2433 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
2434 | // set invalid flag | |
2435 | *pfpsf |= INVALID_EXCEPTION; | |
2436 | // return Integer Indefinite | |
2437 | res = 0x80000000; | |
2438 | BID_RETURN (res); | |
2439 | } | |
2440 | // else cases that can be rounded to a 32-bit int fall through | |
2441 | // to '1 <= q + exp <= 10' | |
2442 | } | |
200359e8 | 2443 | } |
b2a00c89 L |
2444 | } |
2445 | // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31 | |
2446 | // Note: some of the cases tested for above fall through to this point | |
2447 | if ((q + exp) <= 0) { | |
2448 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
2449 | // return 0 | |
2450 | res = 0x00000000; | |
2451 | BID_RETURN (res); | |
2452 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
2453 | // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded | |
2454 | // toward zero to a 32-bit signed integer | |
2455 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
2456 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
2457 | // chop off ind digits from the lower part of C1 | |
2458 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
2459 | tmp64 = C1.w[0]; | |
2460 | if (ind <= 19) { | |
2461 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
2462 | } else { | |
2463 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
2464 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
2465 | } | |
2466 | if (C1.w[0] < tmp64) | |
2467 | C1.w[1]++; | |
2468 | // calculate C* and f* | |
2469 | // C* is actually floor(C*) in this case | |
2470 | // C* and f* need shifting and masking, as shown by | |
2471 | // shiftright128[] and maskhigh128[] | |
2472 | // 1 <= x <= 33 | |
2473 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
2474 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
2475 | // the approximation of 10^(-x) was rounded up to 118 bits | |
2476 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
2477 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2478 | Cstar.w[1] = P256.w[3]; | |
2479 | Cstar.w[0] = P256.w[2]; | |
2480 | fstar.w[3] = 0; | |
2481 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
2482 | fstar.w[1] = P256.w[1]; | |
2483 | fstar.w[0] = P256.w[0]; | |
2484 | } else { // 22 <= ind - 1 <= 33 | |
2485 | Cstar.w[1] = 0; | |
2486 | Cstar.w[0] = P256.w[3]; | |
2487 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
2488 | fstar.w[2] = P256.w[2]; | |
2489 | fstar.w[1] = P256.w[1]; | |
2490 | fstar.w[0] = P256.w[0]; | |
2491 | } | |
2492 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
2493 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
2494 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
2495 | // if floor(C*) is even then C* = floor(C*) - logical right | |
2496 | // shift; C* has p decimal digits, correct by Prop. 1) | |
2497 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
2498 | // shift; C* has p decimal digits, correct by Pr. 1) | |
2499 | // else | |
2500 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
2501 | // correct by Property 1) | |
2502 | // n = C* * 10^(e+x) | |
2503 | ||
2504 | // shift right C* by Ex-128 = shiftright128[ind] | |
2505 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
2506 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2507 | Cstar.w[0] = | |
2508 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
2509 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
2510 | } else { // 22 <= ind - 1 <= 33 | |
2511 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
2512 | } | |
2513 | // determine inexactness of the rounding of C* | |
2514 | // if (0 < f* - 1/2 < 10^(-x)) then | |
2515 | // the result is exact | |
2516 | // else // if (f* - 1/2 > T*) then | |
2517 | // the result is inexact | |
2518 | if (ind - 1 <= 2) { | |
2519 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
2520 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
2521 | if ((tmp64 > ten2mk128trunc[ind - 1].w[1] | |
2522 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
2523 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0]))) { | |
2524 | } // else the result is exact | |
2525 | } else { // the result is inexact; f2* <= 1/2 | |
2526 | is_inexact_gt_midpoint = 1; | |
2527 | } | |
2528 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
2529 | if (fstar.w[3] > 0x0 || | |
2530 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
2531 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
2532 | (fstar.w[1] || fstar.w[0]))) { | |
2533 | // f2* > 1/2 and the result may be exact | |
2534 | // Calculate f2* - 1/2 | |
2535 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
2536 | tmp64A = fstar.w[3]; | |
2537 | if (tmp64 > fstar.w[2]) | |
2538 | tmp64A--; | |
2539 | if (tmp64A || tmp64 | |
2540 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2541 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2542 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
2543 | } // else the result is exact | |
2544 | } else { // the result is inexact; f2* <= 1/2 | |
2545 | is_inexact_gt_midpoint = 1; | |
2546 | } | |
2547 | } else { // if 22 <= ind <= 33 | |
2548 | if (fstar.w[3] > onehalf128[ind - 1] || | |
2549 | (fstar.w[3] == onehalf128[ind - 1] && | |
2550 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
2551 | // f2* > 1/2 and the result may be exact | |
2552 | // Calculate f2* - 1/2 | |
2553 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
2554 | if (tmp64 || fstar.w[2] | |
2555 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2556 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2557 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
2558 | } // else the result is exact | |
2559 | } else { // the result is inexact; f2* <= 1/2 | |
2560 | is_inexact_gt_midpoint = 1; | |
200359e8 L |
2561 | } |
2562 | } | |
b2a00c89 L |
2563 | |
2564 | // if the result was a midpoint it was rounded away from zero, so | |
2565 | // it will need a correction | |
2566 | // check for midpoints | |
2567 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) && | |
2568 | (fstar.w[1] || fstar.w[0]) && | |
2569 | (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] || | |
2570 | (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] && | |
2571 | fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
2572 | // the result is a midpoint; round to nearest | |
2573 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
2574 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
2575 | Cstar.w[0]--; // Cstar.w[0] is now even | |
2576 | is_inexact_gt_midpoint = 0; | |
2577 | } else { // else MP in [ODD, EVEN] | |
2578 | is_midpoint_lt_even = 1; | |
2579 | is_inexact_gt_midpoint = 0; | |
2580 | } | |
200359e8 | 2581 | } |
b2a00c89 L |
2582 | // general correction for RZ |
2583 | if (is_midpoint_lt_even || is_inexact_gt_midpoint) { | |
2584 | Cstar.w[0] = Cstar.w[0] - 1; | |
2585 | } else { | |
2586 | ; // exact, the result is already correct | |
2587 | } | |
2588 | if (x_sign) | |
2589 | res = -Cstar.w[0]; | |
2590 | else | |
2591 | res = Cstar.w[0]; | |
2592 | } else if (exp == 0) { | |
2593 | // 1 <= q <= 10 | |
2594 | // res = +/-C (exact) | |
2595 | if (x_sign) | |
2596 | res = -C1.w[0]; | |
2597 | else | |
2598 | res = C1.w[0]; | |
2599 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
2600 | // res = +/-C * 10^exp (exact) | |
2601 | if (x_sign) | |
2602 | res = -C1.w[0] * ten2k64[exp]; | |
2603 | else | |
2604 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
2605 | } |
2606 | } | |
b2a00c89 L |
2607 | } |
2608 | ||
2609 | BID_RETURN (res); | |
200359e8 L |
2610 | } |
2611 | ||
2612 | /***************************************************************************** | |
2613 | * BID128_to_int32_xint | |
2614 | ****************************************************************************/ | |
2615 | ||
b2a00c89 | 2616 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xint, x) |
200359e8 | 2617 | |
b2a00c89 L |
2618 | int res; |
2619 | UINT64 x_sign; | |
2620 | UINT64 x_exp; | |
2621 | int exp; // unbiased exponent | |
200359e8 | 2622 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
2623 | UINT64 tmp64, tmp64A; |
2624 | BID_UI64DOUBLE tmp1; | |
2625 | unsigned int x_nr_bits; | |
2626 | int q, ind, shift; | |
2627 | UINT128 C1, C; | |
2628 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
2629 | UINT256 fstar; | |
2630 | UINT256 P256; | |
2631 | int is_inexact_gt_midpoint = 0; | |
2632 | int is_midpoint_lt_even = 0; | |
200359e8 L |
2633 | |
2634 | // unpack x | |
b2a00c89 L |
2635 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
2636 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
2637 | C1.w[1] = x.w[1] & MASK_COEFF; | |
2638 | C1.w[0] = x.w[0]; | |
200359e8 L |
2639 | |
2640 | // check for NaN or Infinity | |
b2a00c89 | 2641 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 2642 | // x is special |
b2a00c89 L |
2643 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
2644 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
2645 | // set invalid flag | |
2646 | *pfpsf |= INVALID_EXCEPTION; | |
2647 | // return Integer Indefinite | |
2648 | res = 0x80000000; | |
2649 | } else { // x is QNaN | |
2650 | // set invalid flag | |
2651 | *pfpsf |= INVALID_EXCEPTION; | |
2652 | // return Integer Indefinite | |
2653 | res = 0x80000000; | |
2654 | } | |
2655 | BID_RETURN (res); | |
2656 | } else { // x is not a NaN, so it must be infinity | |
2657 | if (!x_sign) { // x is +inf | |
2658 | // set invalid flag | |
2659 | *pfpsf |= INVALID_EXCEPTION; | |
2660 | // return Integer Indefinite | |
2661 | res = 0x80000000; | |
2662 | } else { // x is -inf | |
2663 | // set invalid flag | |
2664 | *pfpsf |= INVALID_EXCEPTION; | |
2665 | // return Integer Indefinite | |
2666 | res = 0x80000000; | |
200359e8 | 2667 | } |
b2a00c89 L |
2668 | BID_RETURN (res); |
2669 | } | |
2670 | } | |
200359e8 | 2671 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
2672 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
2673 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
2674 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
2675 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
2676 | res = 0x00000000; | |
2677 | BID_RETURN (res); | |
2678 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
2679 | // x is 0 | |
2680 | res = 0x00000000; | |
2681 | BID_RETURN (res); | |
2682 | } else { // x is not special and is not zero | |
2683 | ||
2684 | // q = nr. of decimal digits in x | |
2685 | // determine first the nr. of bits in x | |
2686 | if (C1.w[1] == 0) { | |
2687 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
2688 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
2689 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
2690 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
2691 | x_nr_bits = | |
2692 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2693 | } else { // x < 2^32 | |
2694 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
2695 | x_nr_bits = |
2696 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2697 | } | |
b2a00c89 L |
2698 | } else { // if x < 2^53 |
2699 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 2700 | x_nr_bits = |
b2a00c89 | 2701 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 2702 | } |
b2a00c89 L |
2703 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
2704 | tmp1.d = (double) C1.w[1]; // exact conversion | |
2705 | x_nr_bits = | |
2706 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2707 | } | |
2708 | q = nr_digits[x_nr_bits - 1].digits; | |
2709 | if (q == 0) { | |
2710 | q = nr_digits[x_nr_bits - 1].digits1; | |
2711 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
2712 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
2713 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
2714 | q++; | |
2715 | } | |
2716 | exp = (x_exp >> 49) - 6176; | |
2717 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
2718 | // set invalid flag | |
2719 | *pfpsf |= INVALID_EXCEPTION; | |
2720 | // return Integer Indefinite | |
2721 | res = 0x80000000; | |
2722 | BID_RETURN (res); | |
2723 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
2724 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
2725 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
2726 | // the cases that do not fit are identified here; the ones that fit | |
2727 | // fall through and will be handled with other cases further, | |
2728 | // under '1 <= q + exp <= 10' | |
2729 | if (x_sign) { // if n < 0 and q + exp = 10 | |
2730 | // if n <= -2^31 - 1 then n is too large | |
2731 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1 | |
2732 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34 | |
2733 | if (q <= 11) { | |
2734 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
2735 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
2736 | if (tmp64 >= 0x50000000aull) { | |
2737 | // set invalid flag | |
2738 | *pfpsf |= INVALID_EXCEPTION; | |
2739 | // return Integer Indefinite | |
2740 | res = 0x80000000; | |
2741 | BID_RETURN (res); | |
2742 | } | |
2743 | // else cases that can be rounded to a 32-bit int fall through | |
2744 | // to '1 <= q + exp <= 10' | |
2745 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2746 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a <=> | |
2747 | // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 23 | |
2748 | // (scale 2^31+1 up) | |
2749 | tmp64 = 0x50000000aull; | |
2750 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2751 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2752 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2753 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2754 | } | |
2755 | if (C1.w[1] > C.w[1] | |
2756 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
2757 | // set invalid flag | |
2758 | *pfpsf |= INVALID_EXCEPTION; | |
2759 | // return Integer Indefinite | |
2760 | res = 0x80000000; | |
2761 | BID_RETURN (res); | |
2762 | } | |
2763 | // else cases that can be rounded to a 32-bit int fall through | |
2764 | // to '1 <= q + exp <= 10' | |
2765 | } | |
2766 | } else { // if n > 0 and q + exp = 10 | |
2767 | // if n >= 2^31 then n is too large | |
2768 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31 | |
2769 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=34 | |
2770 | if (q <= 11) { | |
2771 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
2772 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
2773 | if (tmp64 >= 0x500000000ull) { | |
2774 | // set invalid flag | |
2775 | *pfpsf |= INVALID_EXCEPTION; | |
2776 | // return Integer Indefinite | |
2777 | res = 0x80000000; | |
2778 | BID_RETURN (res); | |
2779 | } | |
2780 | // else cases that can be rounded to a 32-bit int fall through | |
2781 | // to '1 <= q + exp <= 10' | |
2782 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
2783 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000 <=> | |
2784 | // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 23 | |
2785 | // (scale 2^31-1/2 up) | |
2786 | tmp64 = 0x500000000ull; | |
2787 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
2788 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
2789 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
2790 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
2791 | } | |
2792 | if (C1.w[1] > C.w[1] | |
2793 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
2794 | // set invalid flag | |
2795 | *pfpsf |= INVALID_EXCEPTION; | |
2796 | // return Integer Indefinite | |
2797 | res = 0x80000000; | |
2798 | BID_RETURN (res); | |
2799 | } | |
2800 | // else cases that can be rounded to a 32-bit int fall through | |
2801 | // to '1 <= q + exp <= 10' | |
200359e8 L |
2802 | } |
2803 | } | |
b2a00c89 L |
2804 | } |
2805 | // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31 | |
2806 | // Note: some of the cases tested for above fall through to this point | |
2807 | if ((q + exp) <= 0) { | |
2808 | // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1) | |
2809 | // set inexact flag | |
2810 | *pfpsf |= INEXACT_EXCEPTION; | |
2811 | // return 0 | |
2812 | res = 0x00000000; | |
2813 | BID_RETURN (res); | |
2814 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
2815 | // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded | |
2816 | // toward zero to a 32-bit signed integer | |
2817 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
2818 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
2819 | // chop off ind digits from the lower part of C1 | |
2820 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
2821 | tmp64 = C1.w[0]; | |
2822 | if (ind <= 19) { | |
2823 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
2824 | } else { | |
2825 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
2826 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
2827 | } | |
2828 | if (C1.w[0] < tmp64) | |
2829 | C1.w[1]++; | |
2830 | // calculate C* and f* | |
2831 | // C* is actually floor(C*) in this case | |
2832 | // C* and f* need shifting and masking, as shown by | |
2833 | // shiftright128[] and maskhigh128[] | |
2834 | // 1 <= x <= 33 | |
2835 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
2836 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
2837 | // the approximation of 10^(-x) was rounded up to 118 bits | |
2838 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
2839 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2840 | Cstar.w[1] = P256.w[3]; | |
2841 | Cstar.w[0] = P256.w[2]; | |
2842 | fstar.w[3] = 0; | |
2843 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
2844 | fstar.w[1] = P256.w[1]; | |
2845 | fstar.w[0] = P256.w[0]; | |
2846 | } else { // 22 <= ind - 1 <= 33 | |
2847 | Cstar.w[1] = 0; | |
2848 | Cstar.w[0] = P256.w[3]; | |
2849 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
2850 | fstar.w[2] = P256.w[2]; | |
2851 | fstar.w[1] = P256.w[1]; | |
2852 | fstar.w[0] = P256.w[0]; | |
2853 | } | |
2854 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
2855 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
2856 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
2857 | // if floor(C*) is even then C* = floor(C*) - logical right | |
2858 | // shift; C* has p decimal digits, correct by Prop. 1) | |
2859 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
2860 | // shift; C* has p decimal digits, correct by Pr. 1) | |
2861 | // else | |
2862 | // C* = floor(C*) (logical right shift; C has p decimal digits, | |
2863 | // correct by Property 1) | |
2864 | // n = C* * 10^(e+x) | |
2865 | ||
2866 | // shift right C* by Ex-128 = shiftright128[ind] | |
2867 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
2868 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
2869 | Cstar.w[0] = | |
2870 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
2871 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
2872 | } else { // 22 <= ind - 1 <= 33 | |
2873 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
2874 | } | |
2875 | // determine inexactness of the rounding of C* | |
2876 | // if (0 < f* - 1/2 < 10^(-x)) then | |
2877 | // the result is exact | |
2878 | // else // if (f* - 1/2 > T*) then | |
2879 | // the result is inexact | |
2880 | if (ind - 1 <= 2) { | |
2881 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
2882 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
2883 | if (tmp64 > ten2mk128trunc[ind - 1].w[1] | |
2884 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
2885 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2886 | // set the inexact flag |
2887 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2888 | } // else the result is exact |
2889 | } else { // the result is inexact; f2* <= 1/2 | |
2890 | // set the inexact flag | |
2891 | *pfpsf |= INEXACT_EXCEPTION; | |
2892 | is_inexact_gt_midpoint = 1; | |
2893 | } | |
2894 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
2895 | if (fstar.w[3] > 0x0 || | |
2896 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
2897 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
2898 | (fstar.w[1] || fstar.w[0]))) { | |
2899 | // f2* > 1/2 and the result may be exact | |
2900 | // Calculate f2* - 1/2 | |
2901 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
2902 | tmp64A = fstar.w[3]; | |
2903 | if (tmp64 > fstar.w[2]) | |
2904 | tmp64A--; | |
2905 | if (tmp64A || tmp64 | |
2906 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2907 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2908 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2909 | // set the inexact flag |
2910 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2911 | } // else the result is exact |
2912 | } else { // the result is inexact; f2* <= 1/2 | |
2913 | // set the inexact flag | |
2914 | *pfpsf |= INEXACT_EXCEPTION; | |
2915 | is_inexact_gt_midpoint = 1; | |
2916 | } | |
2917 | } else { // if 22 <= ind <= 33 | |
2918 | if (fstar.w[3] > onehalf128[ind - 1] || | |
2919 | (fstar.w[3] == onehalf128[ind - 1] && (fstar.w[2] || | |
2920 | fstar.w[1] | |
2921 | || fstar.w[0]))) { | |
2922 | // f2* > 1/2 and the result may be exact | |
2923 | // Calculate f2* - 1/2 | |
2924 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
2925 | if (tmp64 || fstar.w[2] | |
2926 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
2927 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2928 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
2929 | // set the inexact flag |
2930 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2931 | } // else the result is exact |
2932 | } else { // the result is inexact; f2* <= 1/2 | |
2933 | // set the inexact flag | |
2934 | *pfpsf |= INEXACT_EXCEPTION; | |
2935 | is_inexact_gt_midpoint = 1; | |
2936 | } | |
2937 | } | |
2938 | ||
2939 | // if the result was a midpoint it was rounded away from zero, so | |
2940 | // it will need a correction | |
2941 | // check for midpoints | |
2942 | if ((fstar.w[3] == 0) && (fstar.w[2] == 0) | |
2943 | && (fstar.w[1] || fstar.w[0]) | |
2944 | && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1] | |
2945 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
2946 | && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) { | |
2947 | // the result is a midpoint; round to nearest | |
2948 | if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD] | |
2949 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 | |
2950 | Cstar.w[0]--; // Cstar.w[0] is now even | |
2951 | is_inexact_gt_midpoint = 0; | |
2952 | } else { // else MP in [ODD, EVEN] | |
2953 | is_midpoint_lt_even = 1; | |
2954 | is_inexact_gt_midpoint = 0; | |
2955 | } | |
2956 | } | |
2957 | // general correction for RZ | |
2958 | if (is_midpoint_lt_even || is_inexact_gt_midpoint) { | |
2959 | Cstar.w[0] = Cstar.w[0] - 1; | |
2960 | } else { | |
2961 | ; // exact, the result is already correct | |
200359e8 | 2962 | } |
b2a00c89 L |
2963 | if (x_sign) |
2964 | res = -Cstar.w[0]; | |
2965 | else | |
2966 | res = Cstar.w[0]; | |
2967 | } else if (exp == 0) { | |
2968 | // 1 <= q <= 10 | |
2969 | // res = +/-C (exact) | |
2970 | if (x_sign) | |
2971 | res = -C1.w[0]; | |
2972 | else | |
2973 | res = C1.w[0]; | |
2974 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
2975 | // res = +/-C * 10^exp (exact) | |
2976 | if (x_sign) | |
2977 | res = -C1.w[0] * ten2k64[exp]; | |
2978 | else | |
2979 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
2980 | } |
2981 | } | |
b2a00c89 L |
2982 | } |
2983 | ||
2984 | BID_RETURN (res); | |
200359e8 L |
2985 | } |
2986 | ||
2987 | /***************************************************************************** | |
2988 | * BID128_to_int32_rninta | |
2989 | ****************************************************************************/ | |
2990 | ||
b2a00c89 L |
2991 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_rninta, |
2992 | x) | |
200359e8 | 2993 | |
b2a00c89 L |
2994 | int res; |
2995 | UINT64 x_sign; | |
2996 | UINT64 x_exp; | |
2997 | int exp; // unbiased exponent | |
200359e8 | 2998 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
2999 | UINT64 tmp64; |
3000 | BID_UI64DOUBLE tmp1; | |
3001 | unsigned int x_nr_bits; | |
3002 | int q, ind, shift; | |
3003 | UINT128 C1, C; | |
3004 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
3005 | UINT256 P256; | |
200359e8 L |
3006 | |
3007 | // unpack x | |
b2a00c89 L |
3008 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
3009 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
3010 | C1.w[1] = x.w[1] & MASK_COEFF; | |
3011 | C1.w[0] = x.w[0]; | |
200359e8 L |
3012 | |
3013 | // check for NaN or Infinity | |
b2a00c89 | 3014 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 3015 | // x is special |
b2a00c89 L |
3016 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
3017 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
3018 | // set invalid flag | |
3019 | *pfpsf |= INVALID_EXCEPTION; | |
3020 | // return Integer Indefinite | |
3021 | res = 0x80000000; | |
3022 | } else { // x is QNaN | |
3023 | // set invalid flag | |
3024 | *pfpsf |= INVALID_EXCEPTION; | |
3025 | // return Integer Indefinite | |
3026 | res = 0x80000000; | |
3027 | } | |
3028 | BID_RETURN (res); | |
3029 | } else { // x is not a NaN, so it must be infinity | |
3030 | if (!x_sign) { // x is +inf | |
3031 | // set invalid flag | |
3032 | *pfpsf |= INVALID_EXCEPTION; | |
3033 | // return Integer Indefinite | |
3034 | res = 0x80000000; | |
3035 | } else { // x is -inf | |
3036 | // set invalid flag | |
3037 | *pfpsf |= INVALID_EXCEPTION; | |
3038 | // return Integer Indefinite | |
3039 | res = 0x80000000; | |
200359e8 | 3040 | } |
b2a00c89 L |
3041 | BID_RETURN (res); |
3042 | } | |
3043 | } | |
200359e8 | 3044 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
3045 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
3046 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
3047 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
3048 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
3049 | res = 0x00000000; | |
3050 | BID_RETURN (res); | |
3051 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
3052 | // x is 0 | |
3053 | res = 0x00000000; | |
3054 | BID_RETURN (res); | |
3055 | } else { // x is not special and is not zero | |
3056 | ||
3057 | // q = nr. of decimal digits in x | |
3058 | // determine first the nr. of bits in x | |
3059 | if (C1.w[1] == 0) { | |
3060 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
3061 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
3062 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
3063 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
3064 | x_nr_bits = | |
3065 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3066 | } else { // x < 2^32 | |
3067 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
3068 | x_nr_bits = |
3069 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3070 | } | |
b2a00c89 L |
3071 | } else { // if x < 2^53 |
3072 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 3073 | x_nr_bits = |
b2a00c89 | 3074 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 3075 | } |
b2a00c89 L |
3076 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
3077 | tmp1.d = (double) C1.w[1]; // exact conversion | |
3078 | x_nr_bits = | |
3079 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3080 | } | |
3081 | q = nr_digits[x_nr_bits - 1].digits; | |
3082 | if (q == 0) { | |
3083 | q = nr_digits[x_nr_bits - 1].digits1; | |
3084 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
3085 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
3086 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
3087 | q++; | |
3088 | } | |
3089 | exp = (x_exp >> 49) - 6176; | |
3090 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
3091 | // set invalid flag | |
3092 | *pfpsf |= INVALID_EXCEPTION; | |
3093 | // return Integer Indefinite | |
3094 | res = 0x80000000; | |
3095 | BID_RETURN (res); | |
3096 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
3097 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
3098 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
3099 | // the cases that do not fit are identified here; the ones that fit | |
3100 | // fall through and will be handled with other cases further, | |
3101 | // under '1 <= q + exp <= 10' | |
3102 | if (x_sign) { // if n < 0 and q + exp = 10 | |
3103 | // if n <= -2^31 - 1/2 then n is too large | |
3104 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2 | |
3105 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=34 | |
3106 | if (q <= 11) { | |
3107 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
3108 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
3109 | if (tmp64 >= 0x500000005ull) { | |
3110 | // set invalid flag | |
3111 | *pfpsf |= INVALID_EXCEPTION; | |
3112 | // return Integer Indefinite | |
3113 | res = 0x80000000; | |
3114 | BID_RETURN (res); | |
3115 | } | |
3116 | // else cases that can be rounded to a 32-bit int fall through | |
3117 | // to '1 <= q + exp <= 10' | |
3118 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
3119 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005 <=> | |
3120 | // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23 | |
3121 | // (scale 2^31+1/2 up) | |
3122 | tmp64 = 0x500000005ull; | |
3123 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
3124 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
3125 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
3126 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
3127 | } | |
3128 | if (C1.w[1] > C.w[1] | |
3129 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
3130 | // set invalid flag | |
3131 | *pfpsf |= INVALID_EXCEPTION; | |
3132 | // return Integer Indefinite | |
3133 | res = 0x80000000; | |
3134 | BID_RETURN (res); | |
3135 | } | |
3136 | // else cases that can be rounded to a 32-bit int fall through | |
3137 | // to '1 <= q + exp <= 10' | |
3138 | } | |
3139 | } else { // if n > 0 and q + exp = 10 | |
3140 | // if n >= 2^31 - 1/2 then n is too large | |
3141 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2 | |
3142 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34 | |
3143 | if (q <= 11) { | |
3144 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
3145 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
3146 | if (tmp64 >= 0x4fffffffbull) { | |
3147 | // set invalid flag | |
3148 | *pfpsf |= INVALID_EXCEPTION; | |
3149 | // return Integer Indefinite | |
3150 | res = 0x80000000; | |
3151 | BID_RETURN (res); | |
3152 | } | |
3153 | // else cases that can be rounded to a 32-bit int fall through | |
3154 | // to '1 <= q + exp <= 10' | |
3155 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
3156 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=> | |
3157 | // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23 | |
3158 | // (scale 2^31-1/2 up) | |
3159 | tmp64 = 0x4fffffffbull; | |
3160 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
3161 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
3162 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
3163 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
3164 | } | |
3165 | if (C1.w[1] > C.w[1] | |
3166 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
3167 | // set invalid flag | |
3168 | *pfpsf |= INVALID_EXCEPTION; | |
3169 | // return Integer Indefinite | |
3170 | res = 0x80000000; | |
3171 | BID_RETURN (res); | |
3172 | } | |
3173 | // else cases that can be rounded to a 32-bit int fall through | |
3174 | // to '1 <= q + exp <= 10' | |
3175 | } | |
200359e8 | 3176 | } |
b2a00c89 L |
3177 | } |
3178 | // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2 | |
3179 | // Note: some of the cases tested for above fall through to this point | |
3180 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) | |
3181 | // return 0 | |
3182 | res = 0x00000000; | |
3183 | BID_RETURN (res); | |
3184 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) | |
3185 | // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1) | |
3186 | // res = 0 | |
3187 | // else | |
3188 | // res = +/-1 | |
3189 | ind = q - 1; | |
3190 | if (ind <= 18) { // 0 <= ind <= 18 | |
3191 | if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) { | |
3192 | res = 0x00000000; // return 0 | |
3193 | } else if (x_sign) { // n < 0 | |
3194 | res = 0xffffffff; // return -1 | |
3195 | } else { // n > 0 | |
3196 | res = 0x00000001; // return +1 | |
3197 | } | |
3198 | } else { // 19 <= ind <= 33 | |
3199 | if ((C1.w[1] < midpoint128[ind - 19].w[1]) | |
3200 | || ((C1.w[1] == midpoint128[ind - 19].w[1]) | |
3201 | && (C1.w[0] < midpoint128[ind - 19].w[0]))) { | |
3202 | res = 0x00000000; // return 0 | |
3203 | } else if (x_sign) { // n < 0 | |
3204 | res = 0xffffffff; // return -1 | |
3205 | } else { // n > 0 | |
3206 | res = 0x00000001; // return +1 | |
200359e8 L |
3207 | } |
3208 | } | |
b2a00c89 L |
3209 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) |
3210 | // -2^31-1/2 < x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded | |
3211 | // to nearest-away to a 32-bit signed integer | |
3212 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
3213 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
3214 | // chop off ind digits from the lower part of C1 | |
3215 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
3216 | tmp64 = C1.w[0]; | |
3217 | if (ind <= 19) { | |
3218 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
3219 | } else { | |
3220 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
3221 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
3222 | } | |
3223 | if (C1.w[0] < tmp64) | |
3224 | C1.w[1]++; | |
3225 | // calculate C* and f* | |
3226 | // C* is actually floor(C*) in this case | |
3227 | // C* and f* need shifting and masking, as shown by | |
3228 | // shiftright128[] and maskhigh128[] | |
3229 | // 1 <= x <= 33 | |
3230 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
3231 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
3232 | // the approximation of 10^(-x) was rounded up to 118 bits | |
3233 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
3234 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
3235 | Cstar.w[1] = P256.w[3]; | |
3236 | Cstar.w[0] = P256.w[2]; | |
3237 | } else { // 22 <= ind - 1 <= 33 | |
3238 | Cstar.w[1] = 0; | |
3239 | Cstar.w[0] = P256.w[3]; | |
3240 | } | |
3241 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
3242 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
3243 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
3244 | // if floor(C*) is even then C* = floor(C*) - logical right | |
3245 | // shift; C* has p decimal digits, correct by Prop. 1) | |
3246 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
3247 | // shift; C* has p decimal digits, correct by Pr. 1) | |
200359e8 | 3248 | // else |
b2a00c89 L |
3249 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
3250 | // correct by Property 1) | |
3251 | // n = C* * 10^(e+x) | |
3252 | ||
3253 | // shift right C* by Ex-128 = shiftright128[ind] | |
3254 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
3255 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
3256 | Cstar.w[0] = | |
3257 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
3258 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
3259 | } else { // 22 <= ind - 1 <= 33 | |
3260 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
200359e8 | 3261 | } |
b2a00c89 L |
3262 | // if the result was a midpoint, it was already rounded away from zero |
3263 | if (x_sign) | |
3264 | res = -Cstar.w[0]; | |
3265 | else | |
3266 | res = Cstar.w[0]; | |
3267 | // no need to check for midpoints - already rounded away from zero! | |
3268 | } else if (exp == 0) { | |
3269 | // 1 <= q <= 10 | |
3270 | // res = +/-C (exact) | |
3271 | if (x_sign) | |
3272 | res = -C1.w[0]; | |
3273 | else | |
3274 | res = C1.w[0]; | |
3275 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
3276 | // res = +/-C * 10^exp (exact) | |
3277 | if (x_sign) | |
3278 | res = -C1.w[0] * ten2k64[exp]; | |
3279 | else | |
3280 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
3281 | } |
3282 | } | |
b2a00c89 L |
3283 | } |
3284 | ||
3285 | BID_RETURN (res); | |
200359e8 L |
3286 | } |
3287 | ||
3288 | /***************************************************************************** | |
3289 | * BID128_to_int32_xrninta | |
3290 | ****************************************************************************/ | |
3291 | ||
b2a00c89 L |
3292 | BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (int, bid128_to_int32_xrninta, |
3293 | x) | |
200359e8 | 3294 | |
b2a00c89 L |
3295 | int res; |
3296 | UINT64 x_sign; | |
3297 | UINT64 x_exp; | |
3298 | int exp; // unbiased exponent | |
200359e8 | 3299 | // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64) |
b2a00c89 L |
3300 | UINT64 tmp64, tmp64A; |
3301 | BID_UI64DOUBLE tmp1; | |
3302 | unsigned int x_nr_bits; | |
3303 | int q, ind, shift; | |
3304 | UINT128 C1, C; | |
3305 | UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits | |
3306 | UINT256 fstar; | |
3307 | UINT256 P256; | |
200359e8 L |
3308 | |
3309 | // unpack x | |
b2a00c89 L |
3310 | x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
3311 | x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions | |
3312 | C1.w[1] = x.w[1] & MASK_COEFF; | |
3313 | C1.w[0] = x.w[0]; | |
200359e8 L |
3314 | |
3315 | // check for NaN or Infinity | |
b2a00c89 | 3316 | if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { |
200359e8 | 3317 | // x is special |
b2a00c89 L |
3318 | if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN |
3319 | if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN | |
3320 | // set invalid flag | |
3321 | *pfpsf |= INVALID_EXCEPTION; | |
3322 | // return Integer Indefinite | |
3323 | res = 0x80000000; | |
3324 | } else { // x is QNaN | |
3325 | // set invalid flag | |
3326 | *pfpsf |= INVALID_EXCEPTION; | |
3327 | // return Integer Indefinite | |
3328 | res = 0x80000000; | |
200359e8 | 3329 | } |
b2a00c89 L |
3330 | BID_RETURN (res); |
3331 | } else { // x is not a NaN, so it must be infinity | |
3332 | if (!x_sign) { // x is +inf | |
3333 | // set invalid flag | |
3334 | *pfpsf |= INVALID_EXCEPTION; | |
3335 | // return Integer Indefinite | |
3336 | res = 0x80000000; | |
3337 | } else { // x is -inf | |
3338 | // set invalid flag | |
3339 | *pfpsf |= INVALID_EXCEPTION; | |
3340 | // return Integer Indefinite | |
3341 | res = 0x80000000; | |
3342 | } | |
3343 | BID_RETURN (res); | |
3344 | } | |
3345 | } | |
200359e8 | 3346 | // check for non-canonical values (after the check for special values) |
b2a00c89 L |
3347 | if ((C1.w[1] > 0x0001ed09bead87c0ull) |
3348 | || (C1.w[1] == 0x0001ed09bead87c0ull | |
3349 | && (C1.w[0] > 0x378d8e63ffffffffull)) | |
3350 | || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) { | |
3351 | res = 0x00000000; | |
3352 | BID_RETURN (res); | |
3353 | } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { | |
3354 | // x is 0 | |
3355 | res = 0x00000000; | |
3356 | BID_RETURN (res); | |
3357 | } else { // x is not special and is not zero | |
3358 | ||
3359 | // q = nr. of decimal digits in x | |
3360 | // determine first the nr. of bits in x | |
3361 | if (C1.w[1] == 0) { | |
3362 | if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 | |
3363 | // split the 64-bit value in two 32-bit halves to avoid rounding errors | |
3364 | if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 | |
3365 | tmp1.d = (double) (C1.w[0] >> 32); // exact conversion | |
3366 | x_nr_bits = | |
3367 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3368 | } else { // x < 2^32 | |
3369 | tmp1.d = (double) (C1.w[0]); // exact conversion | |
200359e8 L |
3370 | x_nr_bits = |
3371 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3372 | } | |
b2a00c89 L |
3373 | } else { // if x < 2^53 |
3374 | tmp1.d = (double) C1.w[0]; // exact conversion | |
200359e8 | 3375 | x_nr_bits = |
b2a00c89 | 3376 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); |
200359e8 | 3377 | } |
b2a00c89 L |
3378 | } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) |
3379 | tmp1.d = (double) C1.w[1]; // exact conversion | |
3380 | x_nr_bits = | |
3381 | 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
3382 | } | |
3383 | q = nr_digits[x_nr_bits - 1].digits; | |
3384 | if (q == 0) { | |
3385 | q = nr_digits[x_nr_bits - 1].digits1; | |
3386 | if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi | |
3387 | || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi | |
3388 | && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) | |
3389 | q++; | |
3390 | } | |
3391 | exp = (x_exp >> 49) - 6176; | |
3392 | if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits) | |
3393 | // set invalid flag | |
3394 | *pfpsf |= INVALID_EXCEPTION; | |
3395 | // return Integer Indefinite | |
3396 | res = 0x80000000; | |
3397 | BID_RETURN (res); | |
3398 | } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1) | |
3399 | // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2... | |
3400 | // so x rounded to an integer may or may not fit in a signed 32-bit int | |
3401 | // the cases that do not fit are identified here; the ones that fit | |
3402 | // fall through and will be handled with other cases further, | |
3403 | // under '1 <= q + exp <= 10' | |
3404 | if (x_sign) { // if n < 0 and q + exp = 10 | |
3405 | // if n <= -2^31 - 1/2 then n is too large | |
3406 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2 | |
3407 | // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=34 | |
3408 | if (q <= 11) { | |
3409 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
3410 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
3411 | if (tmp64 >= 0x500000005ull) { | |
3412 | // set invalid flag | |
3413 | *pfpsf |= INVALID_EXCEPTION; | |
3414 | // return Integer Indefinite | |
3415 | res = 0x80000000; | |
3416 | BID_RETURN (res); | |
3417 | } | |
3418 | // else cases that can be rounded to a 32-bit int fall through | |
3419 | // to '1 <= q + exp <= 10' | |
3420 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
3421 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005 <=> | |
3422 | // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 23 | |
3423 | // (scale 2^31+1/2 up) | |
3424 | tmp64 = 0x500000005ull; | |
3425 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
3426 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
3427 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
3428 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
3429 | } | |
3430 | if (C1.w[1] > C.w[1] | |
3431 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
3432 | // set invalid flag | |
3433 | *pfpsf |= INVALID_EXCEPTION; | |
3434 | // return Integer Indefinite | |
3435 | res = 0x80000000; | |
3436 | BID_RETURN (res); | |
3437 | } | |
3438 | // else cases that can be rounded to a 32-bit int fall through | |
3439 | // to '1 <= q + exp <= 10' | |
3440 | } | |
3441 | } else { // if n > 0 and q + exp = 10 | |
3442 | // if n >= 2^31 - 1/2 then n is too large | |
3443 | // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2 | |
3444 | // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=34 | |
3445 | if (q <= 11) { | |
3446 | tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int | |
3447 | // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits) | |
3448 | if (tmp64 >= 0x4fffffffbull) { | |
3449 | // set invalid flag | |
3450 | *pfpsf |= INVALID_EXCEPTION; | |
3451 | // return Integer Indefinite | |
3452 | res = 0x80000000; | |
3453 | BID_RETURN (res); | |
3454 | } | |
3455 | // else cases that can be rounded to a 32-bit int fall through | |
3456 | // to '1 <= q + exp <= 10' | |
3457 | } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2 | |
3458 | // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb <=> | |
3459 | // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 23 | |
3460 | // (scale 2^31-1/2 up) | |
3461 | tmp64 = 0x4fffffffbull; | |
3462 | if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits | |
3463 | __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]); | |
3464 | } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits | |
3465 | __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]); | |
3466 | } | |
3467 | if (C1.w[1] > C.w[1] | |
3468 | || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) { | |
3469 | // set invalid flag | |
3470 | *pfpsf |= INVALID_EXCEPTION; | |
3471 | // return Integer Indefinite | |
3472 | res = 0x80000000; | |
3473 | BID_RETURN (res); | |
3474 | } | |
3475 | // else cases that can be rounded to a 32-bit int fall through | |
3476 | // to '1 <= q + exp <= 10' | |
3477 | } | |
200359e8 | 3478 | } |
b2a00c89 L |
3479 | } |
3480 | // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2 | |
3481 | // Note: some of the cases tested for above fall through to this point | |
3482 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) | |
3483 | // set inexact flag | |
3484 | *pfpsf |= INEXACT_EXCEPTION; | |
3485 | // return 0 | |
3486 | res = 0x00000000; | |
3487 | BID_RETURN (res); | |
3488 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) | |
3489 | // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1) | |
3490 | // res = 0 | |
3491 | // else | |
3492 | // res = +/-1 | |
3493 | ind = q - 1; | |
3494 | if (ind <= 18) { // 0 <= ind <= 18 | |
3495 | if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) { | |
3496 | res = 0x00000000; // return 0 | |
3497 | } else if (x_sign) { // n < 0 | |
3498 | res = 0xffffffff; // return -1 | |
3499 | } else { // n > 0 | |
3500 | res = 0x00000001; // return +1 | |
3501 | } | |
3502 | } else { // 19 <= ind <= 33 | |
3503 | if ((C1.w[1] < midpoint128[ind - 19].w[1]) | |
3504 | || ((C1.w[1] == midpoint128[ind - 19].w[1]) | |
3505 | && (C1.w[0] < midpoint128[ind - 19].w[0]))) { | |
3506 | res = 0x00000000; // return 0 | |
3507 | } else if (x_sign) { // n < 0 | |
3508 | res = 0xffffffff; // return -1 | |
3509 | } else { // n > 0 | |
3510 | res = 0x00000001; // return +1 | |
200359e8 L |
3511 | } |
3512 | } | |
b2a00c89 L |
3513 | // set inexact flag |
3514 | *pfpsf |= INEXACT_EXCEPTION; | |
3515 | } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9) | |
3516 | // -2^31-1/2 < x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded | |
3517 | // to nearest-away to a 32-bit signed integer | |
3518 | if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10 | |
3519 | ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x' | |
3520 | // chop off ind digits from the lower part of C1 | |
3521 | // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits | |
3522 | tmp64 = C1.w[0]; | |
3523 | if (ind <= 19) { | |
3524 | C1.w[0] = C1.w[0] + midpoint64[ind - 1]; | |
3525 | } else { | |
3526 | C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0]; | |
3527 | C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1]; | |
3528 | } | |
3529 | if (C1.w[0] < tmp64) | |
3530 | C1.w[1]++; | |
3531 | // calculate C* and f* | |
3532 | // C* is actually floor(C*) in this case | |
3533 | // C* and f* need shifting and masking, as shown by | |
3534 | // shiftright128[] and maskhigh128[] | |
3535 | // 1 <= x <= 33 | |
3536 | // kx = 10^(-x) = ten2mk128[ind - 1] | |
3537 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) | |
3538 | // the approximation of 10^(-x) was rounded up to 118 bits | |
3539 | __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]); | |
3540 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
3541 | Cstar.w[1] = P256.w[3]; | |
3542 | Cstar.w[0] = P256.w[2]; | |
3543 | fstar.w[3] = 0; | |
3544 | fstar.w[2] = P256.w[2] & maskhigh128[ind - 1]; | |
3545 | fstar.w[1] = P256.w[1]; | |
3546 | fstar.w[0] = P256.w[0]; | |
3547 | } else { // 22 <= ind - 1 <= 33 | |
3548 | Cstar.w[1] = 0; | |
3549 | Cstar.w[0] = P256.w[3]; | |
3550 | fstar.w[3] = P256.w[3] & maskhigh128[ind - 1]; | |
3551 | fstar.w[2] = P256.w[2]; | |
3552 | fstar.w[1] = P256.w[1]; | |
3553 | fstar.w[0] = P256.w[0]; | |
3554 | } | |
3555 | // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g. | |
3556 | // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999 | |
3557 | // if (0 < f* < 10^(-x)) then the result is a midpoint | |
3558 | // if floor(C*) is even then C* = floor(C*) - logical right | |
3559 | // shift; C* has p decimal digits, correct by Prop. 1) | |
3560 | // else if floor(C*) is odd C* = floor(C*)-1 (logical right | |
3561 | // shift; C* has p decimal digits, correct by Pr. 1) | |
200359e8 | 3562 | // else |
b2a00c89 L |
3563 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
3564 | // correct by Property 1) | |
3565 | // n = C* * 10^(e+x) | |
3566 | ||
3567 | // shift right C* by Ex-128 = shiftright128[ind] | |
3568 | shift = shiftright128[ind - 1]; // 0 <= shift <= 102 | |
3569 | if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21 | |
3570 | Cstar.w[0] = | |
3571 | (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift)); | |
3572 | // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift); | |
3573 | } else { // 22 <= ind - 1 <= 33 | |
3574 | Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38 | |
3575 | } | |
3576 | // if the result was a midpoint, it was already rounded away from zero | |
3577 | if (x_sign) | |
3578 | res = -Cstar.w[0]; | |
3579 | else | |
3580 | res = Cstar.w[0]; | |
3581 | // determine inexactness of the rounding of C* | |
3582 | // if (0 < f* - 1/2 < 10^(-x)) then | |
3583 | // the result is exact | |
3584 | // else // if (f* - 1/2 > T*) then | |
3585 | // the result is inexact | |
3586 | if (ind - 1 <= 2) { | |
3587 | if (fstar.w[1] > 0x8000000000000000ull || (fstar.w[1] == 0x8000000000000000ull && fstar.w[0] > 0x0ull)) { // f* > 1/2 and the result may be exact | |
3588 | tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2 | |
3589 | if ((tmp64 > ten2mk128trunc[ind - 1].w[1] | |
3590 | || (tmp64 == ten2mk128trunc[ind - 1].w[1] | |
3591 | && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0]))) { | |
200359e8 L |
3592 | // set the inexact flag |
3593 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
3594 | } // else the result is exact |
3595 | } else { // the result is inexact; f2* <= 1/2 | |
3596 | // set the inexact flag | |
3597 | *pfpsf |= INEXACT_EXCEPTION; | |
3598 | } | |
3599 | } else if (ind - 1 <= 21) { // if 3 <= ind <= 21 | |
3600 | if (fstar.w[3] > 0x0 || | |
3601 | (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) || | |
3602 | (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] && | |
3603 | (fstar.w[1] || fstar.w[0]))) { | |
3604 | // f2* > 1/2 and the result may be exact | |
3605 | // Calculate f2* - 1/2 | |
3606 | tmp64 = fstar.w[2] - onehalf128[ind - 1]; | |
3607 | tmp64A = fstar.w[3]; | |
3608 | if (tmp64 > fstar.w[2]) | |
3609 | tmp64A--; | |
3610 | if (tmp64A || tmp64 | |
3611 | || fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
3612 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
3613 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
3614 | // set the inexact flag |
3615 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
3616 | } // else the result is exact |
3617 | } else { // the result is inexact; f2* <= 1/2 | |
3618 | // set the inexact flag | |
3619 | *pfpsf |= INEXACT_EXCEPTION; | |
3620 | } | |
3621 | } else { // if 22 <= ind <= 33 | |
3622 | if (fstar.w[3] > onehalf128[ind - 1] || | |
3623 | (fstar.w[3] == onehalf128[ind - 1] && | |
3624 | (fstar.w[2] || fstar.w[1] || fstar.w[0]))) { | |
3625 | // f2* > 1/2 and the result may be exact | |
3626 | // Calculate f2* - 1/2 | |
3627 | tmp64 = fstar.w[3] - onehalf128[ind - 1]; | |
3628 | if (tmp64 || fstar.w[2] || | |
3629 | fstar.w[1] > ten2mk128trunc[ind - 1].w[1] | |
3630 | || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1] | |
3631 | && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) { | |
200359e8 L |
3632 | // set the inexact flag |
3633 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
3634 | } // else the result is exact |
3635 | } else { // the result is inexact; f2* <= 1/2 | |
3636 | // set the inexact flag | |
3637 | *pfpsf |= INEXACT_EXCEPTION; | |
3638 | } | |
200359e8 | 3639 | } |
b2a00c89 L |
3640 | // no need to check for midpoints - already rounded away from zero! |
3641 | } else if (exp == 0) { | |
3642 | // 1 <= q <= 10 | |
3643 | // res = +/-C (exact) | |
3644 | if (x_sign) | |
3645 | res = -C1.w[0]; | |
3646 | else | |
3647 | res = C1.w[0]; | |
3648 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
3649 | // res = +/-C * 10^exp (exact) | |
3650 | if (x_sign) | |
3651 | res = -C1.w[0] * ten2k64[exp]; | |
3652 | else | |
3653 | res = C1.w[0] * ten2k64[exp]; | |
200359e8 L |
3654 | } |
3655 | } | |
b2a00c89 L |
3656 | } |
3657 | ||
3658 | BID_RETURN (res); | |
200359e8 | 3659 | } |