]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgcc/config/libbid/bid128_to_int32.c
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid128_to_int32.c
CommitLineData
8d9254fc 1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
200359e8
L
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
748086b7 7Software Foundation; either version 3, or (at your option) any later
200359e8
L
8version.
9
200359e8
L
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13for more details.
14
748086b7
JJ
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see 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 30BID128_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
47x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
48x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
49C1.w[1] = x.w[1] & MASK_COEFF;
50C1.w[0] = x.w[0];
200359e8
L
51
52 // check for NaN or Infinity
b2a00c89 53if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 54 // x is special
b2a00c89
L
55if ((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
84if ((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
343BID_RETURN (res);
200359e8
L
344}
345
346/*****************************************************************************
347 * BID128_to_int32_xrnint
348 ****************************************************************************/
349
b2a00c89
L
350BID128_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
368x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
369x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
370C1.w[1] = x.w[1] & MASK_COEFF;
371C1.w[0] = x.w[0];
200359e8
L
372
373 // check for NaN or Infinity
b2a00c89 374if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 375 // x is special
b2a00c89
L
376if ((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
405if ((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
727BID_RETURN (res);
200359e8
L
728}
729
730/*****************************************************************************
731 * BID128_to_int32_floor
732 ****************************************************************************/
733
b2a00c89 734BID128_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
755x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
756x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
757C1.w[1] = x.w[1] & MASK_COEFF;
758C1.w[0] = x.w[0];
200359e8
L
759
760 // check for NaN or Infinity
b2a00c89 761if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 762 // x is special
b2a00c89
L
763if ((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
792if ((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
1100BID_RETURN (res);
200359e8
L
1101}
1102
1103
1104/*****************************************************************************
1105 * BID128_to_int32_xfloor
1106 ****************************************************************************/
1107
b2a00c89
L
1108BID128_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
1130x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1131x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1132C1.w[1] = x.w[1] & MASK_COEFF;
1133C1.w[0] = x.w[0];
200359e8
L
1134
1135 // check for NaN or Infinity
b2a00c89 1136if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1137 // x is special
b2a00c89
L
1138if ((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
1167if ((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
1489BID_RETURN (res);
200359e8
L
1490}
1491
1492/*****************************************************************************
1493 * BID128_to_int32_ceil
1494 ****************************************************************************/
1495
b2a00c89 1496BID128_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
1517x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1518x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1519C1.w[1] = x.w[1] & MASK_COEFF;
1520C1.w[0] = x.w[0];
200359e8
L
1521
1522 // check for NaN or Infinity
b2a00c89 1523if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1524 // x is special
b2a00c89
L
1525if ((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
1554if ((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
1862BID_RETURN (res);
200359e8
L
1863}
1864
1865/*****************************************************************************
1866 * BID128_to_int32_xceil
1867 ****************************************************************************/
1868
b2a00c89 1869BID128_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
1890x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1891x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
1892C1.w[1] = x.w[1] & MASK_COEFF;
1893C1.w[0] = x.w[0];
200359e8
L
1894
1895 // check for NaN or Infinity
b2a00c89 1896if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 1897 // x is special
b2a00c89
L
1898if ((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
1927if ((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
2249BID_RETURN (res);
200359e8
L
2250}
2251
2252/*****************************************************************************
2253 * BID128_to_int32_int
2254 ****************************************************************************/
2255
b2a00c89 2256BID128_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
2275x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2276x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2277C1.w[1] = x.w[1] & MASK_COEFF;
2278C1.w[0] = x.w[0];
200359e8
L
2279
2280 // check for NaN or Infinity
b2a00c89 2281if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2282 // x is special
b2a00c89
L
2283if ((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
2312if ((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
2609BID_RETURN (res);
200359e8
L
2610}
2611
2612/*****************************************************************************
2613 * BID128_to_int32_xint
2614 ****************************************************************************/
2615
b2a00c89 2616BID128_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
2635x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2636x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
2637C1.w[1] = x.w[1] & MASK_COEFF;
2638C1.w[0] = x.w[0];
200359e8
L
2639
2640 // check for NaN or Infinity
b2a00c89 2641if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 2642 // x is special
b2a00c89
L
2643if ((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
2672if ((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
2984BID_RETURN (res);
200359e8
L
2985}
2986
2987/*****************************************************************************
2988 * BID128_to_int32_rninta
2989 ****************************************************************************/
2990
b2a00c89
L
2991BID128_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
3008x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
3009x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
3010C1.w[1] = x.w[1] & MASK_COEFF;
3011C1.w[0] = x.w[0];
200359e8
L
3012
3013 // check for NaN or Infinity
b2a00c89 3014if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 3015 // x is special
b2a00c89
L
3016if ((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
3045if ((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
3285BID_RETURN (res);
200359e8
L
3286}
3287
3288/*****************************************************************************
3289 * BID128_to_int32_xrninta
3290 ****************************************************************************/
3291
b2a00c89
L
3292BID128_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
3310x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
3311x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
3312C1.w[1] = x.w[1] & MASK_COEFF;
3313C1.w[0] = x.w[0];
200359e8
L
3314
3315 // check for NaN or Infinity
b2a00c89 3316if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
200359e8 3317 // x is special
b2a00c89
L
3318if ((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
3347if ((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
3658BID_RETURN (res);
200359e8 3659}