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

/* Copyright (C) 2007-2022 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 fma
 *****************************************************************************
 *
 *  Algorithm description:
 *
 *  if multiplication is guranteed exact (short coefficients)
 *     call the unpacked arg. equivalent of bid64_add(x*y, z)
 *  else 
 *     get full coefficient_x*coefficient_y product
 *     call subroutine to perform addition of 64-bit argument 
 *                                         to 128-bit product
 *
 ****************************************************************************/

#include "bid_inline_add.h"

#if DECIMAL_CALL_BY_REFERENCE
extern void bid64_mul (UINT64 * pres, UINT64 * px,
		       UINT64 *
		       py _RND_MODE_PARAM _EXC_FLAGS_PARAM
		       _EXC_MASKS_PARAM _EXC_INFO_PARAM);
#else

extern UINT64 bid64_mul (UINT64 x,
			 UINT64 y _RND_MODE_PARAM
			 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
			 _EXC_INFO_PARAM);
#endif

#if DECIMAL_CALL_BY_REFERENCE

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

UINT64
bid64_fma (UINT64 x, UINT64 y,
	   UINT64 z _RND_MODE_PARAM _EXC_FLAGS_PARAM
	   _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
  UINT128 P, PU, CT, CZ;
  UINT64 sign_x, sign_y, coefficient_x, coefficient_y, sign_z,
    coefficient_z;
  UINT64 C64, remainder_y, res;
  UINT64 CYh, CY0L, T, valid_x, valid_y, valid_z;
  int_double tempx, tempy;
  int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
    bin_expon_product, rmode;
  int digits_p, bp, final_exponent, exponent_z, digits_z, ez, ey,
    scale_z, uf_status;

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

  valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
  valid_z = unpack_BID64 (&sign_z, &exponent_z, &coefficient_z, z);

  // unpack arguments, check for NaN, Infinity, or 0
  if (!valid_x || !valid_y || !valid_z) {

    if ((y & MASK_NAN) == MASK_NAN) {	// y is NAN
      // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
      // check first for non-canonical NaN payload
      y = y & 0xfe03ffffffffffffull;	// clear G6-G12
      if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
	y = y & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
      }
      if ((y & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
	// set invalid flag
	*pfpsf |= INVALID_EXCEPTION;
	// return quiet (y)
	res = y & 0xfdffffffffffffffull;
      } else {	// y is QNaN
	// return y
	res = y;
	// if z = SNaN or x = SNaN signal invalid exception
	if ((z & MASK_SNAN) == MASK_SNAN
	    || (x & MASK_SNAN) == MASK_SNAN) {
	  // set invalid flag
	  *pfpsf |= INVALID_EXCEPTION;
	}
      }
      BID_RETURN (res)
    } else if ((z & MASK_NAN) == MASK_NAN) {	// z is NAN
      // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
      // check first for non-canonical NaN payload
      z = z & 0xfe03ffffffffffffull;	// clear G6-G12
      if ((z & 0x0003ffffffffffffull) > 999999999999999ull) {
	z = z & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
      }
      if ((z & MASK_SNAN) == MASK_SNAN) {	// z is SNAN
	// set invalid flag
	*pfpsf |= INVALID_EXCEPTION;
	// return quiet (z)
	res = z & 0xfdffffffffffffffull;
      } else {	// z is QNaN
	// return z
	res = z;
	// if x = SNaN signal invalid exception
	if ((x & MASK_SNAN) == MASK_SNAN) {
	  // set invalid flag
	  *pfpsf |= INVALID_EXCEPTION;
	}
      }
      BID_RETURN (res)
    } else if ((x & MASK_NAN) == MASK_NAN) {	// x is NAN
      // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
      // check first for non-canonical NaN payload
      x = x & 0xfe03ffffffffffffull;	// clear G6-G12
      if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
	x = x & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
      }
      if ((x & MASK_SNAN) == MASK_SNAN) {	// x is SNAN
	// set invalid flag
	*pfpsf |= INVALID_EXCEPTION;
	// return quiet (x)
	res = x & 0xfdffffffffffffffull;
      } else {	// x is QNaN
	// return x
	res = x;	// clear out G[6]-G[16]
      }
      BID_RETURN (res)
    }

    if (!valid_x) {
      // x is Inf. or 0

      // x is Infinity?
      if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
	// check if y is 0
	if (!coefficient_y) {
	  // y==0, return NaN
#ifdef SET_STATUS_FLAGS
	  if ((z & 0x7e00000000000000ull) != 0x7c00000000000000ull)
	    __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	  BID_RETURN (0x7c00000000000000ull);
	}
	// test if z is Inf of oposite sign
	if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)
	    && (((x ^ y) ^ z) & 0x8000000000000000ull)) {
	  // return NaN 
#ifdef SET_STATUS_FLAGS
	  __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	  BID_RETURN (0x7c00000000000000ull);
	}
	// otherwise return +/-Inf
	BID_RETURN (((x ^ y) & 0x8000000000000000ull) |
		    0x7800000000000000ull);
      }
      // x is 0
      if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)
	  && ((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {

	if (coefficient_z) {
	  exponent_y = exponent_x - DECIMAL_EXPONENT_BIAS + exponent_y;

	  sign_z = z & 0x8000000000000000ull;

	  if (exponent_y >= exponent_z)
	    BID_RETURN (z);
	  res =
	    add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,
			&rnd_mode, pfpsf);
	  BID_RETURN (res);
	}
      }
    }
    if (!valid_y) {
      // y is Inf. or 0

      // y is Infinity?
      if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
	// check if x is 0
	if (!coefficient_x) {
	  // y==0, return NaN
#ifdef SET_STATUS_FLAGS
	  __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	  BID_RETURN (0x7c00000000000000ull);
	}
	// test if z is Inf of oposite sign
	if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)
	    && (((x ^ y) ^ z) & 0x8000000000000000ull)) {
#ifdef SET_STATUS_FLAGS
	  __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	  // return NaN
	  BID_RETURN (0x7c00000000000000ull);
	}
	// otherwise return +/-Inf
	BID_RETURN (((x ^ y) & 0x8000000000000000ull) |
		    0x7800000000000000ull);
      }
      // y is 0 
      if (((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {

	if (coefficient_z) {
	  exponent_y += exponent_x - DECIMAL_EXPONENT_BIAS;

	  sign_z = z & 0x8000000000000000ull;

	  if (exponent_y >= exponent_z)
	    BID_RETURN (z);
	  res =
	    add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,
			&rnd_mode, pfpsf);
	  BID_RETURN (res);
	}
      }
    }

    if (!valid_z) {
      // y is Inf. or 0

      // test if y is NaN/Inf
      if ((z & 0x7800000000000000ull) == 0x7800000000000000ull) {
	BID_RETURN (coefficient_z & QUIET_MASK64);
      }
      // z is 0, return x*y
      if ((!coefficient_x) || (!coefficient_y)) {
	//0+/-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;
	if (exponent_x <= exponent_z)
	  res = ((UINT64) exponent_x) << 53;
	else
	  res = ((UINT64) exponent_z) << 53;
	if ((sign_x ^ sign_y) == sign_z)
	  res |= sign_z;
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
	else if (rnd_mode == ROUNDING_DOWN)
	  res |= 0x8000000000000000ull;
#endif
#endif
	BID_RETURN (res);
      }
    }
  }

  /* get binary coefficients of x and y */

  //--- 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;
    final_exponent = exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS;
    if ((final_exponent > 0) || (!coefficient_z)) {
      res =
	get_add64 (sign_x ^ sign_y,
		   final_exponent, C64, sign_z, exponent_z, coefficient_z, rnd_mode, pfpsf);
      BID_RETURN (res);
    } else {
      P.w[0] = C64;
      P.w[1] = 0;
      extra_digits = 0;
    }
  } else {
    if (!coefficient_z) {
#if DECIMAL_CALL_BY_REFERENCE
      bid64_mul (&res, px,
		 py _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		 _EXC_INFO_ARG);
#else
      res =
	bid64_mul (x,
		   y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		   _EXC_INFO_ARG);
#endif
      BID_RETURN (res);
    }
    // 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;
  }

  if (((unsigned) final_exponent) >= 3 * 256) {
    if (final_exponent < 0) {
      //--- get number of bits in the coefficients of z  ---
      tempx.d = (double) coefficient_z;
      bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
      // get number of decimal digits in the coeff_x
      digits_z = estimate_decimal_digits[bin_expon_cx];
      if (coefficient_z >= power10_table_128[digits_z].w[0])
	digits_z++;
      // underflow
      if ((final_exponent + 16 < 0)
	  || (exponent_z + digits_z > 33 + final_exponent)) {
	res =
	  BID_normalize (sign_z, exponent_z, coefficient_z,
			 sign_x ^ sign_y, 1, rnd_mode, pfpsf);
	BID_RETURN (res);
      }

      ez = exponent_z + digits_z - 16;
      if (ez < 0)
	ez = 0;
      scale_z = exponent_z - ez;
      coefficient_z *= power10_table_128[scale_z].w[0];
      ey = final_exponent - extra_digits;
      extra_digits = ez - ey;
      if (extra_digits > 33) {
	res =
	  BID_normalize (sign_z, exponent_z, coefficient_z,
			 sign_x ^ sign_y, 1, rnd_mode, pfpsf);
	BID_RETURN (res);
      }
      //else  // extra_digits<=32

      if (extra_digits > 17) {
	CYh = __truncate (P, 16);
	// get remainder
	T = power10_table_128[16].w[0];
	__mul_64x64_to_64 (CY0L, CYh, T);
	remainder_y = P.w[0] - CY0L;

	extra_digits -= 16;
	P.w[0] = CYh;
	P.w[1] = 0;
      } else
	remainder_y = 0;

      // align coeff_x, CYh
      __mul_64x64_to_128 (CZ, coefficient_z,
			  power10_table_128[extra_digits].w[0]);

      if (sign_z == (sign_y ^ sign_x)) {
	__add_128_128 (CT, CZ, P);
	if (__unsigned_compare_ge_128
	    (CT, power10_table_128[16 + extra_digits])) {
	  extra_digits++;
	  ez++;
	}
      } else {
	if (remainder_y && (__unsigned_compare_ge_128 (CZ, P))) {
	  P.w[0]++;
	  if (!P.w[0])
	    P.w[1]++;
	}
	__sub_128_128 (CT, CZ, P);
	if (((SINT64) CT.w[1]) < 0) {
	  sign_z = sign_y ^ sign_x;
	  CT.w[0] = 0 - CT.w[0];
	  CT.w[1] = 0 - CT.w[1];
	  if (CT.w[0])
	    CT.w[1]--;
	} else if(!(CT.w[1]|CT.w[0]))
		sign_z = (rnd_mode!=ROUNDING_DOWN)? 0: 0x8000000000000000ull;
	if (ez
	    &&
	    (__unsigned_compare_gt_128
	     (power10_table_128[15 + extra_digits], CT))) {
	  extra_digits--;
	  ez--;
	}
      }

#ifdef SET_STATUS_FLAGS
      uf_status = 0;
      if ((!ez)
	  &&
	  __unsigned_compare_gt_128 (power10_table_128
				     [extra_digits + 15], CT)) {
	rmode = rnd_mode;
	if (sign_z && (unsigned) (rmode - 1) < 2)
	  rmode = 3 - rmode;
	//__add_128_64(PU, CT, round_const_table[rmode][extra_digits]);
	PU = power10_table_128[extra_digits + 15];
	PU.w[0]--;
	if (__unsigned_compare_gt_128 (PU, CT)
	    || (rmode == ROUNDING_DOWN)
	    || (rmode == ROUNDING_TO_ZERO))
	  uf_status = UNDERFLOW_EXCEPTION;
	else if (extra_digits < 2) {
	  if ((rmode == ROUNDING_UP)) {
	    if (!extra_digits)
	      uf_status = UNDERFLOW_EXCEPTION;
	    else {
	      if (remainder_y && (sign_z != (sign_y ^ sign_x)))
		remainder_y = power10_table_128[16].w[0] - remainder_y;

	      if (power10_table_128[15].w[0] > remainder_y)
		uf_status = UNDERFLOW_EXCEPTION;
	    }
	  } else	// RN or RN_away
	  {
	    if (remainder_y && (sign_z != (sign_y ^ sign_x)))
	      remainder_y = power10_table_128[16].w[0] - remainder_y;

	    if (!extra_digits) {
	      remainder_y += round_const_table[rmode][15];
	      if (remainder_y < power10_table_128[16].w[0])
		uf_status = UNDERFLOW_EXCEPTION;
	    } else {
	      if (remainder_y < round_const_table[rmode][16])
		uf_status = UNDERFLOW_EXCEPTION;
	    }
	  }
	  //__set_status_flags (pfpsf, uf_status);
	}
      }
#endif
      res =
	__bid_full_round64_remainder (sign_z, ez - extra_digits, CT,
				      extra_digits, remainder_y,
				      rnd_mode, pfpsf, uf_status);
      BID_RETURN (res);

    } else {
      if ((sign_z == (sign_x ^ sign_y))
	  || (final_exponent > 3 * 256 + 15)) {
	res =
	  fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
				   1000000000000000ull, rnd_mode,
				   pfpsf);
	BID_RETURN (res);
      }
    }
  }


  if (extra_digits > 0) {
    res =
      get_add128 (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y,
		  final_exponent, P, extra_digits, rnd_mode, pfpsf);
    BID_RETURN (res);
  }
  // go to convert_format and exit
  else {
    C64 = __low_64 (P);

    res =
      get_add64 (sign_x ^ sign_y,
		 exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64, 
		 sign_z, exponent_z, coefficient_z, 
		 rnd_mode, pfpsf);
    BID_RETURN (res);
  }
}
Commit	Line	Data
7adcbafe	1	/* Copyright (C) 2007-2022 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 fma
	26	*****************************************************************************
	27	*
	28	* Algorithm description:
	29	*
	30	* if multiplication is guranteed exact (short coefficients)
b2a00c89	31	* call the unpacked arg. equivalent of bid64_add(x*y, z)
200359e8 L	32	* else
	33	* get full coefficient_x*coefficient_y product
	34	* call subroutine to perform addition of 64-bit argument
	35	* to 128-bit product
	36	*
	37	****************************************************************************/
	38
