/* Copyright (C) 2007-2023 Free Software Foundation, Inc.

This file is part of GCC.

GCC is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.

GCC is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

Under Section 7 of GPL version 3, you are granted additional
permissions described in the GCC Runtime Library Exception, version
3.1, as published by the Free Software Foundation.

You should have received a copy of the GNU General Public License and
a copy of the GCC Runtime Library Exception along with this program;
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
<http://www.gnu.org/licenses/>.  */

/*****************************************************************************
 *
 *    Helper add functions (for fma)
 *
 *    __BID_INLINE__ UINT64 get_add64(
 *        UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 
 *        UINT64 sign_y, int exponent_y, UINT64 coefficient_y, 
 *  					 int rounding_mode)
 *
 *   __BID_INLINE__ UINT64 get_add128(
 *                       UINT64 sign_x, int exponent_x, UINT64 coefficient_x, 
 *                       UINT64 sign_y, int final_exponent_y, UINT128 CY, 
 *                       int extra_digits, int rounding_mode)
 *
 *****************************************************************************
 *
 *  Algorithm description:
 *
 *  get_add64:  same as BID64 add, but arguments are unpacked and there 
 *                                 are no special case checks
 *
 *  get_add128: add 64-bit coefficient to 128-bit product (which contains 
 *                                        16+extra_digits decimal digits), 
 *                         return BID64 result
 *              - the exponents are compared and the two coefficients are 
 *                properly aligned for addition/subtraction
 *              - multiple paths are needed
 *              - final result exponent is calculated and the lower term is
 *                      rounded first if necessary, to avoid manipulating 
 *                      coefficients longer than 128 bits 
 *
 ****************************************************************************/

#ifndef _INLINE_BID_ADD_H_
#define _INLINE_BID_ADD_H_

#include "bid_internal.h"

#define MAX_FORMAT_DIGITS     16
#define DECIMAL_EXPONENT_BIAS 398
#define MASK_BINARY_EXPONENT  0x7ff0000000000000ull
#define BINARY_EXPONENT_BIAS  0x3ff
#define UPPER_EXPON_LIMIT     51

///////////////////////////////////////////////////////////////////////
//
// get_add64() is essentially the same as bid_add(), except that 
//             the arguments are unpacked
//
//////////////////////////////////////////////////////////////////////
__BID_INLINE__ UINT64
get_add64 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
	   UINT64 sign_y, int exponent_y, UINT64 coefficient_y,
	   int rounding_mode, unsigned *fpsc) {
  UINT128 CA, CT, CT_new;
  UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
    rem_a;
  UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp,
    C64_new;
  int_double tempx;
  int exponent_a, exponent_b, diff_dec_expon;
  int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
  unsigned rmode, status;

  // sort arguments by exponent
  if (exponent_x <= exponent_y) {
    sign_a = sign_y;
    exponent_a = exponent_y;
    coefficient_a = coefficient_y;
    sign_b = sign_x;
    exponent_b = exponent_x;
    coefficient_b = coefficient_x;
  } else {
    sign_a = sign_x;
    exponent_a = exponent_x;
    coefficient_a = coefficient_x;
    sign_b = sign_y;
    exponent_b = exponent_y;
    coefficient_b = coefficient_y;
  }

  // exponent difference
  diff_dec_expon = exponent_a - exponent_b;

  /* get binary coefficients of x and y */

  //--- get number of bits in the coefficients of x and y ---

  tempx.d = (double) coefficient_a;
  bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;

  if (!coefficient_a) {
    return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode,
		      fpsc);
  }
  if (diff_dec_expon > MAX_FORMAT_DIGITS) {
    // normalize a to a 16-digit coefficient

    scale_ca = estimate_decimal_digits[bin_expon_ca];
    if (coefficient_a >= power10_table_128[scale_ca].w[0])
      scale_ca++;

    scale_k = 16 - scale_ca;

    coefficient_a *= power10_table_128[scale_k].w[0];

    diff_dec_expon -= scale_k;
    exponent_a -= scale_k;

    /* get binary coefficients of x and y */

    //--- get number of bits in the coefficients of x and y ---
    tempx.d = (double) coefficient_a;
    bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;

    if (diff_dec_expon > MAX_FORMAT_DIGITS) {
#ifdef SET_STATUS_FLAGS
      if (coefficient_b) {
	__set_status_flags (fpsc, INEXACT_EXCEPTION);
      }
#endif

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
      if (((rounding_mode) & 3) && coefficient_b)	// not ROUNDING_TO_NEAREST
      {
	switch (rounding_mode) {
	case ROUNDING_DOWN:
	  if (sign_b) {
	    coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
	    if (coefficient_a < 1000000000000000ull) {
	      exponent_a--;
	      coefficient_a = 9999999999999999ull;
	    } else if (coefficient_a >= 10000000000000000ull) {
	      exponent_a++;
	      coefficient_a = 1000000000000000ull;
	    }
	  }
	  break;
	case ROUNDING_UP:
	  if (!sign_b) {
	    coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
	    if (coefficient_a < 1000000000000000ull) {
	      exponent_a--;
	      coefficient_a = 9999999999999999ull;
	    } else if (coefficient_a >= 10000000000000000ull) {
	      exponent_a++;
	      coefficient_a = 1000000000000000ull;
	    }
	  }
	  break;
	default:	// RZ
	  if (sign_a != sign_b) {
	    coefficient_a--;
	    if (coefficient_a < 1000000000000000ull) {
	      exponent_a--;
	      coefficient_a = 9999999999999999ull;
	    }
	  }
	  break;
	}
      } else
#endif
#endif
	// check special case here
	if ((coefficient_a == 1000000000000000ull)
	    && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
	    && (sign_a ^ sign_b)
	    && (coefficient_b > 5000000000000000ull)) {
	coefficient_a = 9999999999999999ull;
	exponent_a--;
      }

      return get_BID64 (sign_a, exponent_a, coefficient_a,
			rounding_mode, fpsc);
    }
  }
  // test whether coefficient_a*10^(exponent_a-exponent_b)  may exceed 2^62
  if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
    // coefficient_a*10^(exponent_a-exponent_b)<2^63

    // multiply by 10^(exponent_a-exponent_b)
    coefficient_a *= power10_table_128[diff_dec_expon].w[0];

    // sign mask
    sign_b = ((SINT64) sign_b) >> 63;
    // apply sign to coeff. of b
    coefficient_b = (coefficient_b + sign_b) ^ sign_b;

    // apply sign to coefficient a
    sign_a = ((SINT64) sign_a) >> 63;
    coefficient_a = (coefficient_a + sign_a) ^ sign_a;

    coefficient_a += coefficient_b;
    // get sign
    sign_s = ((SINT64) coefficient_a) >> 63;
    coefficient_a = (coefficient_a + sign_s) ^ sign_s;
    sign_s &= 0x8000000000000000ull;

    // coefficient_a < 10^16 ?
    if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
      if (rounding_mode == ROUNDING_DOWN && (!coefficient_a)
	  && sign_a != sign_b)
	sign_s = 0x8000000000000000ull;
