libgcc/config/libbid/bid_inline_add.h

   1 /* Copyright (C) 2007-2022 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 /*****************************************************************************
  25  *
  26  *    Helper add functions (for fma)
  27  *
  28  *    __BID_INLINE__ UINT64 get_add64(
  29  *        UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
  30  *        UINT64 sign_y, int exponent_y, UINT64 coefficient_y,
  31  *                                       int rounding_mode)
  32  *
  33  *   __BID_INLINE__ UINT64 get_add128(
  34  *                       UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
  35  *                       UINT64 sign_y, int final_exponent_y, UINT128 CY,
  36  *                       int extra_digits, int rounding_mode)
  37  *
  38  *****************************************************************************
  39  *
  40  *  Algorithm description:
  41  *
  42  *  get_add64:  same as BID64 add, but arguments are unpacked and there
  43  *                                 are no special case checks
  44  *
  45  *  get_add128: add 64-bit coefficient to 128-bit product (which contains
  46  *                                        16+extra_digits decimal digits),
  47  *                         return BID64 result
  48  *              - the exponents are compared and the two coefficients are
  49  *                properly aligned for addition/subtraction
  50  *              - multiple paths are needed
  51  *              - final result exponent is calculated and the lower term is
  52  *                      rounded first if necessary, to avoid manipulating
  53  *                      coefficients longer than 128 bits
  54  *
  55  ****************************************************************************/
  56
  57 #ifndef _INLINE_BID_ADD_H_
  58 #define _INLINE_BID_ADD_H_
  59
  60 #include "bid_internal.h"
  61
  62 #define MAX_FORMAT_DIGITS     16
  63 #define DECIMAL_EXPONENT_BIAS 398
  64 #define MASK_BINARY_EXPONENT  0x7ff0000000000000ull
  65 #define BINARY_EXPONENT_BIAS  0x3ff
  66 #define UPPER_EXPON_LIMIT     51
  67
  68 ///////////////////////////////////////////////////////////////////////
  69 //
  70 // get_add64() is essentially the same as bid_add(), except that
  71 //             the arguments are unpacked
  72 //
  73 //////////////////////////////////////////////////////////////////////
  74 __BID_INLINE__ UINT64
  75 get_add64 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
  76            UINT64 sign_y, int exponent_y, UINT64 coefficient_y,
  77            int rounding_mode, unsigned *fpsc) {
  78   UINT128 CA, CT, CT_new;
  79   UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
  80     rem_a;
  81   UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp,
  82     C64_new;
  83   int_double tempx;
  84   int exponent_a, exponent_b, diff_dec_expon;
  85   int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
  86   unsigned rmode, status;
  87
  88   // sort arguments by exponent
  89   if (exponent_x <= exponent_y) {
  90     sign_a = sign_y;
  91     exponent_a = exponent_y;
  92     coefficient_a = coefficient_y;
  93     sign_b = sign_x;
  94     exponent_b = exponent_x;
  95     coefficient_b = coefficient_x;
  96   } else {
  97     sign_a = sign_x;
  98     exponent_a = exponent_x;
  99     coefficient_a = coefficient_x;
 100     sign_b = sign_y;
 101     exponent_b = exponent_y;
 102     coefficient_b = coefficient_y;
 103   }
 104
 105   // exponent difference
 106   diff_dec_expon = exponent_a - exponent_b;
 107
 108   /* get binary coefficients of x and y */
 109
 110   //--- get number of bits in the coefficients of x and y ---
 111
 112   tempx.d = (double) coefficient_a;
 113   bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 114
 115   if (!coefficient_a) {
 116     return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode,
 117                       fpsc);
 118   }
 119   if (diff_dec_expon > MAX_FORMAT_DIGITS) {
 120     // normalize a to a 16-digit coefficient
 121
 122     scale_ca = estimate_decimal_digits[bin_expon_ca];
 123     if (coefficient_a >= power10_table_128[scale_ca].w[0])
 124       scale_ca++;
 125
 126     scale_k = 16 - scale_ca;
 127
 128     coefficient_a *= power10_table_128[scale_k].w[0];
 129
 130     diff_dec_expon -= scale_k;
 131     exponent_a -= scale_k;
 132
 133     /* get binary coefficients of x and y */
 134
 135     //--- get number of bits in the coefficients of x and y ---
 136     tempx.d = (double) coefficient_a;
 137     bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 138
 139     if (diff_dec_expon > MAX_FORMAT_DIGITS) {
 140 #ifdef SET_STATUS_FLAGS
 141       if (coefficient_b) {
 142         __set_status_flags (fpsc, INEXACT_EXCEPTION);
 143       }
 144 #endif
 145
 146 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 147 #ifndef IEEE_ROUND_NEAREST
 148       if (((rounding_mode) & 3) && coefficient_b)       // not ROUNDING_TO_NEAREST
 149       {
 150         switch (rounding_mode) {
 151         case ROUNDING_DOWN:
 152           if (sign_b) {
 153             coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
 154             if (coefficient_a < 1000000000000000ull) {
 155               exponent_a--;
 156               coefficient_a = 9999999999999999ull;
 157             } else if (coefficient_a >= 10000000000000000ull) {
 158               exponent_a++;
 159               coefficient_a = 1000000000000000ull;
 160             }
 161           }
 162           break;
 163         case ROUNDING_UP:
 164           if (!sign_b) {
 165             coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
 166             if (coefficient_a < 1000000000000000ull) {
 167               exponent_a--;
 168               coefficient_a = 9999999999999999ull;
 169             } else if (coefficient_a >= 10000000000000000ull) {
 170               exponent_a++;
 171               coefficient_a = 1000000000000000ull;
 172             }
 173           }
 174           break;
 175         default:        // RZ
 176           if (sign_a != sign_b) {
 177             coefficient_a--;
 178             if (coefficient_a < 1000000000000000ull) {
 179               exponent_a--;
 180               coefficient_a = 9999999999999999ull;
 181             }
 182           }
 183           break;
 184         }
 185       } else
 186 #endif
 187 #endif
 188         // check special case here
 189         if ((coefficient_a == 1000000000000000ull)
 190             && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
 191             && (sign_a ^ sign_b)
 192             && (coefficient_b > 5000000000000000ull)) {
 193         coefficient_a = 9999999999999999ull;
 194         exponent_a--;
 195       }
 196
 197       return get_BID64 (sign_a, exponent_a, coefficient_a,
 198                         rounding_mode, fpsc);
 199     }
 200   }
 201   // test whether coefficient_a*10^(exponent_a-exponent_b)  may exceed 2^62
 202   if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
 203     // coefficient_a*10^(exponent_a-exponent_b)<2^63
 204
 205     // multiply by 10^(exponent_a-exponent_b)
 206     coefficient_a *= power10_table_128[diff_dec_expon].w[0];
 207
 208     // sign mask
 209     sign_b = ((SINT64) sign_b) >> 63;
 210     // apply sign to coeff. of b
 211     coefficient_b = (coefficient_b + sign_b) ^ sign_b;
 212
 213     // apply sign to coefficient a
 214     sign_a = ((SINT64) sign_a) >> 63;
 215     coefficient_a = (coefficient_a + sign_a) ^ sign_a;
 216
 217     coefficient_a += coefficient_b;
 218     // get sign
 219     sign_s = ((SINT64) coefficient_a) >> 63;
 220     coefficient_a = (coefficient_a + sign_s) ^ sign_s;
 221     sign_s &= 0x8000000000000000ull;
 222
 223     // coefficient_a < 10^16 ?
 224     if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
 225 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 226 #ifndef IEEE_ROUND_NEAREST
 227       if (rounding_mode == ROUNDING_DOWN && (!coefficient_a)
 228           && sign_a != sign_b)
 229         sign_s = 0x8000000000000000ull;
 230 #endif
 231 #endif
 232       return get_BID64 (sign_s, exponent_b, coefficient_a,
 233                         rounding_mode, fpsc);
 234     }
 235     // otherwise rounding is necessary
 236
 237     // already know coefficient_a<10^19
 238     // coefficient_a < 10^17 ?
 239     if (coefficient_a < power10_table_128[17].w[0])
 240       extra_digits = 1;
 241     else if (coefficient_a < power10_table_128[18].w[0])
 242       extra_digits = 2;
 243     else
 244       extra_digits = 3;
 245
 246 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 247 #ifndef IEEE_ROUND_NEAREST
 248     rmode = rounding_mode;
 249     if (sign_s && (unsigned) (rmode - 1) < 2)
 250       rmode = 3 - rmode;
 251 #else
 252     rmode = 0;
 253 #endif
 254 #else
 255     rmode = 0;
 256 #endif
 257     coefficient_a += round_const_table[rmode][extra_digits];
 258
 259     // get P*(2^M[extra_digits])/10^extra_digits
 260     __mul_64x64_to_128 (CT, coefficient_a,
 261                         reciprocals10_64[extra_digits]);
 262
 263     // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 264     amount = short_recip_scale[extra_digits];
 265     C64 = CT.w[1] >> amount;
 266
 267   } else {
 268     // coefficient_a*10^(exponent_a-exponent_b) is large
 269     sign_s = sign_a;
 270
 271 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 272 #ifndef IEEE_ROUND_NEAREST
 273     rmode = rounding_mode;
 274     if (sign_s && (unsigned) (rmode - 1) < 2)
 275       rmode = 3 - rmode;
 276 #else
 277     rmode = 0;
 278 #endif
 279 #else
 280     rmode = 0;
 281 #endif
 282
 283     // check whether we can take faster path
 284     scale_ca = estimate_decimal_digits[bin_expon_ca];
 285
 286     sign_ab = sign_a ^ sign_b;
 287     sign_ab = ((SINT64) sign_ab) >> 63;
 288
 289     // T1 = 10^(16-diff_dec_expon)
 290     T1 = power10_table_128[16 - diff_dec_expon].w[0];
 291
 292     // get number of digits in coefficient_a
 293     //P_ca = power10_table_128[scale_ca].w[0];
 294     //P_ca_m1 = power10_table_128[scale_ca-1].w[0];
 295     if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
 296       scale_ca++;
 297       //P_ca_m1 = P_ca;
 298       //P_ca = power10_table_128[scale_ca].w[0];
 299     }
 300
 301     scale_k = 16 - scale_ca;
 302
 303     // apply sign
 304     //Ts = (T1 + sign_ab) ^ sign_ab;
 305
 306     // test range of ca
 307     //X = coefficient_a + Ts - P_ca_m1;
 308
 309     // addition
 310     saved_ca = coefficient_a - T1;
 311     coefficient_a =
 312       (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
 313     extra_digits = diff_dec_expon - scale_k;
 314
 315     // apply sign
 316     saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
 317     // add 10^16 and rounding constant
 318     coefficient_b =
 319       saved_cb + 10000000000000000ull +
 320       round_const_table[rmode][extra_digits];
 321
 322     // get P*(2^M[extra_digits])/10^extra_digits
 323     __mul_64x64_to_128 (CT, coefficient_b,
 324                         reciprocals10_64[extra_digits]);
 325
 326     // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 327     amount = short_recip_scale[extra_digits];
 328     C0_64 = CT.w[1] >> amount;
 329
 330     // result coefficient
 331     C64 = C0_64 + coefficient_a;
 332     // filter out difficult (corner) cases
 333     // the following test is equivalent to
 334     // ( (initial_coefficient_a + Ts) < P_ca &&
 335     //     (initial_coefficient_a + Ts) > P_ca_m1 ),
 336     // which ensures the number of digits in coefficient_a does not change
 337     // after adding (the appropriately scaled and rounded) coefficient_b
 338     if ((UINT64) (C64 - 1000000000000000ull - 1) >
 339         9000000000000000ull - 2) {
 340       if (C64 >= 10000000000000000ull) {
 341         // result has more than 16 digits
 342         if (!scale_k) {
 343           // must divide coeff_a by 10
 344           saved_ca = saved_ca + T1;
 345           __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
 346           //reciprocals10_64[1]);
 347           coefficient_a = CA.w[1] >> 1;
 348           rem_a =
 349             saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
 350           coefficient_a = coefficient_a - T1;
 351
 352           saved_cb +=
 353             /*90000000000000000 */ +rem_a *
 354             power10_table_128[diff_dec_expon].w[0];
 355         } else
 356           coefficient_a =
 357             (SINT64) (saved_ca - T1 -
 358                       (T1 << 3)) * (SINT64) power10_table_128[scale_k -
 359                                                               1].w[0];
 360
 361         extra_digits++;
 362         coefficient_b =
 363           saved_cb + 100000000000000000ull +
 364           round_const_table[rmode][extra_digits];
 365
 366         // get P*(2^M[extra_digits])/10^extra_digits
 367         __mul_64x64_to_128 (CT, coefficient_b,
 368                             reciprocals10_64[extra_digits]);
 369
 370         // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 371         amount = short_recip_scale[extra_digits];
 372         C0_64 = CT.w[1] >> amount;
 373
 374         // result coefficient
 375         C64 = C0_64 + coefficient_a;
 376       } else if (C64 <= 1000000000000000ull) {
 377         // less than 16 digits in result
 378         coefficient_a =
 379           (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
 380                                                         1].w[0];
 381         //extra_digits --;
 382         exponent_b--;
 383         coefficient_b =
 384           (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
 385           round_const_table[rmode][extra_digits];
 386
 387         // get P*(2^M[extra_digits])/10^extra_digits
 388         __mul_64x64_to_128 (CT_new, coefficient_b,
 389                             reciprocals10_64[extra_digits]);
 390
 391         // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 392         amount = short_recip_scale[extra_digits];
 393         C0_64 = CT_new.w[1] >> amount;
 394
 395         // result coefficient
 396         C64_new = C0_64 + coefficient_a;
 397         if (C64_new < 10000000000000000ull) {
 398           C64 = C64_new;
 399 #ifdef SET_STATUS_FLAGS
 400           CT = CT_new;
 401 #endif
 402         } else
 403           exponent_b++;
 404       }
 405
 406     }
 407
 408   }
 409
 410 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 411 #ifndef IEEE_ROUND_NEAREST
 412   if (rmode == 0)       //ROUNDING_TO_NEAREST
 413 #endif
 414     if (C64 & 1) {
 415       // check whether fractional part of initial_P/10^extra_digits
 416       // is exactly .5
 417       // this is the same as fractional part of
 418       //      (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 419
 420       // get remainder
 421       remainder_h = CT.w[1] << (64 - amount);
 422
 423       // test whether fractional part is 0
 424       if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
 425         C64--;
 426       }
 427     }
 428 #endif
 429
 430 #ifdef SET_STATUS_FLAGS
 431   status = INEXACT_EXCEPTION;
 432
 433   // get remainder
 434   remainder_h = CT.w[1] << (64 - amount);
 435
 436   switch (rmode) {
 437   case ROUNDING_TO_NEAREST:
 438   case ROUNDING_TIES_AWAY:
 439     // test whether fractional part is 0
 440     if ((remainder_h == 0x8000000000000000ull)
 441         && (CT.w[0] < reciprocals10_64[extra_digits]))
 442       status = EXACT_STATUS;
 443     break;
 444   case ROUNDING_DOWN:
 445   case ROUNDING_TO_ZERO:
 446     if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
 447       status = EXACT_STATUS;
 448     break;
 449   default:
 450     // round up
 451     __add_carry_out (tmp, carry, CT.w[0],
 452                      reciprocals10_64[extra_digits]);
 453     if ((remainder_h >> (64 - amount)) + carry >=
 454         (((UINT64) 1) << amount))
 455       status = EXACT_STATUS;
 456     break;
 457   }
 458   __set_status_flags (fpsc, status);
 459
 460 #endif
 461
 462   return get_BID64 (sign_s, exponent_b + extra_digits, C64,
 463                     rounding_mode, fpsc);
 464 }
 465
 466
 467 ///////////////////////////////////////////////////////////////////
 468 // round 128-bit coefficient and return result in BID64 format
 469 // do not worry about midpoint cases
 470 //////////////////////////////////////////////////////////////////
 471 static UINT64
 472 __bid_simple_round64_sticky (UINT64 sign, int exponent, UINT128 P,
 473                              int extra_digits, int rounding_mode,
 474                              unsigned *fpsc) {
 475   UINT128 Q_high, Q_low, C128;
 476   UINT64 C64;
 477   int amount, rmode;
 478
 479   rmode = rounding_mode;
 480 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 481 #ifndef IEEE_ROUND_NEAREST
 482   if (sign && (unsigned) (rmode - 1) < 2)
 483     rmode = 3 - rmode;
 484 #endif
 485 #endif
 486   __add_128_64 (P, P, round_const_table[rmode][extra_digits]);
 487
 488   // get P*(2^M[extra_digits])/10^extra_digits
 489   __mul_128x128_full (Q_high, Q_low, P,
 490                       reciprocals10_128[extra_digits]);
 491
 492   // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 493   amount = recip_scale[extra_digits];
 494   __shr_128 (C128, Q_high, amount);
 495
 496   C64 = __low_64 (C128);
 497
 498 #ifdef SET_STATUS_FLAGS
 499
 500   __set_status_flags (fpsc, INEXACT_EXCEPTION);
 501
 502 #endif
 503
 504   return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
 505 }
 506
 507 ///////////////////////////////////////////////////////////////////
 508 // round 128-bit coefficient and return result in BID64 format
 509 ///////////////////////////////////////////////////////////////////
 510 static UINT64
 511 __bid_full_round64 (UINT64 sign, int exponent, UINT128 P,
 512                     int extra_digits, int rounding_mode,
 513                     unsigned *fpsc) {
 514   UINT128 Q_high, Q_low, C128, Stemp, PU;
 515   UINT64 remainder_h, C64, carry, CY;
 516   int amount, amount2, rmode, status = 0;
 517
 518   if (exponent < 0) {
 519     if (exponent >= -16 && (extra_digits + exponent < 0)) {
 520       extra_digits = -exponent;
 521 #ifdef SET_STATUS_FLAGS
 522       if (extra_digits > 0) {
 523         rmode = rounding_mode;
 524         if (sign && (unsigned) (rmode - 1) < 2)
 525           rmode = 3 - rmode;
 526         __add_128_128 (PU, P,
 527                        round_const_table_128[rmode][extra_digits]);
 528         if (__unsigned_compare_gt_128
 529             (power10_table_128[extra_digits + 15], PU))
 530           status = UNDERFLOW_EXCEPTION;
 531       }
 532 #endif
 533     }
 534   }
 535
 536   if (extra_digits > 0) {
 537     exponent += extra_digits;
 538     rmode = rounding_mode;
 539     if (sign && (unsigned) (rmode - 1) < 2)
 540       rmode = 3 - rmode;
 541     __add_128_128 (P, P, round_const_table_128[rmode][extra_digits]);
 542
 543     // get P*(2^M[extra_digits])/10^extra_digits
 544     __mul_128x128_full (Q_high, Q_low, P,
 545                         reciprocals10_128[extra_digits]);
 546
 547     // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 548     amount = recip_scale[extra_digits];
 549     __shr_128_long (C128, Q_high, amount);
 550
 551     C64 = __low_64 (C128);
 552
 553 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 554 #ifndef IEEE_ROUND_NEAREST
 555     if (rmode == 0)     //ROUNDING_TO_NEAREST
 556 #endif
 557       if (C64 & 1) {
 558         // check whether fractional part of initial_P/10^extra_digits
 559         // is exactly .5
 560
 561         // get remainder
 562         amount2 = 64 - amount;
 563         remainder_h = 0;
 564         remainder_h--;
 565         remainder_h >>= amount2;
 566         remainder_h = remainder_h & Q_high.w[0];
 567
 568         if (!remainder_h
 569             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 570                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 571                     && Q_low.w[0] <
 572                     reciprocals10_128[extra_digits].w[0]))) {
 573           C64--;
 574         }
 575       }
 576 #endif
 577
 578 #ifdef SET_STATUS_FLAGS
 579     status |= INEXACT_EXCEPTION;
 580
 581     // get remainder
 582     remainder_h = Q_high.w[0] << (64 - amount);
 583
 584     switch (rmode) {
 585     case ROUNDING_TO_NEAREST:
 586     case ROUNDING_TIES_AWAY:
 587       // test whether fractional part is 0
 588       if (remainder_h == 0x8000000000000000ull
 589           && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 590               || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 591                   && Q_low.w[0] <
 592                   reciprocals10_128[extra_digits].w[0])))
 593         status = EXACT_STATUS;
 594       break;
 595     case ROUNDING_DOWN:
 596     case ROUNDING_TO_ZERO:
 597       if (!remainder_h
 598           && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 599               || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 600                   && Q_low.w[0] <
 601                   reciprocals10_128[extra_digits].w[0])))
 602         status = EXACT_STATUS;
 603       break;
 604     default:
 605       // round up
 606       __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
 607                        reciprocals10_128[extra_digits].w[0]);
 608       __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
 609                           reciprocals10_128[extra_digits].w[1], CY);
 610       if ((remainder_h >> (64 - amount)) + carry >=
 611           (((UINT64) 1) << amount))
 612         status = EXACT_STATUS;
 613     }
 614
 615     __set_status_flags (fpsc, status);
 616
 617 #endif
 618   } else {
 619     C64 = P.w[0];
 620     if (!C64) {
 621       sign = 0;
 622       if (rounding_mode == ROUNDING_DOWN)
 623         sign = 0x8000000000000000ull;
 624     }
 625   }
 626   return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
 627 }
 628
 629 /////////////////////////////////////////////////////////////////////////////////
 630 // round 192-bit coefficient (P, remainder_P) and return result in BID64 format
 631 // the lowest 64 bits (remainder_P) are used for midpoint checking only
 632 ////////////////////////////////////////////////////////////////////////////////
 633 static UINT64
 634 __bid_full_round64_remainder (UINT64 sign, int exponent, UINT128 P,
 635                               int extra_digits, UINT64 remainder_P,
 636                               int rounding_mode, unsigned *fpsc,
 637                               unsigned uf_status) {
 638   UINT128 Q_high, Q_low, C128, Stemp;
 639   UINT64 remainder_h, C64, carry, CY;
 640   int amount, amount2, rmode, status = uf_status;
 641
 642   rmode = rounding_mode;
 643   if (sign && (unsigned) (rmode - 1) < 2)
 644     rmode = 3 - rmode;
 645   if (rmode == ROUNDING_UP && remainder_P) {
 646     P.w[0]++;
 647     if (!P.w[0])
 648       P.w[1]++;
 649   }
 650
 651   if (extra_digits) {
 652     __add_128_64 (P, P, round_const_table[rmode][extra_digits]);
 653
 654     // get P*(2^M[extra_digits])/10^extra_digits
 655     __mul_128x128_full (Q_high, Q_low, P,
 656                         reciprocals10_128[extra_digits]);
 657
 658     // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 659     amount = recip_scale[extra_digits];
 660     __shr_128 (C128, Q_high, amount);
 661
 662     C64 = __low_64 (C128);
 663
 664 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 665 #ifndef IEEE_ROUND_NEAREST
 666     if (rmode == 0)     //ROUNDING_TO_NEAREST
 667 #endif
 668       if (!remainder_P && (C64 & 1)) {
 669         // check whether fractional part of initial_P/10^extra_digits
 670         // is exactly .5
 671
 672         // get remainder
 673         amount2 = 64 - amount;
 674         remainder_h = 0;
 675         remainder_h--;
 676         remainder_h >>= amount2;
 677         remainder_h = remainder_h & Q_high.w[0];
 678
 679         if (!remainder_h
 680             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 681                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 682                     && Q_low.w[0] <
 683                     reciprocals10_128[extra_digits].w[0]))) {
 684           C64--;
 685         }
 686       }
 687 #endif
 688
 689 #ifdef SET_STATUS_FLAGS
 690     status |= INEXACT_EXCEPTION;
 691
 692     if (!remainder_P) {
 693       // get remainder
 694       remainder_h = Q_high.w[0] << (64 - amount);
 695
 696       switch (rmode) {
 697       case ROUNDING_TO_NEAREST:
 698       case ROUNDING_TIES_AWAY:
 699         // test whether fractional part is 0
 700         if (remainder_h == 0x8000000000000000ull
 701             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 702                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 703                     && Q_low.w[0] <
 704                     reciprocals10_128[extra_digits].w[0])))
 705           status = EXACT_STATUS;
 706         break;
 707       case ROUNDING_DOWN:
 708       case ROUNDING_TO_ZERO:
 709         if (!remainder_h
 710             && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
 711                 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
 712                     && Q_low.w[0] <
 713                     reciprocals10_128[extra_digits].w[0])))
 714           status = EXACT_STATUS;
 715         break;
 716       default:
 717         // round up
 718         __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
 719                          reciprocals10_128[extra_digits].w[0]);
 720         __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
 721                             reciprocals10_128[extra_digits].w[1], CY);
 722         if ((remainder_h >> (64 - amount)) + carry >=
 723             (((UINT64) 1) << amount))
 724           status = EXACT_STATUS;
 725       }
 726     }
 727     __set_status_flags (fpsc, status);
 728
 729 #endif
 730   } else {
 731     C64 = P.w[0];
 732 #ifdef SET_STATUS_FLAGS
 733     if (remainder_P) {
 734       __set_status_flags (fpsc, uf_status | INEXACT_EXCEPTION);
 735     }
 736 #endif
 737   }
 738
 739   return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode,
 740                     fpsc);
 741 }
 742
 743
 744 ///////////////////////////////////////////////////////////////////
 745 // get P/10^extra_digits
 746 // result fits in 64 bits
 747 ///////////////////////////////////////////////////////////////////
 748 __BID_INLINE__ UINT64
 749 __truncate (UINT128 P, int extra_digits)
 750 // extra_digits <= 16
 751 {
 752   UINT128 Q_high, Q_low, C128;
 753   UINT64 C64;
 754   int amount;
 755
 756   // get P*(2^M[extra_digits])/10^extra_digits
 757   __mul_128x128_full (Q_high, Q_low, P,
 758                       reciprocals10_128[extra_digits]);
 759
 760   // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 761   amount = recip_scale[extra_digits];
 762   __shr_128 (C128, Q_high, amount);
 763
 764   C64 = __low_64 (C128);
 765
 766   return C64;
 767 }
 768
 769
 770 ///////////////////////////////////////////////////////////////////
 771 // return number of decimal digits in 128-bit value X
 772 ///////////////////////////////////////////////////////////////////
 773 __BID_INLINE__ int
 774 __get_dec_digits64 (UINT128 X) {
 775   int_double tempx;
 776   int digits_x, bin_expon_cx;
 777
 778   if (!X.w[1]) {
 779     //--- get number of bits in the coefficients of x and y ---
 780     tempx.d = (double) X.w[0];
 781     bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 782     // get number of decimal digits in the coeff_x
 783     digits_x = estimate_decimal_digits[bin_expon_cx];
 784     if (X.w[0] >= power10_table_128[digits_x].w[0])
 785       digits_x++;
 786     return digits_x;
 787   }
 788   tempx.d = (double) X.w[1];
 789   bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 790   // get number of decimal digits in the coeff_x
 791   digits_x = estimate_decimal_digits[bin_expon_cx + 64];
 792   if (__unsigned_compare_ge_128 (X, power10_table_128[digits_x]))
 793     digits_x++;
 794
 795   return digits_x;
 796 }
 797
 798
 799 ////////////////////////////////////////////////////////////////////////////////
 800 //
 801 // add 64-bit coefficient to 128-bit coefficient, return result in BID64 format
 802 //
 803 ////////////////////////////////////////////////////////////////////////////////
 804 __BID_INLINE__ UINT64
 805 get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
 806             UINT64 sign_y, int final_exponent_y, UINT128 CY,
 807             int extra_digits, int rounding_mode, unsigned *fpsc) {
 808   UINT128 CY_L, CX, FS, F, CT, ST, T2;
 809   UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y;
 810   SINT64 D = 0;
 811   int_double tempx;
 812   int diff_dec_expon, extra_digits2, exponent_y, status;
 813   int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode;
 814
 815   // CY has more than 16 decimal digits
 816
 817   exponent_y = final_exponent_y - extra_digits;
 818
 819 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
 820   rounding_mode = 0;
 821 #endif
 822 #ifdef IEEE_ROUND_NEAREST
 823   rounding_mode = 0;
 824 #endif
 825
 826   if (exponent_x > exponent_y) {
 827     // normalize x
 828     //--- get number of bits in the coefficients of x and y ---
 829     tempx.d = (double) coefficient_x;
 830     bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
 831     // get number of decimal digits in the coeff_x
 832     digits_x = estimate_decimal_digits[bin_expon_cx];
 833     if (coefficient_x >= power10_table_128[digits_x].w[0])
 834       digits_x++;
 835
 836     extra_dx = 16 - digits_x;
 837     coefficient_x *= power10_table_128[extra_dx].w[0];
 838     if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) {
 839       extra_dx++;
 840       coefficient_x = 10000000000000000ull;
 841     }
 842     exponent_x -= extra_dx;
 843
 844     if (exponent_x > exponent_y) {
 845
 846       // exponent_x > exponent_y
 847       diff_dec_expon = exponent_x - exponent_y;
 848
 849       if (exponent_x <= final_exponent_y + 1) {
 850         __mul_64x64_to_128 (CX, coefficient_x,
 851                             power10_table_128[diff_dec_expon].w[0]);
 852
 853         if (sign_x == sign_y) {
 854           __add_128_128 (CT, CY, CX);
 855           if ((exponent_x >
 856                final_exponent_y) /*&& (final_exponent_y>0) */ )
 857             extra_digits++;
 858           if (__unsigned_compare_ge_128
 859               (CT, power10_table_128[16 + extra_digits]))
 860             extra_digits++;
 861         } else {
 862           __sub_128_128 (CT, CY, CX);
 863           if (((SINT64) CT.w[1]) < 0) {
 864             CT.w[0] = 0 - CT.w[0];
 865             CT.w[1] = 0 - CT.w[1];
 866             if (CT.w[0])
 867               CT.w[1]--;
 868             sign_y = sign_x;
 869           } else if (!(CT.w[1] | CT.w[0])) {
 870             sign_y =
 871               (rounding_mode !=
 872                ROUNDING_DOWN) ? 0 : 0x8000000000000000ull;
 873           }
 874           if ((exponent_x + 1 >=
 875                final_exponent_y) /*&& (final_exponent_y>=0) */ ) {
 876             extra_digits = __get_dec_digits64 (CT) - 16;
 877             if (extra_digits <= 0) {
 878               if (!CT.w[0] && rounding_mode == ROUNDING_DOWN)
 879                 sign_y = 0x8000000000000000ull;
 880               return get_BID64 (sign_y, exponent_y, CT.w[0],
 881                                 rounding_mode, fpsc);
 882             }
 883           } else
 884             if (__unsigned_compare_gt_128
 885                 (power10_table_128[15 + extra_digits], CT))
 886             extra_digits--;
 887         }
 888
 889         return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits,
 890                                    rounding_mode, fpsc);
 891       }
 892       // diff_dec2+extra_digits is the number of digits to eliminate from
 893       //                           argument CY
 894       diff_dec2 = exponent_x - final_exponent_y;
 895
 896       if (diff_dec2 >= 17) {
 897 #ifndef IEEE_ROUND_NEAREST
 898 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 899         if ((rounding_mode) & 3) {
 900           switch (rounding_mode) {
 901           case ROUNDING_UP:
 902             if (!sign_y) {
 903               D = ((SINT64) (sign_x ^ sign_y)) >> 63;
 904               D = D + D + 1;
 905               coefficient_x += D;
 906             }
 907             break;
 908           case ROUNDING_DOWN:
 909             if (sign_y) {
 910               D = ((SINT64) (sign_x ^ sign_y)) >> 63;
 911               D = D + D + 1;
 912               coefficient_x += D;
 913             }
 914             break;
 915           case ROUNDING_TO_ZERO:
 916             if (sign_y != sign_x) {
 917               D = 0 - 1;
 918               coefficient_x += D;
 919             }
 920             break;
 921           }
 922           if (coefficient_x < 1000000000000000ull) {
 923             coefficient_x -= D;
 924             coefficient_x =
 925               D + (coefficient_x << 1) + (coefficient_x << 3);
 926             exponent_x--;
 927           }
 928         }
 929 #endif
 930 #endif
 931 #ifdef SET_STATUS_FLAGS
 932         if (CY.w[1] | CY.w[0])
 933           __set_status_flags (fpsc, INEXACT_EXCEPTION);
 934 #endif
 935         return get_BID64 (sign_x, exponent_x, coefficient_x,
 936                           rounding_mode, fpsc);
 937       }
 938       // here exponent_x <= 16+final_exponent_y
 939
 940       // truncate CY to 16 dec. digits
 941       CYh = __truncate (CY, extra_digits);
 942
 943       // get remainder
 944       T = power10_table_128[extra_digits].w[0];
 945       __mul_64x64_to_64 (CY0L, CYh, T);
 946
 947       remainder_y = CY.w[0] - CY0L;
 948
 949       // align coeff_x, CYh
 950       __mul_64x64_to_128 (CX, coefficient_x,
 951                           power10_table_128[diff_dec2].w[0]);
 952
 953       if (sign_x == sign_y) {
 954         __add_128_64 (CT, CX, CYh);
 955         if (__unsigned_compare_ge_128
 956             (CT, power10_table_128[16 + diff_dec2]))
 957           diff_dec2++;
 958       } else {
 959         if (remainder_y)
 960           CYh++;
 961         __sub_128_64 (CT, CX, CYh);
 962         if (__unsigned_compare_gt_128
 963             (power10_table_128[15 + diff_dec2], CT))
 964           diff_dec2--;
 965       }
 966
 967       return __bid_full_round64_remainder (sign_x, final_exponent_y, CT,
 968                                            diff_dec2, remainder_y,
 969                                            rounding_mode, fpsc, 0);
 970     }
 971   }
 972   // Here (exponent_x <= exponent_y)
 973   {
 974     diff_dec_expon = exponent_y - exponent_x;
 975
 976     if (diff_dec_expon > MAX_FORMAT_DIGITS) {
 977       rmode = rounding_mode;
 978
 979       if ((sign_x ^ sign_y)) {
 980         if (!CY.w[0])
 981           CY.w[1]--;
 982         CY.w[0]--;
 983         if (__unsigned_compare_gt_128
 984             (power10_table_128[15 + extra_digits], CY)) {
 985           if (rmode & 3) {
 986             extra_digits--;
 987             final_exponent_y--;
 988           } else {
 989             CY.w[0] = 1000000000000000ull;
 990             CY.w[1] = 0;
 991             extra_digits = 0;
 992           }
 993         }
 994       }
 995       __scale128_10 (CY, CY);
 996       extra_digits++;
 997       CY.w[0] |= 1;
 998
 999       return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY,
1000                                           extra_digits, rmode, fpsc);
1001     }
1002     // apply sign to coeff_x
1003     sign_x ^= sign_y;
1004     sign_x = ((SINT64) sign_x) >> 63;
1005     CX.w[0] = (coefficient_x + sign_x) ^ sign_x;
1006     CX.w[1] = sign_x;
1007
1008     // check whether CY (rounded to 16 digits) and CX have
1009     //                     any digits in the same position
1010     diff_dec2 = final_exponent_y - exponent_x;
1011
1012     if (diff_dec2 <= 17) {
1013       // align CY to 10^ex
1014       S = power10_table_128[diff_dec_expon].w[0];
1015       __mul_64x128_short (CY_L, S, CY);
1016
1017       __add_128_128 (ST, CY_L, CX);
1018       extra_digits2 = __get_dec_digits64 (ST) - 16;
1019       return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2,
1020                                  rounding_mode, fpsc);
1021     }
1022     // truncate CY to 16 dec. digits
1023     CYh = __truncate (CY, extra_digits);
1024
1025     // get remainder
1026     T = power10_table_128[extra_digits].w[0];
1027     __mul_64x64_to_64 (CY0L, CYh, T);
1028
1029     coefficient_y = CY.w[0] - CY0L;
1030     // add rounding constant
1031     rmode = rounding_mode;
1032     if (sign_y && (unsigned) (rmode - 1) < 2)
1033       rmode = 3 - rmode;
1034 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1035 #ifndef IEEE_ROUND_NEAREST
1036     if (!(rmode & 3))   //ROUNDING_TO_NEAREST
1037 #endif
1038 #endif
1039     {
1040       coefficient_y += round_const_table[rmode][extra_digits];
1041     }
1042     // align coefficient_y,  coefficient_x
1043     S = power10_table_128[diff_dec_expon].w[0];
1044     __mul_64x64_to_128 (F, coefficient_y, S);
1045
1046     // fraction
1047     __add_128_128 (FS, F, CX);
1048
1049 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1050 #ifndef IEEE_ROUND_NEAREST
1051     if (rmode == 0)     //ROUNDING_TO_NEAREST
1052 #endif
1053     {
1054       // rounding code, here RN_EVEN
1055       // 10^(extra_digits+diff_dec_expon)
1056       T2 = power10_table_128[diff_dec_expon + extra_digits];
1057       if (__unsigned_compare_gt_128 (FS, T2)
1058           || ((CYh & 1) && __test_equal_128 (FS, T2))) {
1059         CYh++;
1060         __sub_128_128 (FS, FS, T2);
1061       }
1062     }
1063 #endif
1064 #ifndef IEEE_ROUND_NEAREST
1065 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1066     if (rmode == 4)     //ROUNDING_TO_NEAREST
1067 #endif
1068     {
1069       // rounding code, here RN_AWAY
1070       // 10^(extra_digits+diff_dec_expon)
1071       T2 = power10_table_128[diff_dec_expon + extra_digits];
1072       if (__unsigned_compare_ge_128 (FS, T2)) {
1073         CYh++;
1074         __sub_128_128 (FS, FS, T2);
1075       }
1076     }
1077 #endif
1078 #ifndef IEEE_ROUND_NEAREST
1079 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1080     switch (rmode) {
1081     case ROUNDING_DOWN:
1082     case ROUNDING_TO_ZERO:
1083       if ((SINT64) FS.w[1] < 0) {
1084         CYh--;
1085         if (CYh < 1000000000000000ull) {
1086           CYh = 9999999999999999ull;
1087           final_exponent_y--;
1088         }
1089       } else {
1090         T2 = power10_table_128[diff_dec_expon + extra_digits];
1091         if (__unsigned_compare_ge_128 (FS, T2)) {
1092           CYh++;
1093           __sub_128_128 (FS, FS, T2);
1094         }
1095       }
1096       break;
1097     case ROUNDING_UP:
1098       if ((SINT64) FS.w[1] < 0)
1099         break;
1100       T2 = power10_table_128[diff_dec_expon + extra_digits];
1101       if (__unsigned_compare_gt_128 (FS, T2)) {
1102         CYh += 2;
1103         __sub_128_128 (FS, FS, T2);
1104       } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) {
1105         CYh++;
1106         FS.w[1] = FS.w[0] = 0;
1107       } else if (FS.w[1] | FS.w[0])
1108         CYh++;
1109       break;
1110     }
1111 #endif
1112 #endif
1113
1114 #ifdef SET_STATUS_FLAGS
1115     status = INEXACT_EXCEPTION;
1116 #ifndef IEEE_ROUND_NEAREST
1117 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1118     if (!(rmode & 3))
1119 #endif
1120 #endif
1121     {
1122       // RN modes
1123       if ((FS.w[1] ==
1124            round_const_table_128[0][diff_dec_expon + extra_digits].w[1])
1125           && (FS.w[0] ==
1126               round_const_table_128[0][diff_dec_expon +
1127                                        extra_digits].w[0]))
1128         status = EXACT_STATUS;
1129     }
1130 #ifndef IEEE_ROUND_NEAREST
1131 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1132     else if (!FS.w[1] && !FS.w[0])
1133       status = EXACT_STATUS;
1134 #endif
1135 #endif
1136
1137     __set_status_flags (fpsc, status);
1138 #endif
1139
1140     return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode,
1141                       fpsc);
1142   }
1143
1144 }
1145
1146 //////////////////////////////////////////////////////////////////////////
1147 //
1148 //  If coefficient_z is less than 16 digits long, normalize to 16 digits
1149 //
1150 /////////////////////////////////////////////////////////////////////////
1151 static UINT64
1152 BID_normalize (UINT64 sign_z, int exponent_z,
1153                UINT64 coefficient_z, UINT64 round_dir, int round_flag,
1154                int rounding_mode, unsigned *fpsc) {
1155   SINT64 D;
1156   int_double tempx;
1157   int digits_z, bin_expon, scale, rmode;
1158
1159 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1160 #ifndef IEEE_ROUND_NEAREST
1161   rmode = rounding_mode;
1162   if (sign_z && (unsigned) (rmode - 1) < 2)
1163     rmode = 3 - rmode;
1164 #else
1165   if (coefficient_z >= power10_table_128[15].w[0])
1166     return z;
1167 #endif
1168 #endif
1169
1170   //--- get number of bits in the coefficients of x and y ---
1171   tempx.d = (double) coefficient_z;
1172   bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1173   // get number of decimal digits in the coeff_x
1174   digits_z = estimate_decimal_digits[bin_expon];
1175   if (coefficient_z >= power10_table_128[digits_z].w[0])
1176     digits_z++;
1177
1178   scale = 16 - digits_z;
1179   exponent_z -= scale;
1180   if (exponent_z < 0) {
1181     scale += exponent_z;
1182     exponent_z = 0;
1183   }
1184   coefficient_z *= power10_table_128[scale].w[0];
1185
1186 #ifdef SET_STATUS_FLAGS
1187   if (round_flag) {
1188     __set_status_flags (fpsc, INEXACT_EXCEPTION);
1189     if (coefficient_z < 1000000000000000ull)
1190       __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
1191     else if ((coefficient_z == 1000000000000000ull) && !exponent_z
1192              && ((SINT64) (round_dir ^ sign_z) < 0) && round_flag
1193              && (rmode == ROUNDING_DOWN || rmode == ROUNDING_TO_ZERO))
1194       __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
1195   }
1196 #endif
1197
1198 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1199 #ifndef IEEE_ROUND_NEAREST
1200   if (round_flag && (rmode & 3)) {
1201     D = round_dir ^ sign_z;
1202
1203     if (rmode == ROUNDING_UP) {
1204       if (D >= 0)
1205         coefficient_z++;
1206     } else {
1207       if (D < 0)
1208         coefficient_z--;
1209       if (coefficient_z < 1000000000000000ull && exponent_z) {
1210         coefficient_z = 9999999999999999ull;
1211         exponent_z--;
1212       }
1213     }
1214   }
1215 #endif
1216 #endif
1217
1218   return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode,
1219                     fpsc);
1220 }
1221
1222
1223 //////////////////////////////////////////////////////////////////////////
1224 //
1225 //    0*10^ey + cz*10^ez,   ey<ez
1226 //
1227 //////////////////////////////////////////////////////////////////////////
1228
1229 __BID_INLINE__ UINT64
1230 add_zero64 (int exponent_y, UINT64 sign_z, int exponent_z,
1231             UINT64 coefficient_z, unsigned *prounding_mode,
1232             unsigned *fpsc) {
1233   int_double tempx;
1234   int bin_expon, scale_k, scale_cz;
1235   int diff_expon;
1236
1237   diff_expon = exponent_z - exponent_y;
1238
1239   tempx.d = (double) coefficient_z;
1240   bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1241   scale_cz = estimate_decimal_digits[bin_expon];
1242   if (coefficient_z >= power10_table_128[scale_cz].w[0])
1243     scale_cz++;
1244
1245   scale_k = 16 - scale_cz;
1246   if (diff_expon < scale_k)
1247     scale_k = diff_expon;
1248   coefficient_z *= power10_table_128[scale_k].w[0];
1249
1250   return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z,
1251                     *prounding_mode, fpsc);
1252 }
1253 #endif