b2a00c89	39	#include "bid_inline_add.h"
200359e8 L	40
200359e8 L	41	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	42	extern void bid64_mul (UINT64 * pres, UINT64 * px,
200359e8 L	43	UINT64 *
	44	py _RND_MODE_PARAM _EXC_FLAGS_PARAM
	45	_EXC_MASKS_PARAM _EXC_INFO_PARAM);
	46	#else
	47
b2a00c89 L	48	extern UINT64 bid64_mul (UINT64 x,
	49	UINT64 y _RND_MODE_PARAM
	50	_EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	51	_EXC_INFO_PARAM);
200359e8 L	52	#endif
	53
	54	#if DECIMAL_CALL_BY_REFERENCE
	55
	56	void
b2a00c89	57	bid64_fma (UINT64 * pres, UINT64 * px, UINT64 * py,
200359e8 L	58	UINT64 *
	59	pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	60	_EXC_INFO_PARAM) {
	61	UINT64 x, y, z;
	62	#else
	63
	64	UINT64
b2a00c89	65	bid64_fma (UINT64 x, UINT64 y,
200359e8 L	66	UINT64 z _RND_MODE_PARAM _EXC_FLAGS_PARAM
	67	_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
	68	#endif
	69	UINT128 P, PU, CT, CZ;
	70	UINT64 sign_x, sign_y, coefficient_x, coefficient_y, sign_z,
	71	coefficient_z;
	72	UINT64 C64, remainder_y, res;
b2a00c89	73	UINT64 CYh, CY0L, T, valid_x, valid_y, valid_z;
200359e8	74	int_double tempx, tempy;
b2a00c89	75	int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
200359e8 L	76	bin_expon_product, rmode;
	77	int digits_p, bp, final_exponent, exponent_z, digits_z, ez, ey,
	78	scale_z, uf_status;
	79
	80	#if DECIMAL_CALL_BY_REFERENCE
	81	#if !DECIMAL_GLOBAL_ROUNDING
	82	_IDEC_round rnd_mode = *prnd_mode;
	83	#endif
	84	x = *px;
	85	y = *py;
	86	z = *pz;
	87	#endif
	88
b2a00c89 L	89	valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
	90	valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
	91	valid_z = unpack_BID64 (&sign_z, &exponent_z, &coefficient_z, z);
200359e8	92
b2a00c89 L	93	// unpack arguments, check for NaN, Infinity, or 0
	94	if (!valid_x \|\| !valid_y \|\| !valid_z) {
	95
	96	if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
	97	// if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
	98	// check first for non-canonical NaN payload
	99	y = y & 0xfe03ffffffffffffull; // clear G6-G12
	100	if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
	101	y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
200359e8	102	}
b2a00c89 L	103	if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
	104	// set invalid flag
	105	*pfpsf \|= INVALID_EXCEPTION;
	106	// return quiet (y)
	107	res = y & 0xfdffffffffffffffull;
	108	} else { // y is QNaN
	109	// return y
	110	res = y;
	111	// if z = SNaN or x = SNaN signal invalid exception
	112	if ((z & MASK_SNAN) == MASK_SNAN
	113	\|\| (x & MASK_SNAN) == MASK_SNAN) {
	114	// set invalid flag
	115	*pfpsf \|= INVALID_EXCEPTION;
	116	}
200359e8	117	}
b2a00c89 L	118	BID_RETURN (res)
	119	} else if ((z & MASK_NAN) == MASK_NAN) { // z is NAN
	120	// if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
	121	// check first for non-canonical NaN payload
	122	z = z & 0xfe03ffffffffffffull; // clear G6-G12
	123	if ((z & 0x0003ffffffffffffull) > 999999999999999ull) {
	124	z = z & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
200359e8	125	}
b2a00c89 L	126	if ((z & MASK_SNAN) == MASK_SNAN) { // z is SNAN
	127	// set invalid flag
	128	*pfpsf \|= INVALID_EXCEPTION;
	129	// return quiet (z)
	130	res = z & 0xfdffffffffffffffull;
	131	} else { // z is QNaN
	132	// return z
	133	res = z;
	134	// if x = SNaN signal invalid exception
	135	if ((x & MASK_SNAN) == MASK_SNAN) {
	136	// set invalid flag
	137	*pfpsf \|= INVALID_EXCEPTION;
	138	}
200359e8	139	}
b2a00c89 L	140	BID_RETURN (res)
	141	} else if ((x & MASK_NAN) == MASK_NAN) { // x is NAN
	142	// if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
	143	// check first for non-canonical NaN payload
	144	x = x & 0xfe03ffffffffffffull; // clear G6-G12
	145	if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
	146	x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
200359e8	147	}
b2a00c89 L	148	if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN
	149	// set invalid flag
	150	*pfpsf \|= INVALID_EXCEPTION;
	151	// return quiet (x)
	152	res = x & 0xfdffffffffffffffull;
	153	} else { // x is QNaN
	154	// return x
	155	res = x; // clear out G[6]-G[16]
200359e8	156	}
b2a00c89	157	BID_RETURN (res)
200359e8	158	}
200359e8	159
b2a00c89 L	160	if (!valid_x) {
	161	// x is Inf. or 0
	162
	163	// x is Infinity?
	164	if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
	165	// check if y is 0
	166	if (!coefficient_y) {
	167	// y==0, return NaN
200359e8	168	#ifdef SET_STATUS_FLAGS
b2a00c89 L	169	if ((z & 0x7e00000000000000ull) != 0x7c00000000000000ull)
b2a00c89 L	170	__set_status_flags (pfpsf, INVALID_EXCEPTION);
200359e8	171	#endif
b2a00c89 L	172	BID_RETURN (0x7c00000000000000ull);
	173	}
	174	// test if z is Inf of oposite sign
	175	if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)
	176	&& (((x ^ y) ^ z) & 0x8000000000000000ull)) {
	177	// return NaN
200359e8	178	#ifdef SET_STATUS_FLAGS
200359e8 L	179	__set_status_flags (pfpsf, INVALID_EXCEPTION);
200359e8 L	180	#endif
b2a00c89 L	181	BID_RETURN (0x7c00000000000000ull);
	182	}
	183	// otherwise return +/-Inf
	184	BID_RETURN (((x ^ y) & 0x8000000000000000ull) \|
	185	0x7800000000000000ull);