#endif
#endif
      return get_BID64 (sign_s, exponent_b, coefficient_a,
			rounding_mode, fpsc);
    }
    // otherwise rounding is necessary

    // already know coefficient_a<10^19
    // coefficient_a < 10^17 ?
    if (coefficient_a < power10_table_128[17].w[0])
      extra_digits = 1;
    else if (coefficient_a < power10_table_128[18].w[0])
      extra_digits = 2;
    else
      extra_digits = 3;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    rmode = rounding_mode;
    if (sign_s && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
#else
    rmode = 0;
#endif
#else
    rmode = 0;
#endif
    coefficient_a += round_const_table[rmode][extra_digits];

    // get P*(2^M[extra_digits])/10^extra_digits
    __mul_64x64_to_128 (CT, coefficient_a,
			reciprocals10_64[extra_digits]);

    // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
    amount = short_recip_scale[extra_digits];
    C64 = CT.w[1] >> amount;

  } else {
    // coefficient_a*10^(exponent_a-exponent_b) is large
    sign_s = sign_a;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    rmode = rounding_mode;
    if (sign_s && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
#else
    rmode = 0;
#endif
#else
    rmode = 0;
#endif

    // check whether we can take faster path
    scale_ca = estimate_decimal_digits[bin_expon_ca];

    sign_ab = sign_a ^ sign_b;
    sign_ab = ((SINT64) sign_ab) >> 63;

    // T1 = 10^(16-diff_dec_expon)
    T1 = power10_table_128[16 - diff_dec_expon].w[0];

    // get number of digits in coefficient_a
    //P_ca = power10_table_128[scale_ca].w[0];
    //P_ca_m1 = power10_table_128[scale_ca-1].w[0];
    if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
      scale_ca++;
      //P_ca_m1 = P_ca;
      //P_ca = power10_table_128[scale_ca].w[0];
    }

    scale_k = 16 - scale_ca;

    // apply sign
    //Ts = (T1 + sign_ab) ^ sign_ab;

    // test range of ca
    //X = coefficient_a + Ts - P_ca_m1;

    // addition
    saved_ca = coefficient_a - T1;
    coefficient_a =
      (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
    extra_digits = diff_dec_expon - scale_k;

    // apply sign
    saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
    // add 10^16 and rounding constant
    coefficient_b =
      saved_cb + 10000000000000000ull +
      round_const_table[rmode][extra_digits];

    // get P*(2^M[extra_digits])/10^extra_digits
    __mul_64x64_to_128 (CT, coefficient_b,
			reciprocals10_64[extra_digits]);

    // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
    amount = short_recip_scale[extra_digits];
    C0_64 = CT.w[1] >> amount;

    // result coefficient 
    C64 = C0_64 + coefficient_a;
    // filter out difficult (corner) cases
    // the following test is equivalent to 
    // ( (initial_coefficient_a + Ts) < P_ca && 
    //     (initial_coefficient_a + Ts) > P_ca_m1 ), 
    // which ensures the number of digits in coefficient_a does not change 
    // after adding (the appropriately scaled and rounded) coefficient_b
    if ((UINT64) (C64 - 1000000000000000ull - 1) >
	9000000000000000ull - 2) {
      if (C64 >= 10000000000000000ull) {
	// result has more than 16 digits
	if (!scale_k) {
	  // must divide coeff_a by 10
	  saved_ca = saved_ca + T1;
	  __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
	  //reciprocals10_64[1]);
	  coefficient_a = CA.w[1] >> 1;
	  rem_a =
	    saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
	  coefficient_a = coefficient_a - T1;

	  saved_cb +=
	    /*90000000000000000 */ +rem_a *
	    power10_table_128[diff_dec_expon].w[0];
	} else
	  coefficient_a =
	    (SINT64) (saved_ca - T1 -
		      (T1 << 3)) * (SINT64) power10_table_128[scale_k -
							      1].w[0];

	extra_digits++;
	coefficient_b =
	  saved_cb + 100000000000000000ull +
	  round_const_table[rmode][extra_digits];

	// get P*(2^M[extra_digits])/10^extra_digits
	__mul_64x64_to_128 (CT, coefficient_b,
			    reciprocals10_64[extra_digits]);

	// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
	amount = short_recip_scale[extra_digits];
	C0_64 = CT.w[1] >> amount;

	// result coefficient 
	C64 = C0_64 + coefficient_a;
      } else if (C64 <= 1000000000000000ull) {
	// less than 16 digits in result
	coefficient_a =
	  (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
							1].w[0];
	//extra_digits --;
	exponent_b--;
	coefficient_b =
	  (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
	  round_const_table[rmode][extra_digits];

	// get P*(2^M[extra_digits])/10^extra_digits
	__mul_64x64_to_128 (CT_new, coefficient_b,
			    reciprocals10_64[extra_digits]);

	// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
	amount = short_recip_scale[extra_digits];
	C0_64 = CT_new.w[1] >> amount;

	// result coefficient 
	C64_new = C0_64 + coefficient_a;
	if (C64_new < 10000000000000000ull) {
	  C64 = C64_new;
#ifdef SET_STATUS_FLAGS
	  CT = CT_new;
#endif
	} else
	  exponent_b++;
      }

    }

  }

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  if (rmode == 0)	//ROUNDING_TO_NEAREST
#endif
    if (C64 & 1) {
      // check whether fractional part of initial_P/10^extra_digits 
      // is exactly .5
      // this is the same as fractional part of 
      //      (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero

      // get remainder
      remainder_h = CT.w[1] << (64 - amount);

      // test whether fractional part is 0
      if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
	C64--;
      }
    }
#endif

#ifdef SET_STATUS_FLAGS
  status = INEXACT_EXCEPTION;

  // get remainder
  remainder_h = CT.w[1] << (64 - amount);

  switch (rmode) {
  case ROUNDING_TO_NEAREST:
  case ROUNDING_TIES_AWAY:
    // test whether fractional part is 0
    if ((remainder_h == 0x8000000000000000ull)
	&& (CT.w[0] < reciprocals10_64[extra_digits]))
      status = EXACT_STATUS;
    break;
  case ROUNDING_DOWN:
  case ROUNDING_TO_ZERO:
    if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
      status = EXACT_STATUS;
    break;
  default:
    // round up
    __add_carry_out (tmp, carry, CT.w[0],
		     reciprocals10_64[extra_digits]);
    if ((remainder_h >> (64 - amount)) + carry >=
	(((UINT64) 1) << amount))
      status = EXACT_STATUS;
    break;
  }
  __set_status_flags (fpsc, status);

#endif

  return get_BID64 (sign_s, exponent_b + extra_digits, C64,
		    rounding_mode, fpsc);
}


