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