[thirdparty/gcc.git] / libgcc / config / libbid / bid64_mul.c

/* Copyright (C) 2007-2024 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/>.  */

/*****************************************************************************
 *    BID64 multiply
 *****************************************************************************
 *
 *  Algorithm description:
 *
 *  if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed
 *       below 16)
 *      return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias,
 *                     coefficient_x*coefficient_y)
 *  else
 *      get long product: coefficient_x*coefficient_y
 *      determine number of digits to round off (extra_digits)
 *      rounding is performed as a 128x128-bit multiplication by 
 *         2^M[extra_digits]/10^extra_digits, followed by a shift
 *         M[extra_digits] is sufficiently large for required accuracy 
 *
 ****************************************************************************/

#include "bid_internal.h"

#if DECIMAL_CALL_BY_REFERENCE

void
bid64_mul (UINT64 * pres, UINT64 * px,
	   UINT64 *
	   py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	   _EXC_INFO_PARAM) {
  UINT64 x, y;
#else

UINT64
bid64_mul (UINT64 x,
	   UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
	   _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
  UINT128 P, PU, C128, Q_high, Q_low, Stemp;
  UINT64 sign_x, sign_y, coefficient_x, coefficient_y;
  UINT64 C64, remainder_h, carry, CY, res;
  UINT64 valid_x, valid_y;
  int_double tempx, tempy;
  int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
    bin_expon_product;
  int rmode, digits_p, bp, amount, amount2, final_exponent, round_up;
  unsigned status, uf_status;

#if DECIMAL_CALL_BY_REFERENCE
#if !DECIMAL_GLOBAL_ROUNDING
  _IDEC_round rnd_mode = *prnd_mode;
#endif
  x = *px;
  y = *py;
#endif

  valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);

  // unpack arguments, check for NaN or Infinity
  if (!valid_x) {

#ifdef SET_STATUS_FLAGS
    if ((y & SNAN_MASK64) == SNAN_MASK64)	// y is sNaN
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
    // x is Inf. or NaN

    // test if x is NaN
    if ((x & NAN_MASK64) == NAN_MASK64) {
#ifdef SET_STATUS_FLAGS
      if ((x & SNAN_MASK64) == SNAN_MASK64)	// sNaN
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
      BID_RETURN (coefficient_x & QUIET_MASK64);
    }
    // x is Infinity?
    if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
      // check if y is 0
      if (((y & INFINITY_MASK64) != INFINITY_MASK64)
	  && !coefficient_y) {
#ifdef SET_STATUS_FLAGS
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	// y==0 , return NaN
	BID_RETURN (NAN_MASK64);
      }
      // check if y is NaN
      if ((y & NAN_MASK64) == NAN_MASK64)
	// y==NaN , return NaN
	BID_RETURN (coefficient_y & QUIET_MASK64);
      // otherwise return +/-Inf
      BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
    }
    // x is 0
    if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
      if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
	exponent_y = ((UINT32) (y >> 51)) & 0x3ff;
      else
	exponent_y = ((UINT32) (y >> 53)) & 0x3ff;
      sign_y = y & 0x8000000000000000ull;

      exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
      if (exponent_x > DECIMAL_MAX_EXPON_64)
	exponent_x = DECIMAL_MAX_EXPON_64;
      else if (exponent_x < 0)
	exponent_x = 0;
      BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
    }
  }
  if (!valid_y) {
    // y is Inf. or NaN

    // test if y is NaN
    if ((y & NAN_MASK64) == NAN_MASK64) {
#ifdef SET_STATUS_FLAGS
      if ((y & SNAN_MASK64) == SNAN_MASK64)	// sNaN
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
      BID_RETURN (coefficient_y & QUIET_MASK64);
    }
    // y is Infinity?
    if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
      // check if x is 0
      if (!coefficient_x) {
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	// x==0, return NaN
	BID_RETURN (NAN_MASK64);
      }
      // otherwise return +/-Inf
      BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
    }
    // y is 0
    exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
    if (exponent_x > DECIMAL_MAX_EXPON_64)
      exponent_x = DECIMAL_MAX_EXPON_64;
    else if (exponent_x < 0)
      exponent_x = 0;
    BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
  }
  //--- get number of bits in the coefficients of x and y ---
  // version 2 (original)
  tempx.d = (double) coefficient_x;
  bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52);
  tempy.d = (double) coefficient_y;
  bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52);

  // magnitude estimate for coefficient_x*coefficient_y is 
  //        2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
  bin_expon_product = bin_expon_cx + bin_expon_cy;

  // check if coefficient_x*coefficient_y<2^(10*k+3)
  // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
  if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {
    //  easy multiply
    C64 = coefficient_x * coefficient_y;

    res =
      get_BID64_small_mantissa (sign_x ^ sign_y,
				exponent_x + exponent_y -
				DECIMAL_EXPONENT_BIAS, C64, rnd_mode,
				pfpsf);
    BID_RETURN (res);
  } else {
    uf_status = 0;
    // get 128-bit product: coefficient_x*coefficient_y
    __mul_64x64_to_128 (P, coefficient_x, coefficient_y);

    // tighten binary range of P:  leading bit is 2^bp
    // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
    bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;

    __tight_bin_range_128 (bp, P, bin_expon_product);

    // get number of decimal digits in the product
    digits_p = estimate_decimal_digits[bp];
    if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P)))
      digits_p++;	// if power10_table_128[digits_p] <= P

    // determine number of decimal digits to be rounded out
    extra_digits = digits_p - MAX_FORMAT_DIGITS;
    final_exponent =
      exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    rmode = rnd_mode;
    if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