///////////////////////////////////////////////////////////////////
// round 128-bit coefficient and return result in BID64 format
// do not worry about midpoint cases
//////////////////////////////////////////////////////////////////
static UINT64
__bid_simple_round64_sticky (UINT64 sign, int exponent, UINT128 P,
			     int extra_digits, int rounding_mode,
			     unsigned *fpsc) {
  UINT128 Q_high, Q_low, C128;
  UINT64 C64;
  int amount, rmode;

  rmode = rounding_mode;
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  if (sign && (unsigned) (rmode - 1) < 2)
    rmode = 3 - rmode;
#endif
#endif
  __add_128_64 (P, P, round_const_table[rmode][extra_digits]);

  // get P*(2^M[extra_digits])/10^extra_digits
  __mul_128x128_full (Q_high, Q_low, P,
		      reciprocals10_128[extra_digits]);

  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  amount = recip_scale[extra_digits];
  __shr_128 (C128, Q_high, amount);

  C64 = __low_64 (C128);

#ifdef SET_STATUS_FLAGS

  __set_status_flags (fpsc, INEXACT_EXCEPTION);

#endif

  return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
}

///////////////////////////////////////////////////////////////////
// round 128-bit coefficient and return result in BID64 format
///////////////////////////////////////////////////////////////////
static UINT64
__bid_full_round64 (UINT64 sign, int exponent, UINT128 P,
		    int extra_digits, int rounding_mode,
		    unsigned *fpsc) {
  UINT128 Q_high, Q_low, C128, Stemp, PU;
  UINT64 remainder_h, C64, carry, CY;
  int amount, amount2, rmode, status = 0;

  if (exponent < 0) {
    if (exponent >= -16 && (extra_digits + exponent < 0)) {
      extra_digits = -exponent;
#ifdef SET_STATUS_FLAGS
      if (extra_digits > 0) {
	rmode = rounding_mode;
	if (sign && (unsigned) (rmode - 1) < 2)
	  rmode = 3 - rmode;
	__add_128_128 (PU, P,
		       round_const_table_128[rmode][extra_digits]);
	if (__unsigned_compare_gt_128
	    (power10_table_128[extra_digits + 15], PU))
	  status = UNDERFLOW_EXCEPTION;
      }
#endif
    }
  }

  if (extra_digits > 0) {
    exponent += extra_digits;
    rmode = rounding_mode;
    if (sign && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
    __add_128_128 (P, P, round_const_table_128[rmode][extra_digits]);

    // get P*(2^M[extra_digits])/10^extra_digits
    __mul_128x128_full (Q_high, Q_low, P,
			reciprocals10_128[extra_digits]);

    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
    amount = recip_scale[extra_digits];
    __shr_128_long (C128, Q_high, amount);

    C64 = __low_64 (C128);

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    if (rmode == 0)	//ROUNDING_TO_NEAREST
#endif
      if (C64 & 1) {
	// check whether fractional part of initial_P/10^extra_digits 
	// is exactly .5

	// get remainder
	amount2 = 64 - amount;
	remainder_h = 0;
	remainder_h--;
	remainder_h >>= amount2;
	remainder_h = remainder_h & Q_high.w[0];

	if (!remainder_h
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0]))) {
	  C64--;
	}
      }
