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