#else
    rmode = 0;
#endif
#else
    rmode = 0;
#endif

    round_up = 0;
    if (((unsigned) final_exponent) >= 3 * 256) {
      if (final_exponent < 0) {
	// underflow
	if (final_exponent + 16 < 0) {
	  res = sign_x ^ sign_y;
	  __set_status_flags (pfpsf,
			      UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
	  if (rmode == ROUNDING_UP)
	    res |= 1;
	  BID_RETURN (res);
	}

	uf_status = UNDERFLOW_EXCEPTION;
	if (final_exponent == -1) {
	  __add_128_64 (PU, P, round_const_table[rmode][extra_digits]);
	  if (__unsigned_compare_ge_128
	      (PU, power10_table_128[extra_digits + 16]))
	    uf_status = 0;
	}
	extra_digits -= final_exponent;
	final_exponent = 0;

	if (extra_digits > 17) {
	  __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]);

	  amount = recip_scale[16];
	  __shr_128 (P, Q_high, amount);

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

	  extra_digits -= 16;
	  if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1]
			      || (Q_low.w[1] ==
				  reciprocals10_128[16].w[1]
				  && Q_low.w[0] >=
				  reciprocals10_128[16].w[0]))) {
	    round_up = 1;
	    __set_status_flags (pfpsf,
				UNDERFLOW_EXCEPTION |
				INEXACT_EXCEPTION);
	    P.w[0] = (P.w[0] << 3) + (P.w[0] << 1);
	    P.w[0] |= 1;
	    extra_digits++;
	  }
	}
      } else {
	res =
	  fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
				   1000000000000000ull, rnd_mode,
				   pfpsf);
	BID_RETURN (res);
      }
    }


    if (extra_digits > 0) {
      // will divide by 10^(digits_p - 16)

      // add a constant to P, depending on rounding mode
      // 0.5*10^(digits_p - 16) for round-to-nearest
      __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 ((C64 & 1) && !round_up) {
	  // 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 = Q_high.w[0] << (64 - amount);

	  // test whether fractional part is 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 | uf_status;

      // 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 (pfpsf, status);
#endif

      // convert to BID and return
      res =
	fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64,
				 rmode, pfpsf);
      BID_RETURN (res);
    }
    // go to convert_format and exit
    C64 = __low_64 (P);
    res =
      get_BID64 (sign_x ^ sign_y,
		 exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
		 rmode, pfpsf);
    BID_RETURN (res);
  }
}
Commit	Line	Data
a945c346	1	/* Copyright (C) 2007-2024 Free Software Foundation, Inc.
200359e8 L	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
748086b7	7	Software Foundation; either version 3, or (at your option) any later
200359e8 L	8	version.
200359e8 L	9
200359e8 L	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
748086b7 JJ	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/>. */
200359e8 L	23
	24	/*****************************************************************************
	25	* BID64 multiply
	26	*****************************************************************************
	27	*
	28	* Algorithm description:
	29	*
	30	* if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed
	31	* below 16)
	32	* return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias,
	33	* coefficient_x*coefficient_y)
	34	* else
	35	* get long product: coefficient_x*coefficient_y
	36	* determine number of digits to round off (extra_digits)
	37	* rounding is performed as a 128x128-bit multiplication by
	38	* 2^M[extra_digits]/10^extra_digits, followed by a shift
	39	* M[extra_digits] is sufficiently large for required accuracy
	40	*
	41	****************************************************************************/
	42
	43	#include "bid_internal.h"
	44
	45	#if DECIMAL_CALL_BY_REFERENCE
	46
	47	void
