1/* Copyright (C) 2007-2022 Free Software Foundation, Inc. 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 7Software Foundation; either version 3, or (at your option) any later 8version. 9 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 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/>. */ 23 24#include "bid_internal.h" 25 26/***************************************************************************** 27 * BID64_to_int64_rnint 28 ****************************************************************************/ 29 30#if DECIMAL_CALL_BY_REFERENCE 31void 32bid64_to_int64_rnint (SINT64 * pres, UINT64 * px 33 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 34 _EXC_INFO_PARAM) { 35 UINT64 x = *px; 36#else 37SINT64 38bid64_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 263void 264bid64_to_int64_xrnint (SINT64 * pres, UINT64 * px 265 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 266 _EXC_INFO_PARAM) { 267 UINT64 x = *px; 268#else 269SINT64 270bid64_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 536void 537bid64_to_int64_floor (SINT64 * pres, UINT64 * px 538 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 539 _EXC_INFO_PARAM) { 540 UINT64 x = *px; 541#else 542SINT64 543bid64_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 760void 761bid64_to_int64_xfloor (SINT64 * pres, UINT64 * px 762 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 763 _EXC_INFO_PARAM) { 764 UINT64 x = *px; 765#else 766SINT64 767bid64_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 990void 991bid64_to_int64_ceil (SINT64 * pres, UINT64 * px 992 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 993{ 994 UINT64 x = *px; 995#else 996SINT64 997bid64_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 1214void 1215bid64_to_int64_xceil (SINT64 * pres, UINT64 * px 1216 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1217 _EXC_INFO_PARAM) { 1218 UINT64 x = *px; 1219#else 1220SINT64 1221bid64_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 1444void 1445bid64_to_int64_int (SINT64 * pres, UINT64 * px 1446 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 1447 UINT64 x = *px; 1448#else 1449SINT64 1450bid64_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 1638void 1639bid64_to_int64_xint (SINT64 * pres, UINT64 * px 1640 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 1641{ 1642 UINT64 x = *px; 1643#else 1644SINT64 1645bid64_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 1858void 1859bid64_to_int64_rninta (SINT64 * pres, UINT64 * px 1860 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 1861 _EXC_INFO_PARAM) { 1862 UINT64 x = *px; 1863#else 1864SINT64 1865bid64_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 2075void 2076bid64_to_int64_xrninta (SINT64 * pres, UINT64 * px 2077 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 2078 _EXC_INFO_PARAM) { 2079 UINT64 x = *px; 2080#else 2081SINT64 2082bid64_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