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