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 #define BID_128RES 25 #include "bid_internal.h" 26 27 /***************************************************************************** 28 * BID128 nextup 29 ****************************************************************************/ 30 31 #if DECIMAL_CALL_BY_REFERENCE 32 void 33 bid128_nextup (UINT128 * pres, 34 UINT128 * 35 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 36 UINT128 x = *px; 37 #else 38 UINT128 39 bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 40 _EXC_INFO_PARAM) { 41 #endif 42 43 UINT128 res; 44 UINT64 x_sign; 45 UINT64 x_exp; 46 int exp; 47 BID_UI64DOUBLE tmp1; 48 int x_nr_bits; 49 int q1, ind; 50 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64) 51 52 BID_SWAP128 (x); 53 // unpack the argument 54 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 55 C1.w[1] = x.w[1] & MASK_COEFF; 56 C1.w[0] = x.w[0]; 57 58 // check for NaN or Infinity 59 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 60 // x is special 61 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 62 // if x = NaN, then res = Q (x) 63 // check first for non-canonical NaN payload 64 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 65 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) 66 && (x.w[0] > 0x38c15b09ffffffffull))) { 67 x.w[1] = x.w[1] & 0xffffc00000000000ull; 68 x.w[0] = 0x0ull; 69 } 70 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 71 // set invalid flag 72 *pfpsf |= INVALID_EXCEPTION; 73 // return quiet (x) 74 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 75 res.w[0] = x.w[0]; 76 } else { // x is QNaN 77 // return x 78 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 79 res.w[0] = x.w[0]; 80 } 81 } else { // x is not NaN, so it must be infinity 82 if (!x_sign) { // x is +inf 83 res.w[1] = 0x7800000000000000ull; // +inf 84 res.w[0] = 0x0000000000000000ull; 85 } else { // x is -inf 86 res.w[1] = 0xdfffed09bead87c0ull; // -MAXFP = -999...99 * 10^emax 87 res.w[0] = 0x378d8e63ffffffffull; 88 } 89 } 90 BID_RETURN (res); 91 } 92 // check for non-canonical values (treated as zero) 93 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 94 // non-canonical 95 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 96 C1.w[1] = 0; // significand high 97 C1.w[0] = 0; // significand low 98 } else { // G0_G1 != 11 99 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 100 if (C1.w[1] > 0x0001ed09bead87c0ull || 101 (C1.w[1] == 0x0001ed09bead87c0ull 102 && C1.w[0] > 0x378d8e63ffffffffull)) { 103 // x is non-canonical if coefficient is larger than 10^34 -1 104 C1.w[1] = 0; 105 C1.w[0] = 0; 106 } else { // canonical 107 ; 108 } 109 } 110 111 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 112 // x is +/-0 113 res.w[1] = 0x0000000000000000ull; // +1 * 10^emin 114 res.w[0] = 0x0000000000000001ull; 115 } else { // x is not special and is not zero 116 if (x.w[1] == 0x5fffed09bead87c0ull 117 && x.w[0] == 0x378d8e63ffffffffull) { 118 // x = +MAXFP = 999...99 * 10^emax 119 res.w[1] = 0x7800000000000000ull; // +inf 120 res.w[0] = 0x0000000000000000ull; 121 } else if (x.w[1] == 0x8000000000000000ull 122 && x.w[0] == 0x0000000000000001ull) { 123 // x = -MINFP = 1...99 * 10^emin 124 res.w[1] = 0x8000000000000000ull; // -0 125 res.w[0] = 0x0000000000000000ull; 126 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 127 // can add/subtract 1 ulp to the significand 128 129 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34 130 // q1 = nr. of decimal digits in x 131 // determine first the nr. of bits in x 132 if (C1.w[1] == 0) { 133 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 134 // split the 64-bit value in two 32-bit halves to avoid rnd errors 135 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 136 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 137 x_nr_bits = 138 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 139 0x3ff); 140 } else { // x < 2^32 141 tmp1.d = (double) (C1.w[0]); // exact conversion 142 x_nr_bits = 143 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 144 0x3ff); 145 } 146 } else { // if x < 2^53 147 tmp1.d = (double) C1.w[0]; // exact conversion 148 x_nr_bits = 149 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 150 } 151 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 152 tmp1.d = (double) C1.w[1]; // exact conversion 153 x_nr_bits = 154 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 155 } 156 q1 = nr_digits[x_nr_bits - 1].digits; 157 if (q1 == 0) { 158 q1 = nr_digits[x_nr_bits - 1].digits1; 159 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi 160 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi 161 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 162 q1++; 163 } 164 // if q1 < P34 then pad the significand with zeros 165 if (q1 < P34) { 166 exp = (x_exp >> 49) - 6176; 167 if (exp + 6176 > P34 - q1) { 168 ind = P34 - q1; // 1 <= ind <= P34 - 1 169 // pad with P34 - q1 zeros, until exponent = emin 170 // C1 = C1 * 10^ind 171 if (q1 <= 19) { // 64-bit C1 172 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 173 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 174 } else { // 128-bit 10^ind and 64-bit C1 175 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 176 } 177 } else { // C1 is (most likely) 128-bit 178 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely) 179 __mul_128x64_to_128 (C1, ten2k64[ind], C1); 180 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19) 181 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 182 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit) 183 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 184 } 185 } 186 x_exp = x_exp - ((UINT64) ind << 49); 187 } else { // pad with zeros until the exponent reaches emin 188 ind = exp + 6176; 189 // C1 = C1 * 10^ind 190 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33 191 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind 192 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 193 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind 194 __mul_128x64_to_128 (C1, ten2k64[ind], C1); 195 } 196 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 => 197 // 64-bit C1, 128-bit 10^ind 198 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 199 } 200 x_exp = EXP_MIN; 201 } 202 } 203 if (!x_sign) { // x > 0 204 // add 1 ulp (add 1 to the significand) 205 C1.w[0]++; 206 if (C1.w[0] == 0) 207 C1.w[1]++; 208 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34 209 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33 210 C1.w[0] = 0x38c15b0a00000000ull; 211 x_exp = x_exp + EXP_P1; 212 } 213 } else { // x < 0 214 // subtract 1 ulp (subtract 1 from the significand) 215 C1.w[0]--; 216 if (C1.w[0] == 0xffffffffffffffffull) 217 C1.w[1]--; 218 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1 219 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1 220 C1.w[0] = 0x378d8e63ffffffffull; 221 x_exp = x_exp - EXP_P1; 222 } 223 } 224 // assemble the result 225 res.w[1] = x_sign | x_exp | C1.w[1]; 226 res.w[0] = C1.w[0]; 227 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 228 } // end x is not special and is not zero 229 BID_RETURN (res); 230 } 231 232 /***************************************************************************** 233 * BID128 nextdown 234 ****************************************************************************/ 235 236 #if DECIMAL_CALL_BY_REFERENCE 237 void 238 bid128_nextdown (UINT128 * pres, 239 UINT128 * 240 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 241 UINT128 x = *px; 242 #else 243 UINT128 244 bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 245 _EXC_INFO_PARAM) { 246 #endif 247 248 UINT128 res; 249 UINT64 x_sign; 250 UINT64 x_exp; 251 int exp; 252 BID_UI64DOUBLE tmp1; 253 int x_nr_bits; 254 int q1, ind; 255 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64) 256 257 BID_SWAP128 (x); 258 // unpack the argument 259 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative 260 C1.w[1] = x.w[1] & MASK_COEFF; 261 C1.w[0] = x.w[0]; 262 263 // check for NaN or Infinity 264 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) { 265 // x is special 266 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 267 // if x = NaN, then res = Q (x) 268 // check first for non-canonical NaN payload 269 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 270 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) 271 && (x.w[0] > 0x38c15b09ffffffffull))) { 272 x.w[1] = x.w[1] & 0xffffc00000000000ull; 273 x.w[0] = 0x0ull; 274 } 275 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 276 // set invalid flag 277 *pfpsf |= INVALID_EXCEPTION; 278 // return quiet (x) 279 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 280 res.w[0] = x.w[0]; 281 } else { // x is QNaN 282 // return x 283 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 284 res.w[0] = x.w[0]; 285 } 286 } else { // x is not NaN, so it must be infinity 287 if (!x_sign) { // x is +inf 288 res.w[1] = 0x5fffed09bead87c0ull; // +MAXFP = +999...99 * 10^emax 289 res.w[0] = 0x378d8e63ffffffffull; 290 } else { // x is -inf 291 res.w[1] = 0xf800000000000000ull; // -inf 292 res.w[0] = 0x0000000000000000ull; 293 } 294 } 295 BID_RETURN (res); 296 } 297 // check for non-canonical values (treated as zero) 298 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 299 // non-canonical 300 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 301 C1.w[1] = 0; // significand high 302 C1.w[0] = 0; // significand low 303 } else { // G0_G1 != 11 304 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 305 if (C1.w[1] > 0x0001ed09bead87c0ull || 306 (C1.w[1] == 0x0001ed09bead87c0ull 307 && C1.w[0] > 0x378d8e63ffffffffull)) { 308 // x is non-canonical if coefficient is larger than 10^34 -1 309 C1.w[1] = 0; 310 C1.w[0] = 0; 311 } else { // canonical 312 ; 313 } 314 } 315 316 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) { 317 // x is +/-0 318 res.w[1] = 0x8000000000000000ull; // -1 * 10^emin 319 res.w[0] = 0x0000000000000001ull; 320 } else { // x is not special and is not zero 321 if (x.w[1] == 0xdfffed09bead87c0ull 322 && x.w[0] == 0x378d8e63ffffffffull) { 323 // x = -MAXFP = -999...99 * 10^emax 324 res.w[1] = 0xf800000000000000ull; // -inf 325 res.w[0] = 0x0000000000000000ull; 326 } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) { // +MINFP 327 res.w[1] = 0x0000000000000000ull; // +0 328 res.w[0] = 0x0000000000000000ull; 329 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 330 // can add/subtract 1 ulp to the significand 331 332 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34 333 // q1 = nr. of decimal digits in x 334 // determine first the nr. of bits in x 335 if (C1.w[1] == 0) { 336 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53 337 // split the 64-bit value in two 32-bit halves to avoid rnd errors 338 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32 339 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion 340 x_nr_bits = 341 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 342 0x3ff); 343 } else { // x < 2^32 344 tmp1.d = (double) (C1.w[0]); // exact conversion 345 x_nr_bits = 346 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 347 0x3ff); 348 } 349 } else { // if x < 2^53 350 tmp1.d = (double) C1.w[0]; // exact conversion 351 x_nr_bits = 352 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 353 } 354 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1]) 355 tmp1.d = (double) C1.w[1]; // exact conversion 356 x_nr_bits = 357 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); 358 } 359 q1 = nr_digits[x_nr_bits - 1].digits; 360 if (q1 == 0) { 361 q1 = nr_digits[x_nr_bits - 1].digits1; 362 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi 363 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi 364 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo)) 365 q1++; 366 } 367 // if q1 < P then pad the significand with zeros 368 if (q1 < P34) { 369 exp = (x_exp >> 49) - 6176; 370 if (exp + 6176 > P34 - q1) { 371 ind = P34 - q1; // 1 <= ind <= P34 - 1 372 // pad with P34 - q1 zeros, until exponent = emin 373 // C1 = C1 * 10^ind 374 if (q1 <= 19) { // 64-bit C1 375 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 376 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 377 } else { // 128-bit 10^ind and 64-bit C1 378 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 379 } 380 } else { // C1 is (most likely) 128-bit 381 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely) 382 __mul_128x64_to_128 (C1, ten2k64[ind], C1); 383 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19) 384 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 385 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit) 386 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 387 } 388 } 389 x_exp = x_exp - ((UINT64) ind << 49); 390 } else { // pad with zeros until the exponent reaches emin 391 ind = exp + 6176; 392 // C1 = C1 * 10^ind 393 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33 394 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind 395 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]); 396 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind 397 __mul_128x64_to_128 (C1, ten2k64[ind], C1); 398 } 399 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 => 400 // 64-bit C1, 128-bit 10^ind 401 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]); 402 } 403 x_exp = EXP_MIN; 404 } 405 } 406 if (x_sign) { // x < 0 407 // add 1 ulp (add 1 to the significand) 408 C1.w[0]++; 409 if (C1.w[0] == 0) 410 C1.w[1]++; 411 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34 412 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33 413 C1.w[0] = 0x38c15b0a00000000ull; 414 x_exp = x_exp + EXP_P1; 415 } 416 } else { // x > 0 417 // subtract 1 ulp (subtract 1 from the significand) 418 C1.w[0]--; 419 if (C1.w[0] == 0xffffffffffffffffull) 420 C1.w[1]--; 421 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1 422 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1 423 C1.w[0] = 0x378d8e63ffffffffull; 424 x_exp = x_exp - EXP_P1; 425 } 426 } 427 // assemble the result 428 res.w[1] = x_sign | x_exp | C1.w[1]; 429 res.w[0] = C1.w[0]; 430 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp 431 } // end x is not special and is not zero 432 BID_RETURN (res); 433 } 434 435 /***************************************************************************** 436 * BID128 nextafter 437 ****************************************************************************/ 438 439 #if DECIMAL_CALL_BY_REFERENCE 440 void 441 bid128_nextafter (UINT128 * pres, UINT128 * px, 442 UINT128 * 443 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) 444 { 445 UINT128 x = *px; 446 UINT128 y = *py; 447 UINT128 xnswp = *px; 448 UINT128 ynswp = *py; 449 #else 450 UINT128 451 bid128_nextafter (UINT128 x, 452 UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 453 _EXC_INFO_PARAM) { 454 UINT128 xnswp = x; 455 UINT128 ynswp = y; 456 #endif 457 458 UINT128 res; 459 UINT128 tmp1, tmp2, tmp3; 460 FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions 461 int res1, res2; 462 UINT64 x_exp; 463 464 465 BID_SWAP128 (x); 466 BID_SWAP128 (y); 467 // check for NaNs 468 if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) 469 || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) { 470 // x is special or y is special 471 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN 472 // if x = NaN, then res = Q (x) 473 // check first for non-canonical NaN payload 474 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 475 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) 476 && (x.w[0] > 0x38c15b09ffffffffull))) { 477 x.w[1] = x.w[1] & 0xffffc00000000000ull; 478 x.w[0] = 0x0ull; 479 } 480 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 481 // set invalid flag 482 *pfpsf |= INVALID_EXCEPTION; 483 // return quiet (x) 484 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 485 res.w[0] = x.w[0]; 486 } else { // x is QNaN 487 // return x 488 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 489 res.w[0] = x.w[0]; 490 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN 491 // set invalid flag 492 *pfpsf |= INVALID_EXCEPTION; 493 } 494 } 495 BID_RETURN (res) 496 } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN 497 // if x = NaN, then res = Q (x) 498 // check first for non-canonical NaN payload 499 if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) || 500 (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) 501 && (y.w[0] > 0x38c15b09ffffffffull))) { 502 y.w[1] = y.w[1] & 0xffffc00000000000ull; 503 y.w[0] = 0x0ull; 504 } 505 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN 506 // set invalid flag 507 *pfpsf |= INVALID_EXCEPTION; 508 // return quiet (x) 509 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16] 510 res.w[0] = y.w[0]; 511 } else { // x is QNaN 512 // return x 513 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16] 514 res.w[0] = y.w[0]; 515 } 516 BID_RETURN (res) 517 } else { // at least one is infinity 518 if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf 519 x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF); 520 x.w[0] = 0x0ull; 521 } 522 if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf 523 y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF); 524 y.w[0] = 0x0ull; 525 } 526 } 527 } 528 // neither x nor y is NaN 529 530 // if not infinity, check for non-canonical values x (treated as zero) 531 if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf 532 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11 533 // non-canonical 534 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits 535 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp; 536 x.w[0] = 0x0ull; 537 } else { // G0_G1 != 11 538 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits 539 if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull || 540 ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull 541 && x.w[0] > 0x378d8e63ffffffffull)) { 542 // x is non-canonical if coefficient is larger than 10^34 -1 543 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp; 544 x.w[0] = 0x0ull; 545 } else { // canonical 546 ; 547 } 548 } 549 } 550 // no need to check for non-canonical y 551 552 // neither x nor y is NaN 553 tmp_fpsf = *pfpsf; // save fpsf 554 #if DECIMAL_CALL_BY_REFERENCE 555 bid128_quiet_equal (&res1, &xnswp, 556 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 557 _EXC_INFO_ARG); 558 bid128_quiet_greater (&res2, &xnswp, 559 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 560 _EXC_INFO_ARG); 561 #else 562 res1 = 563 bid128_quiet_equal (xnswp, 564 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 565 _EXC_INFO_ARG); 566 res2 = 567 bid128_quiet_greater (xnswp, 568 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 569 _EXC_INFO_ARG); 570 #endif 571 *pfpsf = tmp_fpsf; // restore fpsf 572 573 if (res1) { // x = y 574 // return x with the sign of y 575 res.w[1] = 576 (x.w[1] & 0x7fffffffffffffffull) | (y. 577 w[1] & 0x8000000000000000ull); 578 res.w[0] = x.w[0]; 579 } else if (res2) { // x > y 580 #if DECIMAL_CALL_BY_REFERENCE 581 bid128_nextdown (&res, 582 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 583 _EXC_INFO_ARG); 584 #else 585 res = 586 bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG 587 _EXC_INFO_ARG); 588 #endif 589 BID_SWAP128 (res); 590 } else { // x < y 591 #if DECIMAL_CALL_BY_REFERENCE 592 bid128_nextup (&res, 593 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 594 #else 595 res = 596 bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG); 597 #endif 598 BID_SWAP128 (res); 599 } 600 // if the operand x is finite but the result is infinite, signal 601 // overflow and inexact 602 if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL) 603 && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) { 604 // set the inexact flag 605 *pfpsf |= INEXACT_EXCEPTION; 606 // set the overflow flag 607 *pfpsf |= OVERFLOW_EXCEPTION; 608 } 609 // if the result is in (-10^emin, 10^emin), and is different from the 610 // operand x, signal underflow and inexact 611 tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull; 612 tmp1.w[LOW_128W] = 0x38c15b0a00000000ull; // +100...0[34] * 10^emin 613 tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull; 614 tmp2.w[LOW_128W] = res.w[0]; 615 tmp3.w[HIGH_128W] = res.w[1]; 616 tmp3.w[LOW_128W] = res.w[0]; 617 tmp_fpsf = *pfpsf; // save fpsf 618 #if DECIMAL_CALL_BY_REFERENCE 619 bid128_quiet_greater (&res1, &tmp1, 620 &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG 621 _EXC_INFO_ARG); 622 bid128_quiet_not_equal (&res2, &xnswp, 623 &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG 624 _EXC_INFO_ARG); 625 #else 626 res1 = 627 bid128_quiet_greater (tmp1, 628 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG 629 _EXC_INFO_ARG); 630 res2 = 631 bid128_quiet_not_equal (xnswp, 632 tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG 633 _EXC_INFO_ARG); 634 #endif 635 *pfpsf = tmp_fpsf; // restore fpsf 636 if (res1 && res2) { 637 // set the inexact flag 638 *pfpsf |= INEXACT_EXCEPTION; 639 // set the underflow flag 640 *pfpsf |= UNDERFLOW_EXCEPTION; 641 } 642 BID_RETURN (res); 643 } 644