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