b2a00c89	48	bid64_mul (UINT64 * pres, UINT64 * px,
200359e8 L	49	UINT64 *
	50	py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	51	_EXC_INFO_PARAM) {
	52	UINT64 x, y;
	53	#else
	54
	55	UINT64
b2a00c89	56	bid64_mul (UINT64 x,
200359e8 L	57	UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
	58	_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
	59	#endif
	60	UINT128 P, PU, C128, Q_high, Q_low, Stemp;
	61	UINT64 sign_x, sign_y, coefficient_x, coefficient_y;
	62	UINT64 C64, remainder_h, carry, CY, res;
	63	UINT64 valid_x, valid_y;
	64	int_double tempx, tempy;
b2a00c89	65	int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
200359e8 L	66	bin_expon_product;
	67	int rmode, digits_p, bp, amount, amount2, final_exponent, round_up;
	68	unsigned status, uf_status;
	69
	70	#if DECIMAL_CALL_BY_REFERENCE
	71	#if !DECIMAL_GLOBAL_ROUNDING
	72	_IDEC_round rnd_mode = *prnd_mode;
	73	#endif
	74	x = *px;
	75	y = *py;
	76	#endif
	77
	78	valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
	79	valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
	80
	81	// unpack arguments, check for NaN or Infinity
	82	if (!valid_x) {
	83
	84	#ifdef SET_STATUS_FLAGS
	85	if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
	86	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	87	#endif
	88	// x is Inf. or NaN
	89
	90	// test if x is NaN
	91	if ((x & NAN_MASK64) == NAN_MASK64) {
	92	#ifdef SET_STATUS_FLAGS
	93	if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
	94	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	95	#endif
b2a00c89	96	BID_RETURN (coefficient_x & QUIET_MASK64);
200359e8 L	97	}
	98	// x is Infinity?
	99	if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
	100	// check if y is 0
b2a00c89 L	101	if (((y & INFINITY_MASK64) != INFINITY_MASK64)
b2a00c89 L	102	&& !coefficient_y) {
200359e8 L	103	#ifdef SET_STATUS_FLAGS
	104	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	105	#endif
	106	// y==0 , return NaN
	107	BID_RETURN (NAN_MASK64);
	108	}
	109	// check if y is NaN
	110	if ((y & NAN_MASK64) == NAN_MASK64)
	111	// y==NaN , return NaN
b2a00c89	112	BID_RETURN (coefficient_y & QUIET_MASK64);
200359e8 L	113	// otherwise return +/-Inf
	114	BID_RETURN (((x ^ y) & 0x8000000000000000ull) \| INFINITY_MASK64);
	115	}
	116	// x is 0
	117	if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
	118	if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
	119	exponent_y = ((UINT32) (y >> 51)) & 0x3ff;
	120	else
	121	exponent_y = ((UINT32) (y >> 53)) & 0x3ff;
	122	sign_y = y & 0x8000000000000000ull;
	123
	124	exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
	125	if (exponent_x > DECIMAL_MAX_EXPON_64)
	126	exponent_x = DECIMAL_MAX_EXPON_64;
	127	else if (exponent_x < 0)
	128	exponent_x = 0;
	129	BID_RETURN ((sign_x ^ sign_y) \| (((UINT64) exponent_x) << 53));
	130	}
	131	}
	132	if (!valid_y) {
	133	// y is Inf. or NaN
	134
	135	// test if y is NaN
	136	if ((y & NAN_MASK64) == NAN_MASK64) {
	137	#ifdef SET_STATUS_FLAGS
	138	if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN
	139	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	140	#endif
b2a00c89	141	BID_RETURN (coefficient_y & QUIET_MASK64);
200359e8 L	142	}
	143	// y is Infinity?
	144	if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
	145	// check if x is 0
b2a00c89	146	if (!coefficient_x) {
200359e8 L	147	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	148	// x==0, return NaN
	149	BID_RETURN (NAN_MASK64);
	150	}
	151	// otherwise return +/-Inf
	152	BID_RETURN (((x ^ y) & 0x8000000000000000ull) \| INFINITY_MASK64);
	153	}
	154	// y is 0
	155	exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
	156	if (exponent_x > DECIMAL_MAX_EXPON_64)
	157	exponent_x = DECIMAL_MAX_EXPON_64;
	158	else if (exponent_x < 0)
	159	exponent_x = 0;
	160	BID_RETURN ((sign_x ^ sign_y) \| (((UINT64) exponent_x) << 53));
	161	}
	162	//--- get number of bits in the coefficients of x and y ---
	163	// version 2 (original)
	164	tempx.d = (double) coefficient_x;
	165	bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52);
	166	tempy.d = (double) coefficient_y;
	167	bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52);
	168
	169	// magnitude estimate for coefficient_x*coefficient_y is
	170	// 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
	171	bin_expon_product = bin_expon_cx + bin_expon_cy;
	172
	173	// check if coefficient_xcoefficient_y<2^(10k+3)
	174	// equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
	175	if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {
	176	// easy multiply
	177	C64 = coefficient_x * coefficient_y;
	178
	179	res =
	180	get_BID64_small_mantissa (sign_x ^ sign_y,
	181	exponent_x + exponent_y -
	182	DECIMAL_EXPONENT_BIAS, C64, rnd_mode,
	183	pfpsf);
	184	BID_RETURN (res);
	185	} else {
	186	uf_status = 0;
	187	// get 128-bit product: coefficient_x*coefficient_y
	188	__mul_64x64_to_128 (P, coefficient_x, coefficient_y);
	189
	190	// tighten binary range of P: leading bit is 2^bp
	191	// unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
	192	bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;
	193
	194	__tight_bin_range_128 (bp, P, bin_expon_product);
	195
	196	// get number of decimal digits in the product
b2a00c89 L	197	digits_p = estimate_decimal_digits[bp];
	198	if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P)))
	199	digits_p++; // if power10_table_128[digits_p] <= P
