1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc. 2 3 This file is part of GCC. 4 5 GCC is free software; you can redistribute it and/or modify it under 6 the terms of the GNU General Public License as published by the Free 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11 WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13 for more details. 14 15 Under Section 7 of GPL version 3, you are granted additional 16 permissions described in the GCC Runtime Library Exception, version 17 3.1, as published by the Free Software Foundation. 18 19 You should have received a copy of the GNU General Public License and 20 a copy of the GCC Runtime Library Exception along with this program; 21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22 <http://www.gnu.org/licenses/>. */ 23 24 #include "bid_internal.h" 25 26 /***************************************************************************** 27 * BID64 nextup 28 ****************************************************************************/ 29 30 #if DECIMAL_CALL_BY_REFERENCE 31 void 32 bid64_nextup (UINT64 * pres, 33 UINT64 * 34 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 35 UINT64 x = *px; 36 #else 37 UINT64 38 bid64_nextup (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 39 _EXC_INFO_PARAM) { 40 #endif 41 42 UINT64 res; 43 UINT64 x_sign; 44 UINT64 x_exp; 45 BID_UI64DOUBLE tmp1; 46 int x_nr_bits; 47 int q1, ind; 48 UINT64 C1; // C1 represents x_signif (UINT64) 49 50 // check for NaNs and infinities 51 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 52 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 53 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 54 else 55 x = x & 0xfe03ffffffffffffull; // clear G6-G12 56 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 57 // set invalid flag 58 *pfpsf |= INVALID_EXCEPTION; 59 // return quiet (SNaN) 60 res = x & 0xfdffffffffffffffull; 61 } else { // QNaN 62 res = x; 63 } 64 BID_RETURN (res); 65 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 66 if (!(x & 0x8000000000000000ull)) { // x is +inf 67 res = 0x7800000000000000ull; 68 } else { // x is -inf 69 res = 0xf7fb86f26fc0ffffull; // -MAXFP = -999...99 * 10^emax 70 } 71 BID_RETURN (res); 72 } 73 // unpack the argument 74 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 75 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 76 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 77 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 78 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 79 if (C1 > 9999999999999999ull) { // non-canonical 80 x_exp = 0; 81 C1 = 0; 82 } 83 } else { 84 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 85 C1 = x & MASK_BINARY_SIG1; 86 } 87 88 // check for zeros (possibly from non-canonical values) 89 if (C1 == 0x0ull) { 90 // x is 0 91 res = 0x0000000000000001ull; // MINFP = 1 * 10^emin 92 } else { // x is not special and is not zero 93 if (x == 0x77fb86f26fc0ffffull) { 94 // x = +MAXFP = 999...99 * 10^emax 95 res = 0x7800000000000000ull; // +inf 96 } else if (x == 0x8000000000000001ull) { 97 // x = -MINFP = 1...99 * 10^emin 98 res = 0x8000000000000000ull; // -0 99 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 100 // can add/subtract 1 ulp to the significand 101 102 // Note: we could check here if x >= 10^16 to speed up the case q1 =16 103 // q1 = nr. of decimal digits in x (1 <= q1 <= 54) 104 // determine first the nr. of bits in x 105 if (C1 >= MASK_BINARY_OR2) { // x >= 2^53 106 // split the 64-bit value in two 32-bit halves to avoid rounding errors 107 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 108 tmp1.d = (double) (C1 >> 32); // exact conversion 109 x_nr_bits = 110 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 111 } else { // x < 2^32 112 tmp1.d = (double) C1; // exact conversion 113 x_nr_bits = 114 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 115 } 116 } else { // if x < 2^53 117 tmp1.d = (double) C1; // exact conversion 118 x_nr_bits = 119 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 120 } 121 q1 = nr_digits[x_nr_bits - 1].digits; 122 if (q1 == 0) { 123 q1 = nr_digits[x_nr_bits - 1].digits1; 124 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 125 q1++; 126 } 127 // if q1 < P16 then pad the significand with zeros 128 if (q1 < P16) { 129 if (x_exp > (UINT64) (P16 - q1)) { 130 ind = P16 - q1; // 1 <= ind <= P16 - 1 131 // pad with P16 - q1 zeros, until exponent = emin 132 // C1 = C1 * 10^ind 133 C1 = C1 * ten2k64[ind]; 134 x_exp = x_exp - ind; 135 } else { // pad with zeros until the exponent reaches emin 136 ind = x_exp; 137 C1 = C1 * ten2k64[ind]; 138 x_exp = EXP_MIN; 139 } 140 } 141 if (!x_sign) { // x > 0 142 // add 1 ulp (add 1 to the significand) 143 C1++; 144 if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16 145 C1 = 0x00038d7ea4c68000ull; // C1 = 10^15 146 x_exp++; 147 } 148 // Ok, because MAXFP = 999...99 * 10^emax was caught already 149 } else { // x < 0 150 // subtract 1 ulp (subtract 1 from the significand) 151 C1--; 152 if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1 153 C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1 154 x_exp--; 155 } 156 } 157 // assemble the result 158 // if significand has 54 bits 159 if (C1 & MASK_BINARY_OR2) { 160 res = 161 x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 & 162 MASK_BINARY_SIG2); 163 } else { // significand fits in 53 bits 164 res = x_sign | (x_exp << 53) | C1; 165 } 166 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 167 } // end x is not special and is not zero 168 BID_RETURN (res); 169 } 170 171 /***************************************************************************** 172 * BID64 nextdown 173 ****************************************************************************/ 174 175 #if DECIMAL_CALL_BY_REFERENCE 176 void 177 bid64_nextdown (UINT64 * pres, 178 UINT64 * 179 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 180 UINT64 x = *px; 181 #else 182 UINT64 183 bid64_nextdown (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 184 _EXC_INFO_PARAM) { 185 #endif 186 187 UINT64 res; 188 UINT64 x_sign; 189 UINT64 x_exp; 190 BID_UI64DOUBLE tmp1; 191 int x_nr_bits; 192 int q1, ind; 193 UINT64 C1; // C1 represents x_signif (UINT64) 194 195 // check for NaNs and infinities 196 if ((x & MASK_NAN) == MASK_NAN) { // check for NaN 197 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 198 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 199 else 200 x = x & 0xfe03ffffffffffffull; // clear G6-G12 201 if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN 202 // set invalid flag 203 *pfpsf |= INVALID_EXCEPTION; 204 // return quiet (SNaN) 205 res = x & 0xfdffffffffffffffull; 206 } else { // QNaN 207 res = x; 208 } 209 BID_RETURN (res); 210 } else if ((x & MASK_INF) == MASK_INF) { // check for Infinity 211 if (x & 0x8000000000000000ull) { // x is -inf 212 res = 0xf800000000000000ull; 213 } else { // x is +inf 214 res = 0x77fb86f26fc0ffffull; // +MAXFP = +999...99 * 10^emax 215 } 216 BID_RETURN (res); 217 } 218 // unpack the argument 219 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 220 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] => 221 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 222 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased 223 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; 224 if (C1 > 9999999999999999ull) { // non-canonical 225 x_exp = 0; 226 C1 = 0; 227 } 228 } else { 229 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased 230 C1 = x & MASK_BINARY_SIG1; 231 } 232 233 // check for zeros (possibly from non-canonical values) 234 if (C1 == 0x0ull) { 235 // x is 0 236 res = 0x8000000000000001ull; // -MINFP = -1 * 10^emin 237 } else { // x is not special and is not zero 238 if (x == 0xf7fb86f26fc0ffffull) { 239 // x = -MAXFP = -999...99 * 10^emax 240 res = 0xf800000000000000ull; // -inf 241 } else if (x == 0x0000000000000001ull) { 242 // x = +MINFP = 1...99 * 10^emin 243 res = 0x0000000000000000ull; // -0 244 } else { // -MAXFP + 1ulp <= x <= -MINFP OR MINFP + 1 ulp <= x <= MAXFP 245 // can add/subtract 1 ulp to the significand 246 247 // Note: we could check here if x >= 10^16 to speed up the case q1 =16 248 // q1 = nr. of decimal digits in x (1 <= q1 <= 16) 249 // determine first the nr. of bits in x 250 if (C1 >= 0x0020000000000000ull) { // x >= 2^53 251 // split the 64-bit value in two 32-bit halves to avoid 252 // rounding errors 253 if (C1 >= 0x0000000100000000ull) { // x >= 2^32 254 tmp1.d = (double) (C1 >> 32); // exact conversion 255 x_nr_bits = 256 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 257 } else { // x < 2^32 258 tmp1.d = (double) C1; // exact conversion 259 x_nr_bits = 260 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 261 } 262 } else { // if x < 2^53 263 tmp1.d = (double) C1; // exact conversion 264 x_nr_bits = 265 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 266 } 267 q1 = nr_digits[x_nr_bits - 1].digits; 268 if (q1 == 0) { 269 q1 = nr_digits[x_nr_bits - 1].digits1; 270 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) 271 q1++; 272 } 273 // if q1 < P16 then pad the significand with zeros 274 if (q1 < P16) { 275 if (x_exp > (UINT64) (P16 - q1)) { 276 ind = P16 - q1; // 1 <= ind <= P16 - 1 277 // pad with P16 - q1 zeros, until exponent = emin 278 // C1 = C1 * 10^ind 279 C1 = C1 * ten2k64[ind]; 280 x_exp = x_exp - ind; 281 } else { // pad with zeros until the exponent reaches emin 282 ind = x_exp; 283 C1 = C1 * ten2k64[ind]; 284 x_exp = EXP_MIN; 285 } 286 } 287 if (x_sign) { // x < 0 288 // add 1 ulp (add 1 to the significand) 289 C1++; 290 if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16 291 C1 = 0x00038d7ea4c68000ull; // C1 = 10^15 292 x_exp++; 293 // Ok, because -MAXFP = -999...99 * 10^emax was caught already 294 } 295 } else { // x > 0 296 // subtract 1 ulp (subtract 1 from the significand) 297 C1--; 298 if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1 299 C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1 300 x_exp--; 301 } 302 } 303 // assemble the result 304 // if significand has 54 bits 305 if (C1 & MASK_BINARY_OR2) { 306 res = 307 x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 & 308 MASK_BINARY_SIG2); 309 } else { // significand fits in 53 bits 310 res = x_sign | (x_exp << 53) | C1; 311 } 312 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 313 } // end x is not special and is not zero 314 BID_RETURN (res); 315 } 316 317 /***************************************************************************** 318 * BID64 nextafter 319 ****************************************************************************/ 320 321 #if DECIMAL_CALL_BY_REFERENCE 322 void 323 bid64_nextafter (UINT64 * pres, UINT64 * px, 324 UINT64 * 325 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 326 UINT64 x = *px; 327 UINT64 y = *py; 328 #else 329 UINT64 330 bid64_nextafter (UINT64 x, 331 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 332 _EXC_INFO_PARAM) { 333 #endif 334 335 UINT64 res; 336 UINT64 tmp1, tmp2; 337 FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions 338 int res1, res2; 339 340 // check for NaNs or infinities 341 if (((x & MASK_SPECIAL) == MASK_SPECIAL) || 342 ((y & MASK_SPECIAL) == MASK_SPECIAL)) { 343 // x is NaN or infinity or y is NaN or infinity 344 345 if ((x & MASK_NAN) == MASK_NAN) { // x is NAN 346 if ((x & 0x0003ffffffffffffull) > 999999999999999ull) 347 x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 348 else 349 x = x & 0xfe03ffffffffffffull; // clear G6-G12 350 if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN 351 // set invalid flag 352 *pfpsf |= INVALID_EXCEPTION; 353 // return quiet (x) 354 res = x & 0xfdffffffffffffffull; 355 } else { // x is QNaN 356 if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN 357 // set invalid flag 358 *pfpsf |= INVALID_EXCEPTION; 359 } 360 // return x 361 res = x; 362 } 363 BID_RETURN (res); 364 } else if ((y & MASK_NAN) == MASK_NAN) { // y is NAN 365 if ((y & 0x0003ffffffffffffull) > 999999999999999ull) 366 y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits 367 else 368 y = y & 0xfe03ffffffffffffull; // clear G6-G12 369 if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN 370 // set invalid flag 371 *pfpsf |= INVALID_EXCEPTION; 372 // return quiet (y) 373 res = y & 0xfdffffffffffffffull; 374 } else { // y is QNaN 375 // return y 376 res = y; 377 } 378 BID_RETURN (res); 379 } else { // at least one is infinity 380 if ((x & MASK_ANY_INF) == MASK_INF) { // x = inf 381 x = x & (MASK_SIGN | MASK_INF); 382 } 383 if ((y & MASK_ANY_INF) == MASK_INF) { // y = inf 384 y = y & (MASK_SIGN | MASK_INF); 385 } 386 } 387 } 388 // neither x nor y is NaN 389 390 // if not infinity, check for non-canonical values x (treated as zero) 391 if ((x & MASK_ANY_INF) != MASK_INF) { // x != inf 392 // unpack x 393 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) { 394 // if the steering bits are 11 (condition will be 0), then 395 // the exponent is G[0:w+1] 396 if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) > 397 9999999999999999ull) { 398 // non-canonical 399 x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2); 400 } 401 } else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) x is unch. 402 ; // canonical 403 } 404 } 405 // no need to check for non-canonical y 406 407 // neither x nor y is NaN 408 tmp_fpsf = *pfpsf; // save fpsf 409 #if DECIMAL_CALL_BY_REFERENCE 410 bid64_quiet_equal (&res1, px, 411 py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 412 bid64_quiet_greater (&res2, px, 413 py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 414 #else 415 res1 = 416 bid64_quiet_equal (x, 417 y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 418 res2 = 419 bid64_quiet_greater (x, 420 y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 421 #endif 422 *pfpsf = tmp_fpsf; // restore fpsf 423 if (res1) { // x = y 424 // return x with the sign of y 425 res = (y & 0x8000000000000000ull) | (x & 0x7fffffffffffffffull); 426 } else if (res2) { // x > y 427 #if DECIMAL_CALL_BY_REFERENCE 428 bid64_nextdown (&res, 429 px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 430 #else 431 res = 432 bid64_nextdown (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 433 #endif 434 } else { // x < y 435 #if DECIMAL_CALL_BY_REFERENCE 436 bid64_nextup (&res, px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 437 #else 438 res = bid64_nextup (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 439 #endif 440 } 441 // if the operand x is finite but the result is infinite, signal 442 // overflow and inexact 443 if (((x & MASK_INF) != MASK_INF) && ((res & MASK_INF) == MASK_INF)) { 444 // set the inexact flag 445 *pfpsf |= INEXACT_EXCEPTION; 446 // set the overflow flag 447 *pfpsf |= OVERFLOW_EXCEPTION; 448 } 449 // if the result is in (-10^emin, 10^emin), and is different from the 450 // operand x, signal underflow and inexact 451 tmp1 = 0x00038d7ea4c68000ull; // +100...0[16] * 10^emin 452 tmp2 = res & 0x7fffffffffffffffull; 453 tmp_fpsf = *pfpsf; // save fpsf 454 #if DECIMAL_CALL_BY_REFERENCE 455 bid64_quiet_greater (&res1, &tmp1, 456 &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG 457 _EXC_INFO_ARG); 458 bid64_quiet_not_equal (&res2, &x, 459 &res _EXC_FLAGS_ARG _EXC_MASKS_ARG 460 _EXC_INFO_ARG); 461 #else 462 res1 = 463 bid64_quiet_greater (tmp1, 464 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG 465 _EXC_INFO_ARG); 466 res2 = 467 bid64_quiet_not_equal (x, 468 res _EXC_FLAGS_ARG _EXC_MASKS_ARG 469 _EXC_INFO_ARG); 470 #endif 471 *pfpsf = tmp_fpsf; // restore fpsf 472 if (res1 && res2) { 473 // if (bid64_quiet_greater (tmp1, tmp2, &tmp_fpsf) && 474 // bid64_quiet_not_equal (x, res, &tmp_fpsf)) { 475 // set the inexact flag 476 *pfpsf |= INEXACT_EXCEPTION; 477 // set the underflow flag 478 *pfpsf |= UNDERFLOW_EXCEPTION; 479 } 480 BID_RETURN (res); 481 } 482