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