#endif

#ifdef SET_STATUS_FLAGS
    status |= INEXACT_EXCEPTION;

    // get remainder
    remainder_h = Q_high.w[0] << (64 - amount);

    switch (rmode) {
    case ROUNDING_TO_NEAREST:
    case ROUNDING_TIES_AWAY:
      // test whether fractional part is 0
      if (remainder_h == 0x8000000000000000ull
	  && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
	      || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		  && Q_low.w[0] <
		  reciprocals10_128[extra_digits].w[0])))
	status = EXACT_STATUS;
      break;
    case ROUNDING_DOWN:
    case ROUNDING_TO_ZERO:
      if (!remainder_h
	  && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
	      || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		  && Q_low.w[0] <
		  reciprocals10_128[extra_digits].w[0])))
	status = EXACT_STATUS;
      break;
    default:
      // round up
      __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
		       reciprocals10_128[extra_digits].w[0]);
      __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
			  reciprocals10_128[extra_digits].w[1], CY);
      if ((remainder_h >> (64 - amount)) + carry >=
	  (((UINT64) 1) << amount))
	status = EXACT_STATUS;
    }

    __set_status_flags (fpsc, status);

#endif
  } else {
    C64 = P.w[0];
    if (!C64) {
      sign = 0;
      if (rounding_mode == ROUNDING_DOWN)
	sign = 0x8000000000000000ull;
    }
  }
  return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
}