200359e8	186	}
b2a00c89 L	187	// x is 0
	188	if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)
	189	&& ((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {
	190
	191	if (coefficient_z) {
	192	exponent_y = exponent_x - DECIMAL_EXPONENT_BIAS + exponent_y;
	193
	194	sign_z = z & 0x8000000000000000ull;
	195
	196	if (exponent_y >= exponent_z)
	197	BID_RETURN (z);
	198	res =
	199	add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,
	200	&rnd_mode, pfpsf);
	201	BID_RETURN (res);
	202	}
	203	}
	204	}
	205	if (!valid_y) {
	206	// y is Inf. or 0
	207
	208	// y is Infinity?
	209	if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
	210	// check if x is 0
	211	if (!coefficient_x) {
	212	// y==0, return NaN
200359e8	213	#ifdef SET_STATUS_FLAGS
b2a00c89	214	__set_status_flags (pfpsf, INVALID_EXCEPTION);
200359e8	215	#endif
b2a00c89 L	216	BID_RETURN (0x7c00000000000000ull);
	217	}
	218	// test if z is Inf of oposite sign
	219	if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)
	220	&& (((x ^ y) ^ z) & 0x8000000000000000ull)) {
200359e8	221	#ifdef SET_STATUS_FLAGS
200359e8 L	222	__set_status_flags (pfpsf, INVALID_EXCEPTION);
200359e8 L	223	#endif
b2a00c89 L	224	// return NaN
	225	BID_RETURN (0x7c00000000000000ull);
	226	}
	227	// otherwise return +/-Inf
	228	BID_RETURN (((x ^ y) & 0x8000000000000000ull) \|
	229	0x7800000000000000ull);
