1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc. 2 3 This file is part of GCC. 4 5 GCC is free software; you can redistribute it and/or modify it under 6 the terms of the GNU General Public License as published by the Free 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11 WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13 for more details. 14 15 Under Section 7 of GPL version 3, you are granted additional 16 permissions described in the GCC Runtime Library Exception, version 17 3.1, as published by the Free Software Foundation. 18 19 You should have received a copy of the GNU General Public License and 20 a copy of the GCC Runtime Library Exception along with this program; 21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22 <http://www.gnu.org/licenses/>. */ 23 24 #include "bid_internal.h" 25 26 /***************************************************************************** 27 * BID64_to_int32_rnint 28 ****************************************************************************/ 29 30 #if DECIMAL_CALL_BY_REFERENCE 31 void 32 bid64_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 38 int 39 bid64_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 290 void 291 bid64_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 297 int 298 bid64_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 589 void 590 bid64_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 596 int 597 bid64_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 840 void 841 bid64_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 847 int 848 bid64_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 1097 void 1098 bid64_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 1104 int 1105 bid64_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 1348 void 1349 bid64_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 1355 int 1356 bid64_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 1605 void 1606 bid64_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 1612 int 1613 bid64_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 1827 void 1828 bid64_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 1834 int 1835 bid64_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 2075 void 2076 bid64_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 2082 int 2083 bid64_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 2313 void 2314 bid64_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 2320 int 2321 bid64_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