/////////////////////////////////////////////////////////////////////////////////
// round 192-bit coefficient (P, remainder_P) and return result in BID64 format
// the lowest 64 bits (remainder_P) are used for midpoint checking only
////////////////////////////////////////////////////////////////////////////////
static UINT64
__bid_full_round64_remainder (UINT64 sign, int exponent, UINT128 P,
			      int extra_digits, UINT64 remainder_P,
			      int rounding_mode, unsigned *fpsc,
			      unsigned uf_status) {
  UINT128 Q_high, Q_low, C128, Stemp;
  UINT64 remainder_h, C64, carry, CY;
  int amount, amount2, rmode, status = uf_status;

  rmode = rounding_mode;
  if (sign && (unsigned) (rmode - 1) < 2)
    rmode = 3 - rmode;
  if (rmode == ROUNDING_UP && remainder_P) {
    P.w[0]++;
    if (!P.w[0])
      P.w[1]++;
  }

  if (extra_digits) {
    __add_128_64 (P, P, round_const_table[rmode][extra_digits]);

    // get P*(2^M[extra_digits])/10^extra_digits
    __mul_128x128_full (Q_high, Q_low, P,
			reciprocals10_128[extra_digits]);

    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
    amount = recip_scale[extra_digits];
    __shr_128 (C128, Q_high, amount);

    C64 = __low_64 (C128);

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    if (rmode == 0)	//ROUNDING_TO_NEAREST
#endif
      if (!remainder_P && (C64 & 1)) {
	// check whether fractional part of initial_P/10^extra_digits 
	// is exactly .5

	// get remainder
	amount2 = 64 - amount;
	remainder_h = 0;
	remainder_h--;
	remainder_h >>= amount2;
	remainder_h = remainder_h & Q_high.w[0];

	if (!remainder_h
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0]))) {
	  C64--;
	}
      }
#endif

#ifdef SET_STATUS_FLAGS
    status |= INEXACT_EXCEPTION;

    if (!remainder_P) {
      // get remainder
      remainder_h = Q_high.w[0] << (64 - amount);

      switch (rmode) {
      case ROUNDING_TO_NEAREST:
      case ROUNDING_TIES_AWAY:
	// test whether fractional part is 0
	if (remainder_h == 0x8000000000000000ull
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0])))
	  status = EXACT_STATUS;
	break;
      case ROUNDING_DOWN:
      case ROUNDING_TO_ZERO:
	if (!remainder_h
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0])))
	  status = EXACT_STATUS;
	break;
      default:
	// round up
	__add_carry_out (Stemp.w[0], CY, Q_low.w[0],
			 reciprocals10_128[extra_digits].w[0]);
	__add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
			    reciprocals10_128[extra_digits].w[1], CY);
	if ((remainder_h >> (64 - amount)) + carry >=
	    (((UINT64) 1) << amount))
	  status = EXACT_STATUS;
      }
    }
    __set_status_flags (fpsc, status);

#endif
  } else {
    C64 = P.w[0];
#ifdef SET_STATUS_FLAGS
    if (remainder_P) {
      __set_status_flags (fpsc, uf_status | INEXACT_EXCEPTION);
    }
#endif
  }

  return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode,
		    fpsc);
}


///////////////////////////////////////////////////////////////////
// get P/10^extra_digits
// result fits in 64 bits
///////////////////////////////////////////////////////////////////
__BID_INLINE__ UINT64
__truncate (UINT128 P, int extra_digits)
// extra_digits <= 16
{
  UINT128 Q_high, Q_low, C128;
  UINT64 C64;
  int amount;

  // get P*(2^M[extra_digits])/10^extra_digits
  __mul_128x128_full (Q_high, Q_low, P,
		      reciprocals10_128[extra_digits]);

  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
  amount = recip_scale[extra_digits];
  __shr_128 (C128, Q_high, amount);

  C64 = __low_64 (C128);

  return C64;
}


///////////////////////////////////////////////////////////////////
// return number of decimal digits in 128-bit value X
///////////////////////////////////////////////////////////////////
__BID_INLINE__ int
__get_dec_digits64 (UINT128 X) {
  int_double tempx;
  int digits_x, bin_expon_cx;

  if (!X.w[1]) {
    //--- get number of bits in the coefficients of x and y ---
    tempx.d = (double) X.w[0];
    bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
    // get number of decimal digits in the coeff_x
    digits_x = estimate_decimal_digits[bin_expon_cx];
    if (X.w[0] >= power10_table_128[digits_x].w[0])
      digits_x++;
    return digits_x;
  }
  tempx.d = (double) X.w[1];
  bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
  // get number of decimal digits in the coeff_x
  digits_x = estimate_decimal_digits[bin_expon_cx + 64];
  if (__unsigned_compare_ge_128 (X, power10_table_128[digits_x]))
    digits_x++;

  return digits_x;
}