200359e8	230	}
b2a00c89 L	231	// y is 0
b2a00c89 L	232	if (((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {
200359e8	233
b2a00c89 L	234	if (coefficient_z) {
b2a00c89 L	235	exponent_y += exponent_x - DECIMAL_EXPONENT_BIAS;
200359e8	236
b2a00c89	237	sign_z = z & 0x8000000000000000ull;
200359e8	238
b2a00c89 L	239	if (exponent_y >= exponent_z)
	240	BID_RETURN (z);
	241	res =
	242	add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,
	243	&rnd_mode, pfpsf);
	244	BID_RETURN (res);
	245	}
200359e8 L	246	}
200359e8 L	247	}
200359e8	248
b2a00c89 L	249	if (!valid_z) {
b2a00c89 L	250	// y is Inf. or 0
200359e8	251
b2a00c89 L	252	// test if y is NaN/Inf
	253	if ((z & 0x7800000000000000ull) == 0x7800000000000000ull) {
	254	BID_RETURN (coefficient_z & QUIET_MASK64);
	255	}
	256	// z is 0, return x*y
	257	if ((!coefficient_x) \|\| (!coefficient_y)) {
	258	//0+/-0
	259	exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;
	260	if (exponent_x > DECIMAL_MAX_EXPON_64)
	261	exponent_x = DECIMAL_MAX_EXPON_64;
	262	else if (exponent_x < 0)
	263	exponent_x = 0;
	264	if (exponent_x <= exponent_z)
	265	res = ((UINT64) exponent_x) << 53;
	266	else
	267	res = ((UINT64) exponent_z) << 53;
	268	if ((sign_x ^ sign_y) == sign_z)
	269	res \|= sign_z;
200359e8 L	270	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
200359e8 L	271	#ifndef IEEE_ROUND_NEAREST
b2a00c89 L	272	else if (rnd_mode == ROUNDING_DOWN)
b2a00c89 L	273	res \|= 0x8000000000000000ull;
200359e8 L	274	#endif
200359e8 L	275	#endif
b2a00c89 L	276	BID_RETURN (res);
b2a00c89 L	277	}
200359e8 L	278	}
	279	}
	280
