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