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