200359e8 L	281	/* get binary coefficients of x and y */
	282
	283	//--- get number of bits in the coefficients of x and y ---
	284	// version 2 (original)
	285	tempx.d = (double) coefficient_x;
	286	bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52);
	287
	288	tempy.d = (double) coefficient_y;
	289	bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52);
	290
	291	// magnitude estimate for coefficient_x*coefficient_y is
	292	// 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
	293	bin_expon_product = bin_expon_cx + bin_expon_cy;
	294
	295	// check if coefficient_xcoefficient_y<2^(10k+3)
	296	// equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
	297	if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {
	298	// easy multiply
	299	C64 = coefficient_x * coefficient_y;
	300	final_exponent = exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS;
	301	if ((final_exponent > 0) \|\| (!coefficient_z)) {
	302	res =
b2a00c89 L	303	get_add64 (sign_x ^ sign_y,
b2a00c89 L	304	final_exponent, C64, sign_z, exponent_z, coefficient_z, rnd_mode, pfpsf);
200359e8 L	305	BID_RETURN (res);
	306	} else {
	307	P.w[0] = C64;
	308	P.w[1] = 0;
	309	extra_digits = 0;
	310	}
	311	} else {
	312	if (!coefficient_z) {
	313	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	314	bid64_mul (&res, px,
200359e8 L	315	py _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	316	_EXC_INFO_ARG);
	317	#else
	318	res =
b2a00c89	319	bid64_mul (x,
200359e8 L	320	y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	321	_EXC_INFO_ARG);
	322	#endif
	323	BID_RETURN (res);
	324	}
	325	// get 128-bit product: coefficient_x*coefficient_y
	326	__mul_64x64_to_128 (P, coefficient_x, coefficient_y);
	327
	328	// tighten binary range of P: leading bit is 2^bp
	329	// unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
	330	bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;
	331	__tight_bin_range_128 (bp, P, bin_expon_product);
	332
	333	// get number of decimal digits in the product
b2a00c89 L	334	digits_p = estimate_decimal_digits[bp];
	335	if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P)))
	336	digits_p++; // if power10_table_128[digits_p] <= P
