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