////////////////////////////////////////////////////////////////////////////////
//
// add 64-bit coefficient to 128-bit coefficient, return result in BID64 format
//
////////////////////////////////////////////////////////////////////////////////
__BID_INLINE__ UINT64
get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
	    UINT64 sign_y, int final_exponent_y, UINT128 CY,
	    int extra_digits, int rounding_mode, unsigned *fpsc) {
  UINT128 CY_L, CX, FS, F, CT, ST, T2;
  UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y;
  SINT64 D = 0;
  int_double tempx;
  int diff_dec_expon, extra_digits2, exponent_y, status;
  int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode;

  // CY has more than 16 decimal digits

  exponent_y = final_exponent_y - extra_digits;

#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
  rounding_mode = 0;
#endif
#ifdef IEEE_ROUND_NEAREST
  rounding_mode = 0;
#endif

  if (exponent_x > exponent_y) {
    // normalize x
    //--- get number of bits in the coefficients of x and y ---
    tempx.d = (double) coefficient_x;
    bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
    // get number of decimal digits in the coeff_x
    digits_x = estimate_decimal_digits[bin_expon_cx];
    if (coefficient_x >= power10_table_128[digits_x].w[0])
      digits_x++;

    extra_dx = 16 - digits_x;
    coefficient_x *= power10_table_128[extra_dx].w[0];
    if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) {
      extra_dx++;
      coefficient_x = 10000000000000000ull;
    }
    exponent_x -= extra_dx;

    if (exponent_x > exponent_y) {

      // exponent_x > exponent_y
      diff_dec_expon = exponent_x - exponent_y;

      if (exponent_x <= final_exponent_y + 1) {
	__mul_64x64_to_128 (CX, coefficient_x,
			    power10_table_128[diff_dec_expon].w[0]);

	if (sign_x == sign_y) {
	  __add_128_128 (CT, CY, CX);
	  if ((exponent_x >
	       final_exponent_y) /*&& (final_exponent_y>0) */ )
	    extra_digits++;
	  if (__unsigned_compare_ge_128
	      (CT, power10_table_128[16 + extra_digits]))
	    extra_digits++;
	} else {
	  __sub_128_128 (CT, CY, CX);
	  if (((SINT64) CT.w[1]) < 0) {
	    CT.w[0] = 0 - CT.w[0];
	    CT.w[1] = 0 - CT.w[1];
	    if (CT.w[0])
	      CT.w[1]--;
	    sign_y = sign_x;
	  } else if (!(CT.w[1] | CT.w[0])) {
	    sign_y =
	      (rounding_mode !=
	       ROUNDING_DOWN) ? 0 : 0x8000000000000000ull;
	  }
	  if ((exponent_x + 1 >=
	       final_exponent_y) /*&& (final_exponent_y>=0) */ ) {
	    extra_digits = __get_dec_digits64 (CT) - 16;
	    if (extra_digits <= 0) {
	      if (!CT.w[0] && rounding_mode == ROUNDING_DOWN)
		sign_y = 0x8000000000000000ull;
	      return get_BID64 (sign_y, exponent_y, CT.w[0],
				rounding_mode, fpsc);
	    }
	  } else
	    if (__unsigned_compare_gt_128
		(power10_table_128[15 + extra_digits], CT))
	    extra_digits--;
	}

	return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits,
				   rounding_mode, fpsc);
      }
      // diff_dec2+extra_digits is the number of digits to eliminate from 
      //                           argument CY
      diff_dec2 = exponent_x - final_exponent_y;

      if (diff_dec2 >= 17) {
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
	if ((rounding_mode) & 3) {
	  switch (rounding_mode) {
	  case ROUNDING_UP:
	    if (!sign_y) {
	      D = ((SINT64) (sign_x ^ sign_y)) >> 63;
	      D = D + D + 1;
	      coefficient_x += D;
	    }
	    break;
	  case ROUNDING_DOWN:
	    if (sign_y) {
	      D = ((SINT64) (sign_x ^ sign_y)) >> 63;
	      D = D + D + 1;
	      coefficient_x += D;
	    }
	    break;
	  case ROUNDING_TO_ZERO:
	    if (sign_y != sign_x) {
	      D = 0 - 1;
	      coefficient_x += D;
	    }
	    break;
	  }
	  if (coefficient_x < 1000000000000000ull) {
	    coefficient_x -= D;
	    coefficient_x =
	      D + (coefficient_x << 1) + (coefficient_x << 3);
	    exponent_x--;
	  }
	}
#endif
#endif
#ifdef SET_STATUS_FLAGS
	if (CY.w[1] | CY.w[0])
	  __set_status_flags (fpsc, INEXACT_EXCEPTION);
#endif
	return get_BID64 (sign_x, exponent_x, coefficient_x,
			  rounding_mode, fpsc);
      }
      // here exponent_x <= 16+final_exponent_y

      // truncate CY to 16 dec. digits
      CYh = __truncate (CY, extra_digits);

      // get remainder
      T = power10_table_128[extra_digits].w[0];
      __mul_64x64_to_64 (CY0L, CYh, T);

      remainder_y = CY.w[0] - CY0L;

      // align coeff_x, CYh
      __mul_64x64_to_128 (CX, coefficient_x,
			  power10_table_128[diff_dec2].w[0]);

      if (sign_x == sign_y) {
	__add_128_64 (CT, CX, CYh);
	if (__unsigned_compare_ge_128
	    (CT, power10_table_128[16 + diff_dec2]))
	  diff_dec2++;
      } else {
	if (remainder_y)
	  CYh++;
	__sub_128_64 (CT, CX, CYh);
	if (__unsigned_compare_gt_128
	    (power10_table_128[15 + diff_dec2], CT))
	  diff_dec2--;
      }

      return __bid_full_round64_remainder (sign_x, final_exponent_y, CT,
					   diff_dec2, remainder_y,
					   rounding_mode, fpsc, 0);
    }
  }
  // Here (exponent_x <= exponent_y)
  {
    diff_dec_expon = exponent_y - exponent_x;

    if (diff_dec_expon > MAX_FORMAT_DIGITS) {
      rmode = rounding_mode;

      if ((sign_x ^ sign_y)) {
	if (!CY.w[0])
	  CY.w[1]--;
	CY.w[0]--;
	if (__unsigned_compare_gt_128
	    (power10_table_128[15 + extra_digits], CY)) {
	  if (rmode & 3) {
	    extra_digits--;
	    final_exponent_y--;
	  } else {
	    CY.w[0] = 1000000000000000ull;
	    CY.w[1] = 0;
	    extra_digits = 0;
	  }
	}
      }
      __scale128_10 (CY, CY);
      extra_digits++;
      CY.w[0] |= 1;

      return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY,
					  extra_digits, rmode, fpsc);
    }
    // apply sign to coeff_x
    sign_x ^= sign_y;
    sign_x = ((SINT64) sign_x) >> 63;
    CX.w[0] = (coefficient_x + sign_x) ^ sign_x;
    CX.w[1] = sign_x;

    // check whether CY (rounded to 16 digits) and CX have 
    //                     any digits in the same position
    diff_dec2 = final_exponent_y - exponent_x;

    if (diff_dec2 <= 17) {
      // align CY to 10^ex
      S = power10_table_128[diff_dec_expon].w[0];
      __mul_64x128_short (CY_L, S, CY);

      __add_128_128 (ST, CY_L, CX);
      extra_digits2 = __get_dec_digits64 (ST) - 16;
      return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2,
				 rounding_mode, fpsc);
    }
    // truncate CY to 16 dec. digits
    CYh = __truncate (CY, extra_digits);

    // get remainder
    T = power10_table_128[extra_digits].w[0];
    __mul_64x64_to_64 (CY0L, CYh, T);

    coefficient_y = CY.w[0] - CY0L;
    // add rounding constant
    rmode = rounding_mode;
    if (sign_y && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    if (!(rmode & 3))	//ROUNDING_TO_NEAREST
#endif
#endif
    {
      coefficient_y += round_const_table[rmode][extra_digits];
    }
    // align coefficient_y,  coefficient_x
    S = power10_table_128[diff_dec_expon].w[0];
    __mul_64x64_to_128 (F, coefficient_y, S);

    // fraction
    __add_128_128 (FS, F, CX);

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    if (rmode == 0)	//ROUNDING_TO_NEAREST
#endif
    {
      // rounding code, here RN_EVEN
      // 10^(extra_digits+diff_dec_expon)
      T2 = power10_table_128[diff_dec_expon + extra_digits];
      if (__unsigned_compare_gt_128 (FS, T2)
	  || ((CYh & 1) && __test_equal_128 (FS, T2))) {
	CYh++;
	__sub_128_128 (FS, FS, T2);
      }
    }
#endif
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
    if (rmode == 4)	//ROUNDING_TO_NEAREST
#endif
    {
      // rounding code, here RN_AWAY
      // 10^(extra_digits+diff_dec_expon)
      T2 = power10_table_128[diff_dec_expon + extra_digits];
      if (__unsigned_compare_ge_128 (FS, T2)) {
	CYh++;
	__sub_128_128 (FS, FS, T2);
      }
    }
#endif
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
    switch (rmode) {
    case ROUNDING_DOWN:
    case ROUNDING_TO_ZERO:
      if ((SINT64) FS.w[1] < 0) {
	CYh--;
	if (CYh < 1000000000000000ull) {
	  CYh = 9999999999999999ull;
	  final_exponent_y--;
	}
      } else {
	T2 = power10_table_128[diff_dec_expon + extra_digits];
	if (__unsigned_compare_ge_128 (FS, T2)) {
	  CYh++;
	  __sub_128_128 (FS, FS, T2);
	}
      }
      break;
    case ROUNDING_UP:
      if ((SINT64) FS.w[1] < 0)
	break;
      T2 = power10_table_128[diff_dec_expon + extra_digits];
      if (__unsigned_compare_gt_128 (FS, T2)) {
	CYh += 2;
	__sub_128_128 (FS, FS, T2);
      } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) {
	CYh++;
	FS.w[1] = FS.w[0] = 0;
      } else if (FS.w[1] | FS.w[0])
	CYh++;
      break;
    }
