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