200359e8 L	200
	201	// determine number of decimal digits to be rounded out
	202	extra_digits = digits_p - MAX_FORMAT_DIGITS;
	203	final_exponent =
	204	exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;
	205
	206	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
	207	#ifndef IEEE_ROUND_NEAREST
	208	rmode = rnd_mode;
	209	if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
	210	rmode = 3 - rmode;
	211	#else
	212	rmode = 0;
	213	#endif
	214	#else
	215	rmode = 0;
	216	#endif
	217
	218	round_up = 0;
	219	if (((unsigned) final_exponent) >= 3 * 256) {
	220	if (final_exponent < 0) {
	221	// underflow
	222	if (final_exponent + 16 < 0) {
	223	res = sign_x ^ sign_y;
	224	__set_status_flags (pfpsf,
	225	UNDERFLOW_EXCEPTION \| INEXACT_EXCEPTION);
	226	if (rmode == ROUNDING_UP)
	227	res \|= 1;
	228	BID_RETURN (res);
	229	}
	230
	231	uf_status = UNDERFLOW_EXCEPTION;
	232	if (final_exponent == -1) {
b2a00c89	233	__add_128_64 (PU, P, round_const_table[rmode][extra_digits]);
200359e8	234	if (__unsigned_compare_ge_128
b2a00c89	235	(PU, power10_table_128[extra_digits + 16]))
200359e8 L	236	uf_status = 0;
	237	}
	238	extra_digits -= final_exponent;
	239	final_exponent = 0;
	240
	241	if (extra_digits > 17) {
b2a00c89	242	__mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]);
200359e8	243
b2a00c89	244	amount = recip_scale[16];
200359e8 L	245	__shr_128 (P, Q_high, amount);
	246
	247	// get sticky bits
	248	amount2 = 64 - amount;
	249	remainder_h = 0;
	250	remainder_h--;
	251	remainder_h >>= amount2;
	252	remainder_h = remainder_h & Q_high.w[0];
	253
	254	extra_digits -= 16;
b2a00c89	255	if (remainder_h \|\| (Q_low.w[1] > reciprocals10_128[16].w[1]
200359e8	256	\|\| (Q_low.w[1] ==
b2a00c89	257	reciprocals10_128[16].w[1]
200359e8	258	&& Q_low.w[0] >=
b2a00c89	259	reciprocals10_128[16].w[0]))) {
200359e8 L	260	round_up = 1;
	261	__set_status_flags (pfpsf,
	262	UNDERFLOW_EXCEPTION \|
	263	INEXACT_EXCEPTION);
	264	P.w[0] = (P.w[0] << 3) + (P.w[0] << 1);
	265	P.w[0] \|= 1;
	266	extra_digits++;
	267	}
	268	}
	269	} else {
	270	res =
	271	fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
	272	1000000000000000ull, rnd_mode,
	273	pfpsf);
	274	BID_RETURN (res);
	275	}
	276	}
	277
	278
	279	if (extra_digits > 0) {
	280	// will divide by 10^(digits_p - 16)
	281
	282	// add a constant to P, depending on rounding mode
	283	// 0.5*10^(digits_p - 16) for round-to-nearest
b2a00c89	284	__add_128_64 (P, P, round_const_table[rmode][extra_digits]);
200359e8 L	285
	286	// get P*(2^M[extra_digits])/10^extra_digits
	287	__mul_128x128_full (Q_high, Q_low, P,
b2a00c89	288	reciprocals10_128[extra_digits]);
200359e8 L	289
200359e8 L	290	// now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
b2a00c89	291	amount = recip_scale[extra_digits];
200359e8 L	292	__shr_128 (C128, Q_high, amount);
	293
	294	C64 = __low_64 (C128);
	295
	296	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
	297	#ifndef IEEE_ROUND_NEAREST
	298	if (rmode == 0) //ROUNDING_TO_NEAREST
	299	#endif
	300	if ((C64 & 1) && !round_up) {
	301	// check whether fractional part of initial_P/10^extra_digits
b2a00c89	302	// is exactly .5
200359e8 L	303	// this is the same as fractional part of
	304	// (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
	305
	306	// get remainder
	307	remainder_h = Q_high.w[0] << (64 - amount);
	308
	309	// test whether fractional part is 0
	310	if (!remainder_h
b2a00c89 L	311	&& (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
b2a00c89 L	312	\|\| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
200359e8	313	&& Q_low.w[0] <
b2a00c89	314	reciprocals10_128[extra_digits].w[0]))) {
200359e8 L	315	C64--;
	316	}
	317	}
	318	#endif
	319
	320	#ifdef SET_STATUS_FLAGS
	321	status = INEXACT_EXCEPTION \| uf_status;
	322
	323	// get remainder
	324	remainder_h = Q_high.w[0] << (64 - amount);
	325
	326	switch (rmode) {
	327	case ROUNDING_TO_NEAREST:
	328	case ROUNDING_TIES_AWAY:
	329	// test whether fractional part is 0
	330	if (remainder_h == 0x8000000000000000ull
b2a00c89 L	331	&& (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
b2a00c89 L	332	\|\| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
200359e8	333	&& Q_low.w[0] <
b2a00c89	334	reciprocals10_128[extra_digits].w[0])))
200359e8 L	335	status = EXACT_STATUS;
	336	break;
	337	case ROUNDING_DOWN:
	338	case ROUNDING_TO_ZERO:
	339	if (!remainder_h
b2a00c89 L	340	&& (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
b2a00c89 L	341	\|\| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
200359e8	342	&& Q_low.w[0] <
b2a00c89	343	reciprocals10_128[extra_digits].w[0])))
200359e8 L	344	status = EXACT_STATUS;
	345	break;
	346	default:
	347	// round up
	348	__add_carry_out (Stemp.w[0], CY, Q_low.w[0],
b2a00c89	349	reciprocals10_128[extra_digits].w[0]);
200359e8	350	__add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
b2a00c89	351	reciprocals10_128[extra_digits].w[1], CY);
200359e8 L	352	if ((remainder_h >> (64 - amount)) + carry >=
	353	(((UINT64) 1) << amount))
	354	status = EXACT_STATUS;
	355	}
	356
	357	__set_status_flags (pfpsf, status);
	358	#endif
	359
	360	// convert to BID and return
	361	res =
	362	fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64,
	363	rmode, pfpsf);
	364	BID_RETURN (res);
	365	}
	366	// go to convert_format and exit
	367	C64 = __low_64 (P);
	368	res =
	369	get_BID64 (sign_x ^ sign_y,
	370	exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
	371	rmode, pfpsf);
	372	BID_RETURN (res);
	373	}
	374	}