#endif
#endif

#ifdef SET_STATUS_FLAGS
    status = INEXACT_EXCEPTION;
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
    if (!(rmode & 3))
#endif
#endif
    {
      // RN modes
      if ((FS.w[1] ==
	   round_const_table_128[0][diff_dec_expon + extra_digits].w[1])
	  && (FS.w[0] ==
	      round_const_table_128[0][diff_dec_expon +
				       extra_digits].w[0]))
	status = EXACT_STATUS;
    }
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
    else if (!FS.w[1] && !FS.w[0])
      status = EXACT_STATUS;
#endif
#endif

    __set_status_flags (fpsc, status);
#endif

    return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode,
		      fpsc);
  }

}

//////////////////////////////////////////////////////////////////////////
//
//  If coefficient_z is less than 16 digits long, normalize to 16 digits
//
/////////////////////////////////////////////////////////////////////////
static UINT64
BID_normalize (UINT64 sign_z, int exponent_z,
	       UINT64 coefficient_z, UINT64 round_dir, int round_flag,
	       int rounding_mode, unsigned *fpsc) {
  SINT64 D;
  int_double tempx;
  int digits_z, bin_expon, scale, rmode;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  rmode = rounding_mode;
  if (sign_z && (unsigned) (rmode - 1) < 2)
    rmode = 3 - rmode;
#else
  if (coefficient_z >= power10_table_128[15].w[0])
    return z;
#endif
#endif

  //--- get number of bits in the coefficients of x and y ---
  tempx.d = (double) coefficient_z;
  bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
  // get number of decimal digits in the coeff_x
  digits_z = estimate_decimal_digits[bin_expon];
  if (coefficient_z >= power10_table_128[digits_z].w[0])
    digits_z++;

  scale = 16 - digits_z;
  exponent_z -= scale;
  if (exponent_z < 0) {
    scale += exponent_z;
    exponent_z = 0;
  }
  coefficient_z *= power10_table_128[scale].w[0];

#ifdef SET_STATUS_FLAGS
  if (round_flag) {
    __set_status_flags (fpsc, INEXACT_EXCEPTION);
    if (coefficient_z < 1000000000000000ull)
      __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
    else if ((coefficient_z == 1000000000000000ull) && !exponent_z
	     && ((SINT64) (round_dir ^ sign_z) < 0) && round_flag
	     && (rmode == ROUNDING_DOWN || rmode == ROUNDING_TO_ZERO))
      __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
  }
#endif

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  if (round_flag && (rmode & 3)) {
    D = round_dir ^ sign_z;

    if (rmode == ROUNDING_UP) {
      if (D >= 0)
	coefficient_z++;
    } else {
      if (D < 0)
	coefficient_z--;
      if (coefficient_z < 1000000000000000ull && exponent_z) {
	coefficient_z = 9999999999999999ull;
	exponent_z--;
      }
    }
  }
#endif
#endif

  return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode,
		    fpsc);
}


//////////////////////////////////////////////////////////////////////////
//
//    0*10^ey + cz*10^ez,   ey<ez  
//
//////////////////////////////////////////////////////////////////////////

__BID_INLINE__ UINT64
add_zero64 (int exponent_y, UINT64 sign_z, int exponent_z,
	    UINT64 coefficient_z, unsigned *prounding_mode,
	    unsigned *fpsc) {
  int_double tempx;
  int bin_expon, scale_k, scale_cz;
  int diff_expon;

  diff_expon = exponent_z - exponent_y;

  tempx.d = (double) coefficient_z;
  bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
  scale_cz = estimate_decimal_digits[bin_expon];
  if (coefficient_z >= power10_table_128[scale_cz].w[0])
    scale_cz++;

  scale_k = 16 - scale_cz;
  if (diff_expon < scale_k)
    scale_k = diff_expon;
  coefficient_z *= power10_table_128[scale_k].w[0];

  return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z,
		    *prounding_mode, fpsc);
}
#endif