1 /* Copyright (C) 2007-2023 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
18
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
23
24 #include "bid_internal.h"
25
26 /*****************************************************************************
27 * BID64_to_int32_rnint
28 ****************************************************************************/
29
30 #if DECIMAL_CALL_BY_REFERENCE
31 void
32 bid64_to_int32_rnint (int *pres,
33 UINT64 *
34 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
35 _EXC_INFO_PARAM) {
36 UINT64 x = *px;
37 #else
38 int
39 bid64_to_int32_rnint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
40 _EXC_INFO_PARAM) {
41 #endif
42 int res;
43 UINT64 x_sign;
44 UINT64 x_exp;
45 int exp; // unbiased exponent
46 // Note: C1 represents x_significand (UINT64)
47 UINT64 tmp64;
48 BID_UI64DOUBLE tmp1;
49 unsigned int x_nr_bits;
50 int q, ind, shift;
51 UINT64 C1;
52 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
53 UINT128 fstar;
54 UINT128 P128;
55
56 // check for NaN or Infinity
57 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
58 // set invalid flag
59 *pfpsf |= INVALID_EXCEPTION;
60 // return Integer Indefinite
61 res = 0x80000000;
62 BID_RETURN (res);
63 }
64 // unpack x
65 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
66 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
67 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
68 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
69 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
70 if (C1 > 9999999999999999ull) { // non-canonical
71 x_exp = 0;
72 C1 = 0;
73 }
74 } else {
75 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
76 C1 = x & MASK_BINARY_SIG1;
77 }
78
79 // check for zeros (possibly from non-canonical values)
80 if (C1 == 0x0ull) {
81 // x is 0
82 res = 0x00000000;
83 BID_RETURN (res);
84 }
85 // x is not special and is not zero
86
87 // q = nr. of decimal digits in x (1 <= q <= 54)
88 // determine first the nr. of bits in x
89 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
90 // split the 64-bit value in two 32-bit halves to avoid rounding errors
91 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
92 tmp1.d = (double) (C1 >> 32); // exact conversion
93 x_nr_bits =
94 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
95 } else { // x < 2^32
96 tmp1.d = (double) C1; // exact conversion
97 x_nr_bits =
98 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
99 }
100 } else { // if x < 2^53
101 tmp1.d = (double) C1; // exact conversion
102 x_nr_bits =
103 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
104 }
105 q = nr_digits[x_nr_bits - 1].digits;
106 if (q == 0) {
107 q = nr_digits[x_nr_bits - 1].digits1;
108 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
109 q++;
110 }
111 exp = x_exp - 398; // unbiased exponent
112
113 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
114 // set invalid flag
115 *pfpsf |= INVALID_EXCEPTION;
116 // return Integer Indefinite
117 res = 0x80000000;
118 BID_RETURN (res);
119 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
120 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
121 // so x rounded to an integer may or may not fit in a signed 32-bit int
122 // the cases that do not fit are identified here; the ones that fit
123 // fall through and will be handled with other cases further,
124 // under '1 <= q + exp <= 10'
125 if (x_sign) { // if n < 0 and q + exp = 10
126 // if n < -2^31 - 1/2 then n is too large
127 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
128 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=16
129 // <=> C * 10^(11-q) > 0x500000005, 1<=q<=16
130 if (q <= 11) {
131 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
132 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
133 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
134 if (tmp64 > 0x500000005ull) {
135 // set invalid flag
136 *pfpsf |= INVALID_EXCEPTION;
137 // return Integer Indefinite
138 res = 0x80000000;
139 BID_RETURN (res);
140 }
141 // else cases that can be rounded to a 32-bit int fall through
142 // to '1 <= q + exp <= 10'
143 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
144 // C * 10^(11-q) > 0x500000005 <=>
145 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
146 // (scale 2^31+1/2 up)
147 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
148 tmp64 = 0x500000005ull * ten2k64[q - 11];
149 if (C1 > tmp64) {
150 // set invalid flag
151 *pfpsf |= INVALID_EXCEPTION;
152 // return Integer Indefinite
153 res = 0x80000000;
154 BID_RETURN (res);
155 }
156 // else cases that can be rounded to a 32-bit int fall through
157 // to '1 <= q + exp <= 10'
158 }
159 } else { // if n > 0 and q + exp = 10
160 // if n >= 2^31 - 1/2 then n is too large
161 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
162 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
163 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
164 if (q <= 11) {
165 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
166 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
167 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
168 if (tmp64 >= 0x4fffffffbull) {
169 // set invalid flag
170 *pfpsf |= INVALID_EXCEPTION;
171 // return Integer Indefinite
172 res = 0x80000000;
173 BID_RETURN (res);
174 }
175 // else cases that can be rounded to a 32-bit int fall through
176 // to '1 <= q + exp <= 10'
177 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
178 // C * 10^(11-q) >= 0x4fffffffb <=>
179 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
180 // (scale 2^31-1/2 up)
181 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
182 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
183 if (C1 >= tmp64) {
184 // set invalid flag
185 *pfpsf |= INVALID_EXCEPTION;
186 // return Integer Indefinite
187 res = 0x80000000;
188 BID_RETURN (res);
189 }
190 // else cases that can be rounded to a 32-bit int fall through
191 // to '1 <= q + exp <= 10'
192 }
193 }
194 }
195 // n is not too large to be converted to int32: -2^31 - 1/2 <= n < 2^31 - 1/2
196 // Note: some of the cases tested for above fall through to this point
197 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
198 // return 0
199 res = 0x00000000;
200 BID_RETURN (res);
201 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
202 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
203 // res = 0
204 // else
205 // res = +/-1
206 ind = q - 1;
207 if (C1 <= midpoint64[ind]) {
208 res = 0x00000000; // return 0
209 } else if (x_sign) { // n < 0
210 res = 0xffffffff; // return -1
211 } else { // n > 0
212 res = 0x00000001; // return +1
213 }
214 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
215 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
216 // to nearest to a 32-bit signed integer
217 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
218 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
219 // chop off ind digits from the lower part of C1
220 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
221 C1 = C1 + midpoint64[ind - 1];
222 // calculate C* and f*
223 // C* is actually floor(C*) in this case
224 // C* and f* need shifting and masking, as shown by
225 // shiftright128[] and maskhigh128[]
226 // 1 <= x <= 15
227 // kx = 10^(-x) = ten2mk64[ind - 1]
228 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
229 // the approximation of 10^(-x) was rounded up to 54 bits
230 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
231 Cstar = P128.w[1];
232 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
233 fstar.w[0] = P128.w[0];
234 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
235 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
236 // if (0 < f* < 10^(-x)) then the result is a midpoint
237 // if floor(C*) is even then C* = floor(C*) - logical right
238 // shift; C* has p decimal digits, correct by Prop. 1)
239 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
240 // shift; C* has p decimal digits, correct by Pr. 1)
241 // else
242 // C* = floor(C*) (logical right shift; C has p decimal digits,
243 // correct by Property 1)
244 // n = C* * 10^(e+x)
245
246 // shift right C* by Ex-64 = shiftright128[ind]
247 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
248 Cstar = Cstar >> shift;
249
250 // if the result was a midpoint it was rounded away from zero, so
251 // it will need a correction
252 // check for midpoints
253 if ((fstar.w[1] == 0) && fstar.w[0]
254 && (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
255 // ten2mk128trunc[ind -1].w[1] is identical to
256 // ten2mk128[ind -1].w[1]
257 // the result is a midpoint; round to nearest
258 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
259 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
260 Cstar--; // Cstar is now even
261 } // else MP in [ODD, EVEN]
262 }
263 if (x_sign)
264 res = -Cstar;
265 else
266 res = Cstar;
267 } else if (exp == 0) {
268 // 1 <= q <= 10
269 // res = +/-C (exact)
270 if (x_sign)
271 res = -C1;
272 else
273 res = C1;
274 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
275 // res = +/-C * 10^exp (exact)
276 if (x_sign)
277 res = -C1 * ten2k64[exp];
278 else
279 res = C1 * ten2k64[exp];
280 }
281 }
282 BID_RETURN (res);
283 }
284
285 /*****************************************************************************
286 * BID64_to_int32_xrnint
287 ****************************************************************************/
288
289 #if DECIMAL_CALL_BY_REFERENCE
290 void
291 bid64_to_int32_xrnint (int *pres,
292 UINT64 *
293 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
294 _EXC_INFO_PARAM) {
295 UINT64 x = *px;
296 #else
297 int
298 bid64_to_int32_xrnint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
299 _EXC_INFO_PARAM) {
300 #endif
301 int res;
302 UINT64 x_sign;
303 UINT64 x_exp;
304 int exp; // unbiased exponent
305 // Note: C1 represents x_significand (UINT64)
306 UINT64 tmp64;
307 BID_UI64DOUBLE tmp1;
308 unsigned int x_nr_bits;
309 int q, ind, shift;
310 UINT64 C1;
311 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
312 UINT128 fstar;
313 UINT128 P128;
314
315 // check for NaN or Infinity
316 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
317 // set invalid flag
318 *pfpsf |= INVALID_EXCEPTION;
319 // return Integer Indefinite
320 res = 0x80000000;
321 BID_RETURN (res);
322 }
323 // unpack x
324 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
325 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
326 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
327 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
328 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
329 if (C1 > 9999999999999999ull) { // non-canonical
330 x_exp = 0;
331 C1 = 0;
332 }
333 } else {
334 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
335 C1 = x & MASK_BINARY_SIG1;
336 }
337
338 // check for zeros (possibly from non-canonical values)
339 if (C1 == 0x0ull) {
340 // x is 0
341 res = 0x00000000;
342 BID_RETURN (res);
343 }
344 // x is not special and is not zero
345
346 // q = nr. of decimal digits in x (1 <= q <= 54)
347 // determine first the nr. of bits in x
348 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
349 // split the 64-bit value in two 32-bit halves to avoid rounding errors
350 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
351 tmp1.d = (double) (C1 >> 32); // exact conversion
352 x_nr_bits =
353 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
354 } else { // x < 2^32
355 tmp1.d = (double) C1; // exact conversion
356 x_nr_bits =
357 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
358 }
359 } else { // if x < 2^53
360 tmp1.d = (double) C1; // exact conversion
361 x_nr_bits =
362 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
363 }
364 q = nr_digits[x_nr_bits - 1].digits;
365 if (q == 0) {
366 q = nr_digits[x_nr_bits - 1].digits1;
367 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
368 q++;
369 }
370 exp = x_exp - 398; // unbiased exponent
371
372 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
373 // set invalid flag
374 *pfpsf |= INVALID_EXCEPTION;
375 // return Integer Indefinite
376 res = 0x80000000;
377 BID_RETURN (res);
378 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
379 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
380 // so x rounded to an integer may or may not fit in a signed 32-bit int
381 // the cases that do not fit are identified here; the ones that fit
382 // fall through and will be handled with other cases further,
383 // under '1 <= q + exp <= 10'
384 if (x_sign) { // if n < 0 and q + exp = 10
385 // if n < -2^31 - 1/2 then n is too large
386 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
387 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=16
388 // <=> C * 10^(11-q) > 0x500000005, 1<=q<=16
389 if (q <= 11) {
390 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
391 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
392 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
393 if (tmp64 > 0x500000005ull) {
394 // set invalid flag
395 *pfpsf |= INVALID_EXCEPTION;
396 // return Integer Indefinite
397 res = 0x80000000;
398 BID_RETURN (res);
399 }
400 // else cases that can be rounded to a 32-bit int fall through
401 // to '1 <= q + exp <= 10'
402 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
403 // C * 10^(11-q) > 0x500000005 <=>
404 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
405 // (scale 2^31+1/2 up)
406 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
407 tmp64 = 0x500000005ull * ten2k64[q - 11];
408 if (C1 > tmp64) {
409 // set invalid flag
410 *pfpsf |= INVALID_EXCEPTION;
411 // return Integer Indefinite
412 res = 0x80000000;
413 BID_RETURN (res);
414 }
415 // else cases that can be rounded to a 32-bit int fall through
416 // to '1 <= q + exp <= 10'
417 }
418 } else { // if n > 0 and q + exp = 10
419 // if n >= 2^31 - 1/2 then n is too large
420 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
421 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
422 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
423 if (q <= 11) {
424 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
425 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
426 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
427 if (tmp64 >= 0x4fffffffbull) {
428 // set invalid flag
429 *pfpsf |= INVALID_EXCEPTION;
430 // return Integer Indefinite
431 res = 0x80000000;
432 BID_RETURN (res);
433 }
434 // else cases that can be rounded to a 32-bit int fall through
435 // to '1 <= q + exp <= 10'
436 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
437 // C * 10^(11-q) >= 0x4fffffffb <=>
438 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
439 // (scale 2^31-1/2 up)
440 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
441 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
442 if (C1 >= tmp64) {
443 // set invalid flag
444 *pfpsf |= INVALID_EXCEPTION;
445 // return Integer Indefinite
446 res = 0x80000000;
447 BID_RETURN (res);
448 }
449 // else cases that can be rounded to a 32-bit int fall through
450 // to '1 <= q + exp <= 10'
451 }
452 }
453 }
454 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
455 // Note: some of the cases tested for above fall through to this point
456 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
457 // set inexact flag
458 *pfpsf |= INEXACT_EXCEPTION;
459 // return 0
460 res = 0x00000000;
461 BID_RETURN (res);
462 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
463 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
464 // res = 0
465 // else
466 // res = +/-1
467 ind = q - 1;
468 if (C1 <= midpoint64[ind]) {
469 res = 0x00000000; // return 0
470 } else if (x_sign) { // n < 0
471 res = 0xffffffff; // return -1
472 } else { // n > 0
473 res = 0x00000001; // return +1
474 }
475 // set inexact flag
476 *pfpsf |= INEXACT_EXCEPTION;
477 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
478 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
479 // to nearest to a 32-bit signed integer
480 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
481 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
482 // chop off ind digits from the lower part of C1
483 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
484 C1 = C1 + midpoint64[ind - 1];
485 // calculate C* and f*
486 // C* is actually floor(C*) in this case
487 // C* and f* need shifting and masking, as shown by
488 // shiftright128[] and maskhigh128[]
489 // 1 <= x <= 15
490 // kx = 10^(-x) = ten2mk64[ind - 1]
491 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
492 // the approximation of 10^(-x) was rounded up to 54 bits
493 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
494 Cstar = P128.w[1];
495 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
496 fstar.w[0] = P128.w[0];
497 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
498 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
499 // if (0 < f* < 10^(-x)) then the result is a midpoint
500 // if floor(C*) is even then C* = floor(C*) - logical right
501 // shift; C* has p decimal digits, correct by Prop. 1)
502 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
503 // shift; C* has p decimal digits, correct by Pr. 1)
504 // else
505 // C* = floor(C*) (logical right shift; C has p decimal digits,
506 // correct by Property 1)
507 // n = C* * 10^(e+x)
508
509 // shift right C* by Ex-64 = shiftright128[ind]
510 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
511 Cstar = Cstar >> shift;
512 // determine inexactness of the rounding of C*
513 // if (0 < f* - 1/2 < 10^(-x)) then
514 // the result is exact
515 // else // if (f* - 1/2 > T*) then
516 // the result is inexact
517 if (ind - 1 <= 2) {
518 if (fstar.w[0] > 0x8000000000000000ull) {
519 // f* > 1/2 and the result may be exact
520 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
521 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
522 // ten2mk128trunc[ind -1].w[1] is identical to
523 // ten2mk128[ind -1].w[1]
524 // set the inexact flag
525 *pfpsf |= INEXACT_EXCEPTION;
526 } // else the result is exact
527 } else { // the result is inexact; f2* <= 1/2
528 // set the inexact flag
529 *pfpsf |= INEXACT_EXCEPTION;
530 }
531 } else { // if 3 <= ind - 1 <= 14
532 if (fstar.w[1] > onehalf128[ind - 1] ||
533 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
534 // f2* > 1/2 and the result may be exact
535 // Calculate f2* - 1/2
536 tmp64 = fstar.w[1] - onehalf128[ind - 1];
537 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
538 // ten2mk128trunc[ind -1].w[1] is identical to
539 // ten2mk128[ind -1].w[1]
540 // set the inexact flag
541 *pfpsf |= INEXACT_EXCEPTION;
542 } // else the result is exact
543 } else { // the result is inexact; f2* <= 1/2
544 // set the inexact flag
545 *pfpsf |= INEXACT_EXCEPTION;
546 }
547 }
548
549 // if the result was a midpoint it was rounded away from zero, so
550 // it will need a correction
551 // check for midpoints
552 if ((fstar.w[1] == 0) && fstar.w[0]
553 && (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
554 // ten2mk128trunc[ind -1].w[1] is identical to
555 // ten2mk128[ind -1].w[1]
556 // the result is a midpoint; round to nearest
557 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
558 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
559 Cstar--; // Cstar is now even
560 } // else MP in [ODD, EVEN]
561 }
562 if (x_sign)
563 res = -Cstar;
564 else
565 res = Cstar;
566 } else if (exp == 0) {
567 // 1 <= q <= 10
568 // res = +/-C (exact)
569 if (x_sign)
570 res = -C1;
571 else
572 res = C1;
573 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
574 // res = +/-C * 10^exp (exact)
575 if (x_sign)
576 res = -C1 * ten2k64[exp];
577 else
578 res = C1 * ten2k64[exp];
579 }
580 }
581 BID_RETURN (res);
582 }
583
584 /*****************************************************************************
585 * BID64_to_int32_floor
586 ****************************************************************************/
587
588 #if DECIMAL_CALL_BY_REFERENCE
589 void
590 bid64_to_int32_floor (int *pres,
591 UINT64 *
592 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
593 _EXC_INFO_PARAM) {
594 UINT64 x = *px;
595 #else
596 int
597 bid64_to_int32_floor (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
598 _EXC_INFO_PARAM) {
599 #endif
600 int res;
601 UINT64 x_sign;
602 UINT64 x_exp;
603 int exp; // unbiased exponent
604 // Note: C1 represents x_significand (UINT64)
605 UINT64 tmp64;
606 BID_UI64DOUBLE tmp1;
607 unsigned int x_nr_bits;
608 int q, ind, shift;
609 UINT64 C1;
610 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
611 UINT128 fstar;
612 UINT128 P128;
613
614 // check for NaN or Infinity
615 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
616 // set invalid flag
617 *pfpsf |= INVALID_EXCEPTION;
618 // return Integer Indefinite
619 res = 0x80000000;
620 BID_RETURN (res);
621 }
622 // unpack x
623 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
624 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
625 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
626 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
627 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
628 if (C1 > 9999999999999999ull) { // non-canonical
629 x_exp = 0;
630 C1 = 0;
631 }
632 } else {
633 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
634 C1 = x & MASK_BINARY_SIG1;
635 }
636
637 // check for zeros (possibly from non-canonical values)
638 if (C1 == 0x0ull) {
639 // x is 0
640 res = 0x00000000;
641 BID_RETURN (res);
642 }
643 // x is not special and is not zero
644
645 // q = nr. of decimal digits in x (1 <= q <= 54)
646 // determine first the nr. of bits in x
647 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
648 // split the 64-bit value in two 32-bit halves to avoid rounding errors
649 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
650 tmp1.d = (double) (C1 >> 32); // exact conversion
651 x_nr_bits =
652 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
653 } else { // x < 2^32
654 tmp1.d = (double) C1; // exact conversion
655 x_nr_bits =
656 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
657 }
658 } else { // if x < 2^53
659 tmp1.d = (double) C1; // exact conversion
660 x_nr_bits =
661 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
662 }
663 q = nr_digits[x_nr_bits - 1].digits;
664 if (q == 0) {
665 q = nr_digits[x_nr_bits - 1].digits1;
666 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
667 q++;
668 }
669 exp = x_exp - 398; // unbiased exponent
670
671 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
672 // set invalid flag
673 *pfpsf |= INVALID_EXCEPTION;
674 // return Integer Indefinite
675 res = 0x80000000;
676 BID_RETURN (res);
677 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
678 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
679 // so x rounded to an integer may or may not fit in a signed 32-bit int
680 // the cases that do not fit are identified here; the ones that fit
681 // fall through and will be handled with other cases further,
682 // under '1 <= q + exp <= 10'
683 if (x_sign) { // if n < 0 and q + exp = 10
684 // if n < -2^31 then n is too large
685 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
686 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=16
687 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
688 if (q <= 11) {
689 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
690 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
691 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
692 if (tmp64 > 0x500000000ull) {
693 // set invalid flag
694 *pfpsf |= INVALID_EXCEPTION;
695 // return Integer Indefinite
696 res = 0x80000000;
697 BID_RETURN (res);
698 }
699 // else cases that can be rounded to a 32-bit int fall through
700 // to '1 <= q + exp <= 10'
701 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
702 // C * 10^(11-q) > 0x500000000 <=>
703 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
704 // (scale 2^31+1 up)
705 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
706 tmp64 = 0x500000000ull * ten2k64[q - 11];
707 if (C1 > tmp64) {
708 // set invalid flag
709 *pfpsf |= INVALID_EXCEPTION;
710 // return Integer Indefinite
711 res = 0x80000000;
712 BID_RETURN (res);
713 }
714 // else cases that can be rounded to a 32-bit int fall through
715 // to '1 <= q + exp <= 10'
716 }
717 } else { // if n > 0 and q + exp = 10
718 // if n >= 2^31 then n is too large
719 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
720 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
721 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
722 if (q <= 11) {
723 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
724 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
725 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
726 if (tmp64 >= 0x500000000ull) {
727 // set invalid flag
728 *pfpsf |= INVALID_EXCEPTION;
729 // return Integer Indefinite
730 res = 0x80000000;
731 BID_RETURN (res);
732 }
733 // else cases that can be rounded to a 32-bit int fall through
734 // to '1 <= q + exp <= 10'
735 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
736 // C * 10^(11-q) >= 0x500000000 <=>
737 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
738 // (scale 2^31-1 up)
739 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
740 tmp64 = 0x500000000ull * ten2k64[q - 11];
741 if (C1 >= tmp64) {
742 // set invalid flag
743 *pfpsf |= INVALID_EXCEPTION;
744 // return Integer Indefinite
745 res = 0x80000000;
746 BID_RETURN (res);
747 }
748 // else cases that can be rounded to a 32-bit int fall through
749 // to '1 <= q + exp <= 10'
750 }
751 }
752 }
753 // n is not too large to be converted to int32: -2^31 <= n < 2^31
754 // Note: some of the cases tested for above fall through to this point
755 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
756 // return -1 or 0
757 if (x_sign)
758 res = 0xffffffff;
759 else
760 res = 0x00000000;
761 BID_RETURN (res);
762 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
763 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
764 // to nearest to a 32-bit signed integer
765 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
766 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
767 // chop off ind digits from the lower part of C1
768 // C1 fits in 64 bits
769 // calculate C* and f*
770 // C* is actually floor(C*) in this case
771 // C* and f* need shifting and masking, as shown by
772 // shiftright128[] and maskhigh128[]
773 // 1 <= x <= 15
774 // kx = 10^(-x) = ten2mk64[ind - 1]
775 // C* = C1 * 10^(-x)
776 // the approximation of 10^(-x) was rounded up to 54 bits
777 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
778 Cstar = P128.w[1];
779 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
780 fstar.w[0] = P128.w[0];
781 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
782 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
783 // C* = floor(C*) (logical right shift; C has p decimal digits,
784 // correct by Property 1)
785 // n = C* * 10^(e+x)
786
787 // shift right C* by Ex-64 = shiftright128[ind]
788 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
789 Cstar = Cstar >> shift;
790 // determine inexactness of the rounding of C*
791 // if (0 < f* < 10^(-x)) then
792 // the result is exact
793 // else // if (f* > T*) then
794 // the result is inexact
795 if (ind - 1 <= 2) {
796 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
797 // ten2mk128trunc[ind -1].w[1] is identical to
798 // ten2mk128[ind -1].w[1]
799 if (x_sign) { // negative and inexact
800 Cstar++;
801 }
802 } // else the result is exact
803 } else { // if 3 <= ind - 1 <= 14
804 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
805 // ten2mk128trunc[ind -1].w[1] is identical to
806 // ten2mk128[ind -1].w[1]
807 if (x_sign) { // negative and inexact
808 Cstar++;
809 }
810 } // else the result is exact
811 }
812
813 if (x_sign)
814 res = -Cstar;
815 else
816 res = Cstar;
817 } else if (exp == 0) {
818 // 1 <= q <= 10
819 // res = +/-C (exact)
820 if (x_sign)
821 res = -C1;
822 else
823 res = C1;
824 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
825 // res = +/-C * 10^exp (exact)
826 if (x_sign)
827 res = -C1 * ten2k64[exp];
828 else
829 res = C1 * ten2k64[exp];
830 }
831 }
832 BID_RETURN (res);
833 }
834
835 /*****************************************************************************
836 * BID64_to_int32_xfloor
837 ****************************************************************************/
838
839 #if DECIMAL_CALL_BY_REFERENCE
840 void
841 bid64_to_int32_xfloor (int *pres,
842 UINT64 *
843 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
844 _EXC_INFO_PARAM) {
845 UINT64 x = *px;
846 #else
847 int
848 bid64_to_int32_xfloor (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
849 _EXC_INFO_PARAM) {
850 #endif
851 int res;
852 UINT64 x_sign;
853 UINT64 x_exp;
854 int exp; // unbiased exponent
855 // Note: C1 represents x_significand (UINT64)
856 UINT64 tmp64;
857 BID_UI64DOUBLE tmp1;
858 unsigned int x_nr_bits;
859 int q, ind, shift;
860 UINT64 C1;
861 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
862 UINT128 fstar;
863 UINT128 P128;
864
865 // check for NaN or Infinity
866 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
867 // set invalid flag
868 *pfpsf |= INVALID_EXCEPTION;
869 // return Integer Indefinite
870 res = 0x80000000;
871 BID_RETURN (res);
872 }
873 // unpack x
874 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
875 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
876 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
877 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
878 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
879 if (C1 > 9999999999999999ull) { // non-canonical
880 x_exp = 0;
881 C1 = 0;
882 }
883 } else {
884 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
885 C1 = x & MASK_BINARY_SIG1;
886 }
887
888 // check for zeros (possibly from non-canonical values)
889 if (C1 == 0x0ull) {
890 // x is 0
891 res = 0x00000000;
892 BID_RETURN (res);
893 }
894 // x is not special and is not zero
895
896 // q = nr. of decimal digits in x (1 <= q <= 54)
897 // determine first the nr. of bits in x
898 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
899 // split the 64-bit value in two 32-bit halves to avoid rounding errors
900 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
901 tmp1.d = (double) (C1 >> 32); // exact conversion
902 x_nr_bits =
903 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
904 } else { // x < 2^32
905 tmp1.d = (double) C1; // exact conversion
906 x_nr_bits =
907 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
908 }
909 } else { // if x < 2^53
910 tmp1.d = (double) C1; // exact conversion
911 x_nr_bits =
912 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
913 }
914 q = nr_digits[x_nr_bits - 1].digits;
915 if (q == 0) {
916 q = nr_digits[x_nr_bits - 1].digits1;
917 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
918 q++;
919 }
920 exp = x_exp - 398; // unbiased exponent
921
922 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
923 // set invalid flag
924 *pfpsf |= INVALID_EXCEPTION;
925 // return Integer Indefinite
926 res = 0x80000000;
927 BID_RETURN (res);
928 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
929 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
930 // so x rounded to an integer may or may not fit in a signed 32-bit int
931 // the cases that do not fit are identified here; the ones that fit
932 // fall through and will be handled with other cases further,
933 // under '1 <= q + exp <= 10'
934 if (x_sign) { // if n < 0 and q + exp = 10
935 // if n < -2^31 then n is too large
936 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
937 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=16
938 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
939 if (q <= 11) {
940 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
941 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
942 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
943 if (tmp64 > 0x500000000ull) {
944 // set invalid flag
945 *pfpsf |= INVALID_EXCEPTION;
946 // return Integer Indefinite
947 res = 0x80000000;
948 BID_RETURN (res);
949 }
950 // else cases that can be rounded to a 32-bit int fall through
951 // to '1 <= q + exp <= 10'
952 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
953 // C * 10^(11-q) > 0x500000000 <=>
954 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
955 // (scale 2^31+1 up)
956 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
957 tmp64 = 0x500000000ull * ten2k64[q - 11];
958 if (C1 > tmp64) {
959 // set invalid flag
960 *pfpsf |= INVALID_EXCEPTION;
961 // return Integer Indefinite
962 res = 0x80000000;
963 BID_RETURN (res);
964 }
965 // else cases that can be rounded to a 32-bit int fall through
966 // to '1 <= q + exp <= 10'
967 }
968 } else { // if n > 0 and q + exp = 10
969 // if n >= 2^31 then n is too large
970 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
971 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
972 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
973 if (q <= 11) {
974 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
975 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
976 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
977 if (tmp64 >= 0x500000000ull) {
978 // set invalid flag
979 *pfpsf |= INVALID_EXCEPTION;
980 // return Integer Indefinite
981 res = 0x80000000;
982 BID_RETURN (res);
983 }
984 // else cases that can be rounded to a 32-bit int fall through
985 // to '1 <= q + exp <= 10'
986 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
987 // C * 10^(11-q) >= 0x500000000 <=>
988 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
989 // (scale 2^31-1 up)
990 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
991 tmp64 = 0x500000000ull * ten2k64[q - 11];
992 if (C1 >= tmp64) {
993 // set invalid flag
994 *pfpsf |= INVALID_EXCEPTION;
995 // return Integer Indefinite
996 res = 0x80000000;
997 BID_RETURN (res);
998 }
999 // else cases that can be rounded to a 32-bit int fall through
1000 // to '1 <= q + exp <= 10'
1001 }
1002 }
1003 }
1004 // n is not too large to be converted to int32: -2^31 <= n < 2^31
1005 // Note: some of the cases tested for above fall through to this point
1006 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1007 // set inexact flag
1008 *pfpsf |= INEXACT_EXCEPTION;
1009 // return -1 or 0
1010 if (x_sign)
1011 res = 0xffffffff;
1012 else
1013 res = 0x00000000;
1014 BID_RETURN (res);
1015 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1016 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
1017 // to nearest to a 32-bit signed integer
1018 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1019 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1020 // chop off ind digits from the lower part of C1
1021 // C1 fits in 64 bits
1022 // calculate C* and f*
1023 // C* is actually floor(C*) in this case
1024 // C* and f* need shifting and masking, as shown by
1025 // shiftright128[] and maskhigh128[]
1026 // 1 <= x <= 15
1027 // kx = 10^(-x) = ten2mk64[ind - 1]
1028 // C* = C1 * 10^(-x)
1029 // the approximation of 10^(-x) was rounded up to 54 bits
1030 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1031 Cstar = P128.w[1];
1032 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1033 fstar.w[0] = P128.w[0];
1034 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1035 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1036 // C* = floor(C*) (logical right shift; C has p decimal digits,
1037 // correct by Property 1)
1038 // n = C* * 10^(e+x)
1039
1040 // shift right C* by Ex-64 = shiftright128[ind]
1041 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1042 Cstar = Cstar >> shift;
1043 // determine inexactness of the rounding of C*
1044 // if (0 < f* < 10^(-x)) then
1045 // the result is exact
1046 // else // if (f* > T*) then
1047 // the result is inexact
1048 if (ind - 1 <= 2) {
1049 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1050 // ten2mk128trunc[ind -1].w[1] is identical to
1051 // ten2mk128[ind -1].w[1]
1052 if (x_sign) { // negative and inexact
1053 Cstar++;
1054 }
1055 // set the inexact flag
1056 *pfpsf |= INEXACT_EXCEPTION;
1057 } // else the result is exact
1058 } else { // if 3 <= ind - 1 <= 14
1059 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1060 // ten2mk128trunc[ind -1].w[1] is identical to
1061 // ten2mk128[ind -1].w[1]
1062 if (x_sign) { // negative and inexact
1063 Cstar++;
1064 }
1065 // set the inexact flag
1066 *pfpsf |= INEXACT_EXCEPTION;
1067 } // else the result is exact
1068 }
1069
1070 if (x_sign)
1071 res = -Cstar;
1072 else
1073 res = Cstar;
1074 } else if (exp == 0) {
1075 // 1 <= q <= 10
1076 // res = +/-C (exact)
1077 if (x_sign)
1078 res = -C1;
1079 else
1080 res = C1;
1081 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1082 // res = +/-C * 10^exp (exact)
1083 if (x_sign)
1084 res = -C1 * ten2k64[exp];
1085 else
1086 res = C1 * ten2k64[exp];
1087 }
1088 }
1089 BID_RETURN (res);
1090 }
1091
1092 /*****************************************************************************
1093 * BID64_to_int32_ceil
1094 ****************************************************************************/
1095
1096 #if DECIMAL_CALL_BY_REFERENCE
1097 void
1098 bid64_to_int32_ceil (int *pres,
1099 UINT64 *
1100 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1101 _EXC_INFO_PARAM) {
1102 UINT64 x = *px;
1103 #else
1104 int
1105 bid64_to_int32_ceil (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1106 _EXC_INFO_PARAM) {
1107 #endif
1108 int res;
1109 UINT64 x_sign;
1110 UINT64 x_exp;
1111 int exp; // unbiased exponent
1112 // Note: C1 represents x_significand (UINT64)
1113 UINT64 tmp64;
1114 BID_UI64DOUBLE tmp1;
1115 unsigned int x_nr_bits;
1116 int q, ind, shift;
1117 UINT64 C1;
1118 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1119 UINT128 fstar;
1120 UINT128 P128;
1121
1122 // check for NaN or Infinity
1123 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1124 // set invalid flag
1125 *pfpsf |= INVALID_EXCEPTION;
1126 // return Integer Indefinite
1127 res = 0x80000000;
1128 BID_RETURN (res);
1129 }
1130 // unpack x
1131 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1132 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1133 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1134 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1135 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1136 if (C1 > 9999999999999999ull) { // non-canonical
1137 x_exp = 0;
1138 C1 = 0;
1139 }
1140 } else {
1141 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1142 C1 = x & MASK_BINARY_SIG1;
1143 }
1144
1145 // check for zeros (possibly from non-canonical values)
1146 if (C1 == 0x0ull) {
1147 // x is 0
1148 res = 0x00000000;
1149 BID_RETURN (res);
1150 }
1151 // x is not special and is not zero
1152
1153 // q = nr. of decimal digits in x (1 <= q <= 54)
1154 // determine first the nr. of bits in x
1155 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1156 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1157 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1158 tmp1.d = (double) (C1 >> 32); // exact conversion
1159 x_nr_bits =
1160 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1161 } else { // x < 2^32
1162 tmp1.d = (double) C1; // exact conversion
1163 x_nr_bits =
1164 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1165 }
1166 } else { // if x < 2^53
1167 tmp1.d = (double) C1; // exact conversion
1168 x_nr_bits =
1169 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1170 }
1171 q = nr_digits[x_nr_bits - 1].digits;
1172 if (q == 0) {
1173 q = nr_digits[x_nr_bits - 1].digits1;
1174 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1175 q++;
1176 }
1177 exp = x_exp - 398; // unbiased exponent
1178
1179 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1180 // set invalid flag
1181 *pfpsf |= INVALID_EXCEPTION;
1182 // return Integer Indefinite
1183 res = 0x80000000;
1184 BID_RETURN (res);
1185 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1186 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1187 // so x rounded to an integer may or may not fit in a signed 32-bit int
1188 // the cases that do not fit are identified here; the ones that fit
1189 // fall through and will be handled with other cases further,
1190 // under '1 <= q + exp <= 10'
1191 if (x_sign) { // if n < 0 and q + exp = 10
1192 // if n <= -2^31 - 1 then n is too large
1193 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1194 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1195 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1196 if (q <= 11) {
1197 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1198 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1199 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1200 if (tmp64 >= 0x50000000aull) {
1201 // set invalid flag
1202 *pfpsf |= INVALID_EXCEPTION;
1203 // return Integer Indefinite
1204 res = 0x80000000;
1205 BID_RETURN (res);
1206 }
1207 // else cases that can be rounded to a 32-bit int fall through
1208 // to '1 <= q + exp <= 10'
1209 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1210 // C * 10^(11-q) >= 0x50000000a <=>
1211 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1212 // (scale 2^31+1 up)
1213 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1214 tmp64 = 0x50000000aull * ten2k64[q - 11];
1215 if (C1 >= tmp64) {
1216 // set invalid flag
1217 *pfpsf |= INVALID_EXCEPTION;
1218 // return Integer Indefinite
1219 res = 0x80000000;
1220 BID_RETURN (res);
1221 }
1222 // else cases that can be rounded to a 32-bit int fall through
1223 // to '1 <= q + exp <= 10'
1224 }
1225 } else { // if n > 0 and q + exp = 10
1226 // if n > 2^31 - 1 then n is too large
1227 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
1228 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=16
1229 // <=> C * 10^(11-q) > 0x4fffffff6, 1<=q<=16
1230 if (q <= 11) {
1231 // Note: C * 10^(11-q) has 10 or 11 digits; 0x4fffffff6 has 11 digits
1232 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1233 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1234 if (tmp64 > 0x4fffffff6ull) {
1235 // set invalid flag
1236 *pfpsf |= INVALID_EXCEPTION;
1237 // return Integer Indefinite
1238 res = 0x80000000;
1239 BID_RETURN (res);
1240 }
1241 // else cases that can be rounded to a 32-bit int fall through
1242 // to '1 <= q + exp <= 10'
1243 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1244 // C * 10^(11-q) > 0x4fffffff6 <=>
1245 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1246 // (scale 2^31-1 up)
1247 // Note: 0x4fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1248 tmp64 = 0x4fffffff6ull * ten2k64[q - 11];
1249 if (C1 > tmp64) {
1250 // set invalid flag
1251 *pfpsf |= INVALID_EXCEPTION;
1252 // return Integer Indefinite
1253 res = 0x80000000;
1254 BID_RETURN (res);
1255 }
1256 // else cases that can be rounded to a 32-bit int fall through
1257 // to '1 <= q + exp <= 10'
1258 }
1259 }
1260 }
1261 // n is not too large to be converted to int32: -2^31 - 1 < n <= 2^31 - 1
1262 // Note: some of the cases tested for above fall through to this point
1263 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1264 // return 0 or 1
1265 if (x_sign)
1266 res = 0x00000000;
1267 else
1268 res = 0x00000001;
1269 BID_RETURN (res);
1270 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1271 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
1272 // to nearest to a 32-bit signed integer
1273 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1274 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1275 // chop off ind digits from the lower part of C1
1276 // C1 fits in 64 bits
1277 // calculate C* and f*
1278 // C* is actually floor(C*) in this case
1279 // C* and f* need shifting and masking, as shown by
1280 // shiftright128[] and maskhigh128[]
1281 // 1 <= x <= 15
1282 // kx = 10^(-x) = ten2mk64[ind - 1]
1283 // C* = C1 * 10^(-x)
1284 // the approximation of 10^(-x) was rounded up to 54 bits
1285 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1286 Cstar = P128.w[1];
1287 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1288 fstar.w[0] = P128.w[0];
1289 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1290 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1291 // C* = floor(C*) (logical right shift; C has p decimal digits,
1292 // correct by Property 1)
1293 // n = C* * 10^(e+x)
1294
1295 // shift right C* by Ex-64 = shiftright128[ind]
1296 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1297 Cstar = Cstar >> shift;
1298 // determine inexactness of the rounding of C*
1299 // if (0 < f* < 10^(-x)) then
1300 // the result is exact
1301 // else // if (f* > T*) then
1302 // the result is inexact
1303 if (ind - 1 <= 2) {
1304 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1305 // ten2mk128trunc[ind -1].w[1] is identical to
1306 // ten2mk128[ind -1].w[1]
1307 if (!x_sign) { // positive and inexact
1308 Cstar++;
1309 }
1310 } // else the result is exact
1311 } else { // if 3 <= ind - 1 <= 14
1312 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1313 // ten2mk128trunc[ind -1].w[1] is identical to
1314 // ten2mk128[ind -1].w[1]
1315 if (!x_sign) { // positive and inexact
1316 Cstar++;
1317 }
1318 } // else the result is exact
1319 }
1320
1321 if (x_sign)
1322 res = -Cstar;
1323 else
1324 res = Cstar;
1325 } else if (exp == 0) {
1326 // 1 <= q <= 10
1327 // res = +/-C (exact)
1328 if (x_sign)
1329 res = -C1;
1330 else
1331 res = C1;
1332 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1333 // res = +/-C * 10^exp (exact)
1334 if (x_sign)
1335 res = -C1 * ten2k64[exp];
1336 else
1337 res = C1 * ten2k64[exp];
1338 }
1339 }
1340 BID_RETURN (res);
1341 }
1342
1343 /*****************************************************************************
1344 * BID64_to_int32_xceil
1345 ****************************************************************************/
1346
1347 #if DECIMAL_CALL_BY_REFERENCE
1348 void
1349 bid64_to_int32_xceil (int *pres,
1350 UINT64 *
1351 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1352 _EXC_INFO_PARAM) {
1353 UINT64 x = *px;
1354 #else
1355 int
1356 bid64_to_int32_xceil (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1357 _EXC_INFO_PARAM) {
1358 #endif
1359 int res;
1360 UINT64 x_sign;
1361 UINT64 x_exp;
1362 int exp; // unbiased exponent
1363 // Note: C1 represents x_significand (UINT64)
1364 UINT64 tmp64;
1365 BID_UI64DOUBLE tmp1;
1366 unsigned int x_nr_bits;
1367 int q, ind, shift;
1368 UINT64 C1;
1369 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1370 UINT128 fstar;
1371 UINT128 P128;
1372
1373 // check for NaN or Infinity
1374 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1375 // set invalid flag
1376 *pfpsf |= INVALID_EXCEPTION;
1377 // return Integer Indefinite
1378 res = 0x80000000;
1379 BID_RETURN (res);
1380 }
1381 // unpack x
1382 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1383 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1384 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1385 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1386 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1387 if (C1 > 9999999999999999ull) { // non-canonical
1388 x_exp = 0;
1389 C1 = 0;
1390 }
1391 } else {
1392 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1393 C1 = x & MASK_BINARY_SIG1;
1394 }
1395
1396 // check for zeros (possibly from non-canonical values)
1397 if (C1 == 0x0ull) {
1398 // x is 0
1399 res = 0x00000000;
1400 BID_RETURN (res);
1401 }
1402 // x is not special and is not zero
1403
1404 // q = nr. of decimal digits in x (1 <= q <= 54)
1405 // determine first the nr. of bits in x
1406 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1407 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1408 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1409 tmp1.d = (double) (C1 >> 32); // exact conversion
1410 x_nr_bits =
1411 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1412 } else { // x < 2^32
1413 tmp1.d = (double) C1; // exact conversion
1414 x_nr_bits =
1415 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1416 }
1417 } else { // if x < 2^53
1418 tmp1.d = (double) C1; // exact conversion
1419 x_nr_bits =
1420 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1421 }
1422 q = nr_digits[x_nr_bits - 1].digits;
1423 if (q == 0) {
1424 q = nr_digits[x_nr_bits - 1].digits1;
1425 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1426 q++;
1427 }
1428 exp = x_exp - 398; // unbiased exponent
1429
1430 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1431 // set invalid flag
1432 *pfpsf |= INVALID_EXCEPTION;
1433 // return Integer Indefinite
1434 res = 0x80000000;
1435 BID_RETURN (res);
1436 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1437 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1438 // so x rounded to an integer may or may not fit in a signed 32-bit int
1439 // the cases that do not fit are identified here; the ones that fit
1440 // fall through and will be handled with other cases further,
1441 // under '1 <= q + exp <= 10'
1442 if (x_sign) { // if n < 0 and q + exp = 10
1443 // if n <= -2^31 - 1 then n is too large
1444 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1445 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1446 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1447 if (q <= 11) {
1448 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1449 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1450 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1451 if (tmp64 >= 0x50000000aull) {
1452 // set invalid flag
1453 *pfpsf |= INVALID_EXCEPTION;
1454 // return Integer Indefinite
1455 res = 0x80000000;
1456 BID_RETURN (res);
1457 }
1458 // else cases that can be rounded to a 32-bit int fall through
1459 // to '1 <= q + exp <= 10'
1460 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1461 // C * 10^(11-q) >= 0x50000000a <=>
1462 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1463 // (scale 2^31+1 up)
1464 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1465 tmp64 = 0x50000000aull * ten2k64[q - 11];
1466 if (C1 >= tmp64) {
1467 // set invalid flag
1468 *pfpsf |= INVALID_EXCEPTION;
1469 // return Integer Indefinite
1470 res = 0x80000000;
1471 BID_RETURN (res);
1472 }
1473 // else cases that can be rounded to a 32-bit int fall through
1474 // to '1 <= q + exp <= 10'
1475 }
1476 } else { // if n > 0 and q + exp = 10
1477 // if n > 2^31 - 1 then n is too large
1478 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
1479 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=16
1480 // <=> C * 10^(11-q) > 0x4fffffff6, 1<=q<=16
1481 if (q <= 11) {
1482 // Note: C * 10^(11-q) has 10 or 11 digits; 0x4fffffff6 has 11 digits
1483 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1484 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1485 if (tmp64 > 0x4fffffff6ull) {
1486 // set invalid flag
1487 *pfpsf |= INVALID_EXCEPTION;
1488 // return Integer Indefinite
1489 res = 0x80000000;
1490 BID_RETURN (res);
1491 }
1492 // else cases that can be rounded to a 32-bit int fall through
1493 // to '1 <= q + exp <= 10'
1494 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1495 // C * 10^(11-q) > 0x4fffffff6 <=>
1496 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1497 // (scale 2^31-1 up)
1498 // Note: 0x4fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1499 tmp64 = 0x4fffffff6ull * ten2k64[q - 11];
1500 if (C1 > tmp64) {
1501 // set invalid flag
1502 *pfpsf |= INVALID_EXCEPTION;
1503 // return Integer Indefinite
1504 res = 0x80000000;
1505 BID_RETURN (res);
1506 }
1507 // else cases that can be rounded to a 32-bit int fall through
1508 // to '1 <= q + exp <= 10'
1509 }
1510 }
1511 }
1512 // n is not too large to be converted to int32: -2^31 - 1 < n <= 2^31 - 1
1513 // Note: some of the cases tested for above fall through to this point
1514 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1515 // set inexact flag
1516 *pfpsf |= INEXACT_EXCEPTION;
1517 // return 0 or 1
1518 if (x_sign)
1519 res = 0x00000000;
1520 else
1521 res = 0x00000001;
1522 BID_RETURN (res);
1523 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1524 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
1525 // to nearest to a 32-bit signed integer
1526 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1527 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1528 // chop off ind digits from the lower part of C1
1529 // C1 fits in 64 bits
1530 // calculate C* and f*
1531 // C* is actually floor(C*) in this case
1532 // C* and f* need shifting and masking, as shown by
1533 // shiftright128[] and maskhigh128[]
1534 // 1 <= x <= 15
1535 // kx = 10^(-x) = ten2mk64[ind - 1]
1536 // C* = C1 * 10^(-x)
1537 // the approximation of 10^(-x) was rounded up to 54 bits
1538 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1539 Cstar = P128.w[1];
1540 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1541 fstar.w[0] = P128.w[0];
1542 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1543 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1544 // C* = floor(C*) (logical right shift; C has p decimal digits,
1545 // correct by Property 1)
1546 // n = C* * 10^(e+x)
1547
1548 // shift right C* by Ex-64 = shiftright128[ind]
1549 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1550 Cstar = Cstar >> shift;
1551 // determine inexactness of the rounding of C*
1552 // if (0 < f* < 10^(-x)) then
1553 // the result is exact
1554 // else // if (f* > T*) then
1555 // the result is inexact
1556 if (ind - 1 <= 2) {
1557 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1558 // ten2mk128trunc[ind -1].w[1] is identical to
1559 // ten2mk128[ind -1].w[1]
1560 if (!x_sign) { // positive and inexact
1561 Cstar++;
1562 }
1563 // set the inexact flag
1564 *pfpsf |= INEXACT_EXCEPTION;
1565 } // else the result is exact
1566 } else { // if 3 <= ind - 1 <= 14
1567 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1568 // ten2mk128trunc[ind -1].w[1] is identical to
1569 // ten2mk128[ind -1].w[1]
1570 if (!x_sign) { // positive and inexact
1571 Cstar++;
1572 }
1573 // set the inexact flag
1574 *pfpsf |= INEXACT_EXCEPTION;
1575 } // else the result is exact
1576 }
1577
1578 if (x_sign)
1579 res = -Cstar;
1580 else
1581 res = Cstar;
1582 } else if (exp == 0) {
1583 // 1 <= q <= 10
1584 // res = +/-C (exact)
1585 if (x_sign)
1586 res = -C1;
1587 else
1588 res = C1;
1589 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1590 // res = +/-C * 10^exp (exact)
1591 if (x_sign)
1592 res = -C1 * ten2k64[exp];
1593 else
1594 res = C1 * ten2k64[exp];
1595 }
1596 }
1597 BID_RETURN (res);
1598 }
1599
1600 /*****************************************************************************
1601 * BID64_to_int32_int
1602 ****************************************************************************/
1603
1604 #if DECIMAL_CALL_BY_REFERENCE
1605 void
1606 bid64_to_int32_int (int *pres,
1607 UINT64 *
1608 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1609 _EXC_INFO_PARAM) {
1610 UINT64 x = *px;
1611 #else
1612 int
1613 bid64_to_int32_int (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1614 _EXC_INFO_PARAM) {
1615 #endif
1616 int res;
1617 UINT64 x_sign;
1618 UINT64 x_exp;
1619 int exp; // unbiased exponent
1620 // Note: C1 represents x_significand (UINT64)
1621 UINT64 tmp64;
1622 BID_UI64DOUBLE tmp1;
1623 unsigned int x_nr_bits;
1624 int q, ind, shift;
1625 UINT64 C1;
1626 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1627 UINT128 P128;
1628
1629 // check for NaN or Infinity
1630 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1631 // set invalid flag
1632 *pfpsf |= INVALID_EXCEPTION;
1633 // return Integer Indefinite
1634 res = 0x80000000;
1635 BID_RETURN (res);
1636 }
1637 // unpack x
1638 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1639 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1640 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1641 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1642 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1643 if (C1 > 9999999999999999ull) { // non-canonical
1644 x_exp = 0;
1645 C1 = 0;
1646 }
1647 } else {
1648 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1649 C1 = x & MASK_BINARY_SIG1;
1650 }
1651
1652 // check for zeros (possibly from non-canonical values)
1653 if (C1 == 0x0ull) {
1654 // x is 0
1655 res = 0x00000000;
1656 BID_RETURN (res);
1657 }
1658 // x is not special and is not zero
1659
1660 // q = nr. of decimal digits in x (1 <= q <= 54)
1661 // determine first the nr. of bits in x
1662 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1663 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1664 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1665 tmp1.d = (double) (C1 >> 32); // exact conversion
1666 x_nr_bits =
1667 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1668 } else { // x < 2^32
1669 tmp1.d = (double) C1; // exact conversion
1670 x_nr_bits =
1671 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1672 }
1673 } else { // if x < 2^53
1674 tmp1.d = (double) C1; // exact conversion
1675 x_nr_bits =
1676 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1677 }
1678 q = nr_digits[x_nr_bits - 1].digits;
1679 if (q == 0) {
1680 q = nr_digits[x_nr_bits - 1].digits1;
1681 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1682 q++;
1683 }
1684 exp = x_exp - 398; // unbiased exponent
1685
1686 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1687 // set invalid flag
1688 *pfpsf |= INVALID_EXCEPTION;
1689 // return Integer Indefinite
1690 res = 0x80000000;
1691 BID_RETURN (res);
1692 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1693 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1694 // so x rounded to an integer may or may not fit in a signed 32-bit int
1695 // the cases that do not fit are identified here; the ones that fit
1696 // fall through and will be handled with other cases further,
1697 // under '1 <= q + exp <= 10'
1698 if (x_sign) { // if n < 0 and q + exp = 10
1699 // if n <= -2^31 - 1 then n is too large
1700 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1701 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1702 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1703 if (q <= 11) {
1704 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1705 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1706 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1707 if (tmp64 >= 0x50000000aull) {
1708 // set invalid flag
1709 *pfpsf |= INVALID_EXCEPTION;
1710 // return Integer Indefinite
1711 res = 0x80000000;
1712 BID_RETURN (res);
1713 }
1714 // else cases that can be rounded to a 32-bit int fall through
1715 // to '1 <= q + exp <= 10'
1716 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1717 // C * 10^(11-q) >= 0x50000000a <=>
1718 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1719 // (scale 2^31+1 up)
1720 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1721 tmp64 = 0x50000000aull * ten2k64[q - 11];
1722 if (C1 >= tmp64) {
1723 // set invalid flag
1724 *pfpsf |= INVALID_EXCEPTION;
1725 // return Integer Indefinite
1726 res = 0x80000000;
1727 BID_RETURN (res);
1728 }
1729 // else cases that can be rounded to a 32-bit int fall through
1730 // to '1 <= q + exp <= 10'
1731 }
1732 } else { // if n > 0 and q + exp = 10
1733 // if n >= 2^31 then n is too large
1734 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
1735 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
1736 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
1737 if (q <= 11) {
1738 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
1739 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1740 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1741 if (tmp64 >= 0x500000000ull) {
1742 // set invalid flag
1743 *pfpsf |= INVALID_EXCEPTION;
1744 // return Integer Indefinite
1745 res = 0x80000000;
1746 BID_RETURN (res);
1747 }
1748 // else cases that can be rounded to a 32-bit int fall through
1749 // to '1 <= q + exp <= 10'
1750 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1751 // C * 10^(11-q) >= 0x500000000 <=>
1752 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
1753 // (scale 2^31-1 up)
1754 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
1755 tmp64 = 0x500000000ull * ten2k64[q - 11];
1756 if (C1 >= tmp64) {
1757 // set invalid flag
1758 *pfpsf |= INVALID_EXCEPTION;
1759 // return Integer Indefinite
1760 res = 0x80000000;
1761 BID_RETURN (res);
1762 }
1763 // else cases that can be rounded to a 32-bit int fall through
1764 // to '1 <= q + exp <= 10'
1765 }
1766 }
1767 }
1768 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
1769 // Note: some of the cases tested for above fall through to this point
1770 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1771 // return 0
1772 res = 0x00000000;
1773 BID_RETURN (res);
1774 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1775 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
1776 // to nearest to a 32-bit signed integer
1777 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1778 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1779 // chop off ind digits from the lower part of C1
1780 // C1 fits in 64 bits
1781 // calculate C* and f*
1782 // C* is actually floor(C*) in this case
1783 // C* and f* need shifting and masking, as shown by
1784 // shiftright128[] and maskhigh128[]
1785 // 1 <= x <= 15
1786 // kx = 10^(-x) = ten2mk64[ind - 1]
1787 // C* = C1 * 10^(-x)
1788 // the approximation of 10^(-x) was rounded up to 54 bits
1789 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1790 Cstar = P128.w[1];
1791 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1792 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1793 // C* = floor(C*) (logical right shift; C has p decimal digits,
1794 // correct by Property 1)
1795 // n = C* * 10^(e+x)
1796
1797 // shift right C* by Ex-64 = shiftright128[ind]
1798 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1799 Cstar = Cstar >> shift;
1800 if (x_sign)
1801 res = -Cstar;
1802 else
1803 res = Cstar;
1804 } else if (exp == 0) {
1805 // 1 <= q <= 10
1806 // res = +/-C (exact)
1807 if (x_sign)
1808 res = -C1;
1809 else
1810 res = C1;
1811 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1812 // res = +/-C * 10^exp (exact)
1813 if (x_sign)
1814 res = -C1 * ten2k64[exp];
1815 else
1816 res = C1 * ten2k64[exp];
1817 }
1818 }
1819 BID_RETURN (res);
1820 }
1821
1822 /*****************************************************************************
1823 * BID64_to_int32_xint
1824 ****************************************************************************/
1825
1826 #if DECIMAL_CALL_BY_REFERENCE
1827 void
1828 bid64_to_int32_xint (int *pres,
1829 UINT64 *
1830 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1831 _EXC_INFO_PARAM) {
1832 UINT64 x = *px;
1833 #else
1834 int
1835 bid64_to_int32_xint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1836 _EXC_INFO_PARAM) {
1837 #endif
1838 int res;
1839 UINT64 x_sign;
1840 UINT64 x_exp;
1841 int exp; // unbiased exponent
1842 // Note: C1 represents x_significand (UINT64)
1843 UINT64 tmp64;
1844 BID_UI64DOUBLE tmp1;
1845 unsigned int x_nr_bits;
1846 int q, ind, shift;
1847 UINT64 C1;
1848 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1849 UINT128 fstar;
1850 UINT128 P128;
1851
1852 // check for NaN or Infinity
1853 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1854 // set invalid flag
1855 *pfpsf |= INVALID_EXCEPTION;
1856 // return Integer Indefinite
1857 res = 0x80000000;
1858 BID_RETURN (res);
1859 }
1860 // unpack x
1861 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1862 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1863 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1864 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1865 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1866 if (C1 > 9999999999999999ull) { // non-canonical
1867 x_exp = 0;
1868 C1 = 0;
1869 }
1870 } else {
1871 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1872 C1 = x & MASK_BINARY_SIG1;
1873 }
1874
1875 // check for zeros (possibly from non-canonical values)
1876 if (C1 == 0x0ull) {
1877 // x is 0
1878 res = 0x00000000;
1879 BID_RETURN (res);
1880 }
1881 // x is not special and is not zero
1882
1883 // q = nr. of decimal digits in x (1 <= q <= 54)
1884 // determine first the nr. of bits in x
1885 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1886 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1887 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1888 tmp1.d = (double) (C1 >> 32); // exact conversion
1889 x_nr_bits =
1890 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1891 } else { // x < 2^32
1892 tmp1.d = (double) C1; // exact conversion
1893 x_nr_bits =
1894 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1895 }
1896 } else { // if x < 2^53
1897 tmp1.d = (double) C1; // exact conversion
1898 x_nr_bits =
1899 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1900 }
1901 q = nr_digits[x_nr_bits - 1].digits;
1902 if (q == 0) {
1903 q = nr_digits[x_nr_bits - 1].digits1;
1904 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1905 q++;
1906 }
1907 exp = x_exp - 398; // unbiased exponent
1908
1909 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1910 // set invalid flag
1911 *pfpsf |= INVALID_EXCEPTION;
1912 // return Integer Indefinite
1913 res = 0x80000000;
1914 BID_RETURN (res);
1915 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1916 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1917 // so x rounded to an integer may or may not fit in a signed 32-bit int
1918 // the cases that do not fit are identified here; the ones that fit
1919 // fall through and will be handled with other cases further,
1920 // under '1 <= q + exp <= 10'
1921 if (x_sign) { // if n < 0 and q + exp = 10
1922 // if n <= -2^31 - 1 then n is too large
1923 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1924 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1925 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1926 if (q <= 11) {
1927 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1928 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1929 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1930 if (tmp64 >= 0x50000000aull) {
1931 // set invalid flag
1932 *pfpsf |= INVALID_EXCEPTION;
1933 // return Integer Indefinite
1934 res = 0x80000000;
1935 BID_RETURN (res);
1936 }
1937 // else cases that can be rounded to a 32-bit int fall through
1938 // to '1 <= q + exp <= 10'
1939 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1940 // C * 10^(11-q) >= 0x50000000a <=>
1941 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1942 // (scale 2^31+1 up)
1943 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1944 tmp64 = 0x50000000aull * ten2k64[q - 11];
1945 if (C1 >= tmp64) {
1946 // set invalid flag
1947 *pfpsf |= INVALID_EXCEPTION;
1948 // return Integer Indefinite
1949 res = 0x80000000;
1950 BID_RETURN (res);
1951 }
1952 // else cases that can be rounded to a 32-bit int fall through
1953 // to '1 <= q + exp <= 10'
1954 }
1955 } else { // if n > 0 and q + exp = 10
1956 // if n >= 2^31 then n is too large
1957 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
1958 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
1959 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
1960 if (q <= 11) {
1961 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
1962 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1963 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1964 if (tmp64 >= 0x500000000ull) {
1965 // set invalid flag
1966 *pfpsf |= INVALID_EXCEPTION;
1967 // return Integer Indefinite
1968 res = 0x80000000;
1969 BID_RETURN (res);
1970 }
1971 // else cases that can be rounded to a 32-bit int fall through
1972 // to '1 <= q + exp <= 10'
1973 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1974 // C * 10^(11-q) >= 0x500000000 <=>
1975 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
1976 // (scale 2^31-1 up)
1977 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
1978 tmp64 = 0x500000000ull * ten2k64[q - 11];
1979 if (C1 >= tmp64) {
1980 // set invalid flag
1981 *pfpsf |= INVALID_EXCEPTION;
1982 // return Integer Indefinite
1983 res = 0x80000000;
1984 BID_RETURN (res);
1985 }
1986 // else cases that can be rounded to a 32-bit int fall through
1987 // to '1 <= q + exp <= 10'
1988 }
1989 }
1990 }
1991 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
1992 // Note: some of the cases tested for above fall through to this point
1993 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1994 // set inexact flag
1995 *pfpsf |= INEXACT_EXCEPTION;
1996 // return 0
1997 res = 0x00000000;
1998 BID_RETURN (res);
1999 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2000 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
2001 // to nearest to a 32-bit signed integer
2002 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2003 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2004 // chop off ind digits from the lower part of C1
2005 // C1 fits in 64 bits
2006 // calculate C* and f*
2007 // C* is actually floor(C*) in this case
2008 // C* and f* need shifting and masking, as shown by
2009 // shiftright128[] and maskhigh128[]
2010 // 1 <= x <= 15
2011 // kx = 10^(-x) = ten2mk64[ind - 1]
2012 // C* = C1 * 10^(-x)
2013 // the approximation of 10^(-x) was rounded up to 54 bits
2014 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2015 Cstar = P128.w[1];
2016 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
2017 fstar.w[0] = P128.w[0];
2018 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2019 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2020 // C* = floor(C*) (logical right shift; C has p decimal digits,
2021 // correct by Property 1)
2022 // n = C* * 10^(e+x)
2023
2024 // shift right C* by Ex-64 = shiftright128[ind]
2025 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2026 Cstar = Cstar >> shift;
2027 // determine inexactness of the rounding of C*
2028 // if (0 < f* < 10^(-x)) then
2029 // the result is exact
2030 // else // if (f* > T*) then
2031 // the result is inexact
2032 if (ind - 1 <= 2) {
2033 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2034 // ten2mk128trunc[ind -1].w[1] is identical to
2035 // ten2mk128[ind -1].w[1]
2036 // set the inexact flag
2037 *pfpsf |= INEXACT_EXCEPTION;
2038 } // else the result is exact
2039 } else { // if 3 <= ind - 1 <= 14
2040 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2041 // ten2mk128trunc[ind -1].w[1] is identical to
2042 // ten2mk128[ind -1].w[1]
2043 // set the inexact flag
2044 *pfpsf |= INEXACT_EXCEPTION;
2045 } // else the result is exact
2046 }
2047
2048 if (x_sign)
2049 res = -Cstar;
2050 else
2051 res = Cstar;
2052 } else if (exp == 0) {
2053 // 1 <= q <= 10
2054 // res = +/-C (exact)
2055 if (x_sign)
2056 res = -C1;
2057 else
2058 res = C1;
2059 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2060 // res = +/-C * 10^exp (exact)
2061 if (x_sign)
2062 res = -C1 * ten2k64[exp];
2063 else
2064 res = C1 * ten2k64[exp];
2065 }
2066 }
2067 BID_RETURN (res);
2068 }
2069
2070 /*****************************************************************************
2071 * BID64_to_int32_rninta
2072 ****************************************************************************/
2073
2074 #if DECIMAL_CALL_BY_REFERENCE
2075 void
2076 bid64_to_int32_rninta (int *pres,
2077 UINT64 *
2078 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2079 _EXC_INFO_PARAM) {
2080 UINT64 x = *px;
2081 #else
2082 int
2083 bid64_to_int32_rninta (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2084 _EXC_INFO_PARAM) {
2085 #endif
2086 int res;
2087 UINT64 x_sign;
2088 UINT64 x_exp;
2089 int exp; // unbiased exponent
2090 // Note: C1 represents x_significand (UINT64)
2091 UINT64 tmp64;
2092 BID_UI64DOUBLE tmp1;
2093 unsigned int x_nr_bits;
2094 int q, ind, shift;
2095 UINT64 C1;
2096 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
2097 UINT128 P128;
2098
2099 // check for NaN or Infinity
2100 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
2101 // set invalid flag
2102 *pfpsf |= INVALID_EXCEPTION;
2103 // return Integer Indefinite
2104 res = 0x80000000;
2105 BID_RETURN (res);
2106 }
2107 // unpack x
2108 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2109 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2110 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
2111 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
2112 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
2113 if (C1 > 9999999999999999ull) { // non-canonical
2114 x_exp = 0;
2115 C1 = 0;
2116 }
2117 } else {
2118 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
2119 C1 = x & MASK_BINARY_SIG1;
2120 }
2121
2122 // check for zeros (possibly from non-canonical values)
2123 if (C1 == 0x0ull) {
2124 // x is 0
2125 res = 0x00000000;
2126 BID_RETURN (res);
2127 }
2128 // x is not special and is not zero
2129
2130 // q = nr. of decimal digits in x (1 <= q <= 54)
2131 // determine first the nr. of bits in x
2132 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
2133 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2134 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
2135 tmp1.d = (double) (C1 >> 32); // exact conversion
2136 x_nr_bits =
2137 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2138 } else { // x < 2^32
2139 tmp1.d = (double) C1; // exact conversion
2140 x_nr_bits =
2141 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2142 }
2143 } else { // if x < 2^53
2144 tmp1.d = (double) C1; // exact conversion
2145 x_nr_bits =
2146 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2147 }
2148 q = nr_digits[x_nr_bits - 1].digits;
2149 if (q == 0) {
2150 q = nr_digits[x_nr_bits - 1].digits1;
2151 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
2152 q++;
2153 }
2154 exp = x_exp - 398; // unbiased exponent
2155
2156 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2157 // set invalid flag
2158 *pfpsf |= INVALID_EXCEPTION;
2159 // return Integer Indefinite
2160 res = 0x80000000;
2161 BID_RETURN (res);
2162 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2163 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2164 // so x rounded to an integer may or may not fit in a signed 32-bit int
2165 // the cases that do not fit are identified here; the ones that fit
2166 // fall through and will be handled with other cases further,
2167 // under '1 <= q + exp <= 10'
2168 if (x_sign) { // if n < 0 and q + exp = 10
2169 // if n <= -2^31 - 1/2 then n is too large
2170 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
2171 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=16
2172 // <=> C * 10^(11-q) >= 0x500000005, 1<=q<=16
2173 if (q <= 11) {
2174 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2175 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2176 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2177 if (tmp64 >= 0x500000005ull) {
2178 // set invalid flag
2179 *pfpsf |= INVALID_EXCEPTION;
2180 // return Integer Indefinite
2181 res = 0x80000000;
2182 BID_RETURN (res);
2183 }
2184 // else cases that can be rounded to a 32-bit int fall through
2185 // to '1 <= q + exp <= 10'
2186 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2187 // C * 10^(11-q) >= 0x500000005 <=>
2188 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
2189 // (scale 2^31+1/2 up)
2190 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
2191 tmp64 = 0x500000005ull * ten2k64[q - 11];
2192 if (C1 >= tmp64) {
2193 // set invalid flag
2194 *pfpsf |= INVALID_EXCEPTION;
2195 // return Integer Indefinite
2196 res = 0x80000000;
2197 BID_RETURN (res);
2198 }
2199 // else cases that can be rounded to a 32-bit int fall through
2200 // to '1 <= q + exp <= 10'
2201 }
2202 } else { // if n > 0 and q + exp = 10
2203 // if n >= 2^31 - 1/2 then n is too large
2204 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
2205 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
2206 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
2207 if (q <= 11) {
2208 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2209 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2210 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2211 if (tmp64 >= 0x4fffffffbull) {
2212 // set invalid flag
2213 *pfpsf |= INVALID_EXCEPTION;
2214 // return Integer Indefinite
2215 res = 0x80000000;
2216 BID_RETURN (res);
2217 }
2218 // else cases that can be rounded to a 32-bit int fall through
2219 // to '1 <= q + exp <= 10'
2220 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2221 // C * 10^(11-q) >= 0x4fffffffb <=>
2222 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
2223 // (scale 2^31-1/2 up)
2224 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
2225 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
2226 if (C1 >= tmp64) {
2227 // set invalid flag
2228 *pfpsf |= INVALID_EXCEPTION;
2229 // return Integer Indefinite
2230 res = 0x80000000;
2231 BID_RETURN (res);
2232 }
2233 // else cases that can be rounded to a 32-bit int fall through
2234 // to '1 <= q + exp <= 10'
2235 }
2236 }
2237 }
2238 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
2239 // Note: some of the cases tested for above fall through to this point
2240 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2241 // return 0
2242 res = 0x00000000;
2243 BID_RETURN (res);
2244 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2245 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2246 // res = 0
2247 // else
2248 // res = +/-1
2249 ind = q - 1;
2250 if (C1 < midpoint64[ind]) {
2251 res = 0x00000000; // return 0
2252 } else if (x_sign) { // n < 0
2253 res = 0xffffffff; // return -1
2254 } else { // n > 0
2255 res = 0x00000001; // return +1
2256 }
2257 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2258 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
2259 // to nearest away to a 32-bit signed integer
2260 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2261 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2262 // chop off ind digits from the lower part of C1
2263 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2264 C1 = C1 + midpoint64[ind - 1];
2265 // calculate C* and f*
2266 // C* is actually floor(C*) in this case
2267 // C* and f* need shifting and masking, as shown by
2268 // shiftright128[] and maskhigh128[]
2269 // 1 <= x <= 15
2270 // kx = 10^(-x) = ten2mk64[ind - 1]
2271 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2272 // the approximation of 10^(-x) was rounded up to 54 bits
2273 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2274 Cstar = P128.w[1];
2275 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2276 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2277 // C* = floor(C*)-1 (logical right shift; C* has p decimal digits,
2278 // correct by Pr. 1)
2279 // n = C* * 10^(e+x)
2280
2281 // shift right C* by Ex-64 = shiftright128[ind]
2282 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2283 Cstar = Cstar >> shift;
2284
2285 // if the result was a midpoint it was rounded away from zero
2286 if (x_sign)
2287 res = -Cstar;
2288 else
2289 res = Cstar;
2290 } else if (exp == 0) {
2291 // 1 <= q <= 10
2292 // res = +/-C (exact)
2293 if (x_sign)
2294 res = -C1;
2295 else
2296 res = C1;
2297 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2298 // res = +/-C * 10^exp (exact)
2299 if (x_sign)
2300 res = -C1 * ten2k64[exp];
2301 else
2302 res = C1 * ten2k64[exp];
2303 }
2304 }
2305 BID_RETURN (res);
2306 }
2307
2308 /*****************************************************************************
2309 * BID64_to_int32_xrninta
2310 ****************************************************************************/
2311
2312 #if DECIMAL_CALL_BY_REFERENCE
2313 void
2314 bid64_to_int32_xrninta (int *pres,
2315 UINT64 *
2316 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2317 _EXC_INFO_PARAM) {
2318 UINT64 x = *px;
2319 #else
2320 int
2321 bid64_to_int32_xrninta (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2322 _EXC_INFO_PARAM) {
2323 #endif
2324 int res;
2325 UINT64 x_sign;
2326 UINT64 x_exp;
2327 int exp; // unbiased exponent
2328 // Note: C1 represents x_significand (UINT64)
2329 UINT64 tmp64;
2330 BID_UI64DOUBLE tmp1;
2331 unsigned int x_nr_bits;
2332 int q, ind, shift;
2333 UINT64 C1;
2334 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
2335 UINT128 fstar;
2336 UINT128 P128;
2337
2338 // check for NaN or Infinity
2339 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
2340 // set invalid flag
2341 *pfpsf |= INVALID_EXCEPTION;
2342 // return Integer Indefinite
2343 res = 0x80000000;
2344 BID_RETURN (res);
2345 }
2346 // unpack x
2347 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2348 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2349 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
2350 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
2351 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
2352 if (C1 > 9999999999999999ull) { // non-canonical
2353 x_exp = 0;
2354 C1 = 0;
2355 }
2356 } else {
2357 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
2358 C1 = x & MASK_BINARY_SIG1;
2359 }
2360
2361 // check for zeros (possibly from non-canonical values)
2362 if (C1 == 0x0ull) {
2363 // x is 0
2364 res = 0x00000000;
2365 BID_RETURN (res);
2366 }
2367 // x is not special and is not zero
2368
2369 // q = nr. of decimal digits in x (1 <= q <= 54)
2370 // determine first the nr. of bits in x
2371 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
2372 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2373 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
2374 tmp1.d = (double) (C1 >> 32); // exact conversion
2375 x_nr_bits =
2376 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2377 } else { // x < 2^32
2378 tmp1.d = (double) C1; // exact conversion
2379 x_nr_bits =
2380 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2381 }
2382 } else { // if x < 2^53
2383 tmp1.d = (double) C1; // exact conversion
2384 x_nr_bits =
2385 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2386 }
2387 q = nr_digits[x_nr_bits - 1].digits;
2388 if (q == 0) {
2389 q = nr_digits[x_nr_bits - 1].digits1;
2390 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
2391 q++;
2392 }
2393 exp = x_exp - 398; // unbiased exponent
2394
2395 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2396 // set invalid flag
2397 *pfpsf |= INVALID_EXCEPTION;
2398 // return Integer Indefinite
2399 res = 0x80000000;
2400 BID_RETURN (res);
2401 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2402 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2403 // so x rounded to an integer may or may not fit in a signed 32-bit int
2404 // the cases that do not fit are identified here; the ones that fit
2405 // fall through and will be handled with other cases further,
2406 // under '1 <= q + exp <= 10'
2407 if (x_sign) { // if n < 0 and q + exp = 10
2408 // if n <= -2^31 - 1/2 then n is too large
2409 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
2410 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=16
2411 // <=> C * 10^(11-q) >= 0x500000005, 1<=q<=16
2412 if (q <= 11) {
2413 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2414 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2415 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2416 if (tmp64 >= 0x500000005ull) {
2417 // set invalid flag
2418 *pfpsf |= INVALID_EXCEPTION;
2419 // return Integer Indefinite
2420 res = 0x80000000;
2421 BID_RETURN (res);
2422 }
2423 // else cases that can be rounded to a 32-bit int fall through
2424 // to '1 <= q + exp <= 10'
2425 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2426 // C * 10^(11-q) >= 0x500000005 <=>
2427 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
2428 // (scale 2^31+1/2 up)
2429 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
2430 tmp64 = 0x500000005ull * ten2k64[q - 11];
2431 if (C1 >= tmp64) {
2432 // set invalid flag
2433 *pfpsf |= INVALID_EXCEPTION;
2434 // return Integer Indefinite
2435 res = 0x80000000;
2436 BID_RETURN (res);
2437 }
2438 // else cases that can be rounded to a 32-bit int fall through
2439 // to '1 <= q + exp <= 10'
2440 }
2441 } else { // if n > 0 and q + exp = 10
2442 // if n >= 2^31 - 1/2 then n is too large
2443 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
2444 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
2445 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
2446 if (q <= 11) {
2447 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2448 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2449 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2450 if (tmp64 >= 0x4fffffffbull) {
2451 // set invalid flag
2452 *pfpsf |= INVALID_EXCEPTION;
2453 // return Integer Indefinite
2454 res = 0x80000000;
2455 BID_RETURN (res);
2456 }
2457 // else cases that can be rounded to a 32-bit int fall through
2458 // to '1 <= q + exp <= 10'
2459 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2460 // C * 10^(11-q) >= 0x4fffffffb <=>
2461 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
2462 // (scale 2^31-1/2 up)
2463 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
2464 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
2465 if (C1 >= tmp64) {
2466 // set invalid flag
2467 *pfpsf |= INVALID_EXCEPTION;
2468 // return Integer Indefinite
2469 res = 0x80000000;
2470 BID_RETURN (res);
2471 }
2472 // else cases that can be rounded to a 32-bit int fall through
2473 // to '1 <= q + exp <= 10'
2474 }
2475 }
2476 }
2477 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
2478 // Note: some of the cases tested for above fall through to this point
2479 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2480 // set inexact flag
2481 *pfpsf |= INEXACT_EXCEPTION;
2482 // return 0
2483 res = 0x00000000;
2484 BID_RETURN (res);
2485 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2486 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2487 // res = 0
2488 // else
2489 // res = +/-1
2490 ind = q - 1;
2491 if (C1 < midpoint64[ind]) {
2492 res = 0x00000000; // return 0
2493 } else if (x_sign) { // n < 0
2494 res = 0xffffffff; // return -1
2495 } else { // n > 0
2496 res = 0x00000001; // return +1
2497 }
2498 // set inexact flag
2499 *pfpsf |= INEXACT_EXCEPTION;
2500 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2501 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
2502 // to nearest away to a 32-bit signed integer
2503 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2504 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2505 // chop off ind digits from the lower part of C1
2506 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2507 C1 = C1 + midpoint64[ind - 1];
2508 // calculate C* and f*
2509 // C* is actually floor(C*) in this case
2510 // C* and f* need shifting and masking, as shown by
2511 // shiftright128[] and maskhigh128[]
2512 // 1 <= x <= 15
2513 // kx = 10^(-x) = ten2mk64[ind - 1]
2514 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2515 // the approximation of 10^(-x) was rounded up to 54 bits
2516 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2517 Cstar = P128.w[1];
2518 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
2519 fstar.w[0] = P128.w[0];
2520 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2521 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2522 // C* = floor(C*)-1 (logical right shift; C* has p decimal digits,
2523 // correct by Pr. 1)
2524 // n = C* * 10^(e+x)
2525
2526 // shift right C* by Ex-64 = shiftright128[ind]
2527 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2528 Cstar = Cstar >> shift;
2529 // determine inexactness of the rounding of C*
2530 // if (0 < f* - 1/2 < 10^(-x)) then
2531 // the result is exact
2532 // else // if (f* - 1/2 > T*) then
2533 // the result is inexact
2534 if (ind - 1 <= 2) {
2535 if (fstar.w[0] > 0x8000000000000000ull) {
2536 // f* > 1/2 and the result may be exact
2537 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
2538 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
2539 // ten2mk128trunc[ind -1].w[1] is identical to
2540 // ten2mk128[ind -1].w[1]
2541 // set the inexact flag
2542 *pfpsf |= INEXACT_EXCEPTION;
2543 } // else the result is exact
2544 } else { // the result is inexact; f2* <= 1/2
2545 // set the inexact flag
2546 *pfpsf |= INEXACT_EXCEPTION;
2547 }
2548 } else { // if 3 <= ind - 1 <= 14
2549 if (fstar.w[1] > onehalf128[ind - 1] ||
2550 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
2551 // f2* > 1/2 and the result may be exact
2552 // Calculate f2* - 1/2
2553 tmp64 = fstar.w[1] - onehalf128[ind - 1];
2554 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2555 // ten2mk128trunc[ind -1].w[1] is identical to
2556 // ten2mk128[ind -1].w[1]
2557 // set the inexact flag
2558 *pfpsf |= INEXACT_EXCEPTION;
2559 } // else the result is exact
2560 } else { // the result is inexact; f2* <= 1/2
2561 // set the inexact flag
2562 *pfpsf |= INEXACT_EXCEPTION;
2563 }
2564 }
2565
2566 // if the result was a midpoint it was rounded away from zero
2567 if (x_sign)
2568 res = -Cstar;
2569 else
2570 res = Cstar;
2571 } else if (exp == 0) {
2572 // 1 <= q <= 10
2573 // res = +/-C (exact)
2574 if (x_sign)
2575 res = -C1;
2576 else
2577 res = C1;
2578 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2579 // res = +/-C * 10^exp (exact)
2580 if (x_sign)
2581 res = -C1 * ten2k64[exp];
2582 else
2583 res = C1 * ten2k64[exp];
2584 }
2585 }
2586 BID_RETURN (res);
2587 }