200359e8 L	337
	338	// determine number of decimal digits to be rounded out
	339	extra_digits = digits_p - MAX_FORMAT_DIGITS;
	340	final_exponent =
	341	exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;
	342	}
	343
	344	if (((unsigned) final_exponent) >= 3 * 256) {
	345	if (final_exponent < 0) {
	346	//--- get number of bits in the coefficients of z ---
	347	tempx.d = (double) coefficient_z;
	348	bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
	349	// get number of decimal digits in the coeff_x
b2a00c89 L	350	digits_z = estimate_decimal_digits[bin_expon_cx];
b2a00c89 L	351	if (coefficient_z >= power10_table_128[digits_z].w[0])
200359e8 L	352	digits_z++;
	353	// underflow
	354	if ((final_exponent + 16 < 0)
	355	\|\| (exponent_z + digits_z > 33 + final_exponent)) {
b2a00c89 L	356	res =
	357	BID_normalize (sign_z, exponent_z, coefficient_z,
	358	sign_x ^ sign_y, 1, rnd_mode, pfpsf);
200359e8 L	359	BID_RETURN (res);
	360	}
	361
	362	ez = exponent_z + digits_z - 16;
	363	if (ez < 0)
	364	ez = 0;
	365	scale_z = exponent_z - ez;
b2a00c89	366	coefficient_z *= power10_table_128[scale_z].w[0];
200359e8 L	367	ey = final_exponent - extra_digits;
	368	extra_digits = ez - ey;
	369	if (extra_digits > 33) {
b2a00c89 L	370	res =
b2a00c89 L	371	BID_normalize (sign_z, exponent_z, coefficient_z,
200359e8	372	sign_x ^ sign_y, 1, rnd_mode, pfpsf);
200359e8 L	373	BID_RETURN (res);
	374	}
	375	//else // extra_digits<=32
	376
	377	if (extra_digits > 17) {
	378	CYh = __truncate (P, 16);
	379	// get remainder
b2a00c89	380	T = power10_table_128[16].w[0];
200359e8 L	381	__mul_64x64_to_64 (CY0L, CYh, T);
	382	remainder_y = P.w[0] - CY0L;
	383
	384	extra_digits -= 16;
	385	P.w[0] = CYh;
	386	P.w[1] = 0;
	387	} else
	388	remainder_y = 0;
	389
	390	// align coeff_x, CYh
	391	__mul_64x64_to_128 (CZ, coefficient_z,
b2a00c89	392	power10_table_128[extra_digits].w[0]);
200359e8 L	393
	394	if (sign_z == (sign_y ^ sign_x)) {
	395	__add_128_128 (CT, CZ, P);
	396	if (__unsigned_compare_ge_128
b2a00c89	397	(CT, power10_table_128[16 + extra_digits])) {
200359e8 L	398	extra_digits++;
	399	ez++;
	400	}
	401	} else {
	402	if (remainder_y && (__unsigned_compare_ge_128 (CZ, P))) {
	403	P.w[0]++;
	404	if (!P.w[0])
	405	P.w[1]++;
	406	}
	407	__sub_128_128 (CT, CZ, P);
	408	if (((SINT64) CT.w[1]) < 0) {
	409	sign_z = sign_y ^ sign_x;
	410	CT.w[0] = 0 - CT.w[0];
	411	CT.w[1] = 0 - CT.w[1];
	412	if (CT.w[0])
	413	CT.w[1]--;
b2a00c89 L	414	} else if(!(CT.w[1]\|CT.w[0]))
b2a00c89 L	415	sign_z = (rnd_mode!=ROUNDING_DOWN)? 0: 0x8000000000000000ull;
200359e8 L	416	if (ez
	417	&&
	418	(__unsigned_compare_gt_128
b2a00c89	419	(power10_table_128[15 + extra_digits], CT))) {
200359e8 L	420	extra_digits--;
	421	ez--;
	422	}
	423	}
	424
	425	#ifdef SET_STATUS_FLAGS
	426	uf_status = 0;
	427	if ((!ez)
	428	&&
b2a00c89	429	__unsigned_compare_gt_128 (power10_table_128
200359e8 L	430	[extra_digits + 15], CT)) {
	431	rmode = rnd_mode;
	432	if (sign_z && (unsigned) (rmode - 1) < 2)
	433	rmode = 3 - rmode;
b2a00c89 L	434	//__add_128_64(PU, CT, round_const_table[rmode][extra_digits]);
b2a00c89 L	435	PU = power10_table_128[extra_digits + 15];
200359e8 L	436	PU.w[0]--;
	437	if (__unsigned_compare_gt_128 (PU, CT)
	438	\|\| (rmode == ROUNDING_DOWN)
	439	\|\| (rmode == ROUNDING_TO_ZERO))
	440	uf_status = UNDERFLOW_EXCEPTION;
	441	else if (extra_digits < 2) {
	442	if ((rmode == ROUNDING_UP)) {
	443	if (!extra_digits)
	444	uf_status = UNDERFLOW_EXCEPTION;
	445	else {
	446	if (remainder_y && (sign_z != (sign_y ^ sign_x)))
b2a00c89	447	remainder_y = power10_table_128[16].w[0] - remainder_y;
200359e8	448
b2a00c89	449	if (power10_table_128[15].w[0] > remainder_y)
200359e8 L	450	uf_status = UNDERFLOW_EXCEPTION;
	451	}
	452	} else // RN or RN_away
	453	{
	454	if (remainder_y && (sign_z != (sign_y ^ sign_x)))
b2a00c89	455	remainder_y = power10_table_128[16].w[0] - remainder_y;
200359e8 L	456
200359e8 L	457	if (!extra_digits) {
b2a00c89 L	458	remainder_y += round_const_table[rmode][15];
b2a00c89 L	459	if (remainder_y < power10_table_128[16].w[0])
200359e8 L	460	uf_status = UNDERFLOW_EXCEPTION;
200359e8 L	461	} else {
b2a00c89	462	if (remainder_y < round_const_table[rmode][16])
200359e8 L	463	uf_status = UNDERFLOW_EXCEPTION;
	464	}
	465	}
	466	//__set_status_flags (pfpsf, uf_status);
	467	}
	468	}
	469	#endif
	470	res =
	471	__bid_full_round64_remainder (sign_z, ez - extra_digits, CT,
	472	extra_digits, remainder_y,
	473	rnd_mode, pfpsf, uf_status);
	474	BID_RETURN (res);
	475
	476	} else {
	477	if ((sign_z == (sign_x ^ sign_y))
	478	\|\| (final_exponent > 3 * 256 + 15)) {
	479	res =
	480	fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
	481	1000000000000000ull, rnd_mode,
	482	pfpsf);
	483	BID_RETURN (res);
	484	}
	485	}
	486	}
	487
	488
	489	if (extra_digits > 0) {
	490	res =
	491	get_add128 (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y,
	492	final_exponent, P, extra_digits, rnd_mode, pfpsf);
	493	BID_RETURN (res);
	494	}
	495	// go to convert_format and exit
	496	else {
	497	C64 = __low_64 (P);
	498
	499	res =
b2a00c89 L	500	get_add64 (sign_x ^ sign_y,
	501	exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
	502	sign_z, exponent_z, coefficient_z,
200359e8 L	503	rnd_mode, pfpsf);
	504	BID_RETURN (res);
	505	}
	506	}