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

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

#include "bid_internal.h"

/*****************************************************************************
 *  BID64 nextup
 ****************************************************************************/

#if DECIMAL_CALL_BY_REFERENCE
void
bid64_nextup (UINT64 * pres,
	      UINT64 *
	      px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  UINT64 x = *px;
#else
UINT64
bid64_nextup (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
#endif

  UINT64 res;
  UINT64 x_sign;
  UINT64 x_exp;
  BID_UI64DOUBLE tmp1;
  int x_nr_bits;
  int q1, ind;
  UINT64 C1;			// C1 represents x_signif (UINT64)

  // check for NaNs and infinities
  if ((x & MASK_NAN) == MASK_NAN) {	// check for NaN
    if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
      x = x & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
    else
      x = x & 0xfe03ffffffffffffull;	// clear G6-G12
    if ((x & MASK_SNAN) == MASK_SNAN) {	// SNaN
      // set invalid flag
      *pfpsf |= INVALID_EXCEPTION;
      // return quiet (SNaN)
      res = x & 0xfdffffffffffffffull;
    } else {	// QNaN
      res = x;
    }
    BID_RETURN (res);
  } else if ((x & MASK_INF) == MASK_INF) {	// check for Infinity
    if (!(x & 0x8000000000000000ull)) {	// x is +inf
      res = 0x7800000000000000ull;
    } else {	// x is -inf
      res = 0xf7fb86f26fc0ffffull;	// -MAXFP = -999...99 * 10^emax
    }
    BID_RETURN (res);
  }
  // unpack the argument
  x_sign = x & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
  // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
  if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
    x_exp = (x & MASK_BINARY_EXPONENT2) >> 51;	// biased
    C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
    if (C1 > 9999999999999999ull) {	// non-canonical
      x_exp = 0;
      C1 = 0;
    }
  } else {
    x_exp = (x & MASK_BINARY_EXPONENT1) >> 53;	// biased
    C1 = x & MASK_BINARY_SIG1;
  }

  // check for zeros (possibly from non-canonical values)
  if (C1 == 0x0ull) {
    // x is 0
    res = 0x0000000000000001ull;	// MINFP = 1 * 10^emin
  } else {	// x is not special and is not zero
    if (x == 0x77fb86f26fc0ffffull) {
      // x = +MAXFP = 999...99 * 10^emax
      res = 0x7800000000000000ull;	// +inf
    } else if (x == 0x8000000000000001ull) {
      // x = -MINFP = 1...99 * 10^emin
      res = 0x8000000000000000ull;	// -0
    } else {	// -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
      // can add/subtract 1 ulp to the significand

      // Note: we could check here if x >= 10^16 to speed up the case q1 =16 
      // q1 = nr. of decimal digits in x (1 <= q1 <= 54)
      //  determine first the nr. of bits in x
      if (C1 >= MASK_BINARY_OR2) {	// x >= 2^53
	// split the 64-bit value in two 32-bit halves to avoid rounding errors
	if (C1 >= 0x0000000100000000ull) {	// x >= 2^32
	  tmp1.d = (double) (C1 >> 32);	// exact conversion
	  x_nr_bits =
	    33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	} else {	// x < 2^32
	  tmp1.d = (double) C1;	// exact conversion
	  x_nr_bits =
	    1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	}
      } else {	// if x < 2^53
	tmp1.d = (double) C1;	// exact conversion
	x_nr_bits =
	  1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
      }
      q1 = nr_digits[x_nr_bits - 1].digits;
      if (q1 == 0) {
	q1 = nr_digits[x_nr_bits - 1].digits1;
	if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
	  q1++;
      }
      // if q1 < P16 then pad the significand with zeros
      if (q1 < P16) {
	if (x_exp > (UINT64) (P16 - q1)) {
	  ind = P16 - q1;	// 1 <= ind <= P16 - 1
	  // pad with P16 - q1 zeros, until exponent = emin
	  // C1 = C1 * 10^ind
	  C1 = C1 * ten2k64[ind];
	  x_exp = x_exp - ind;
	} else {	// pad with zeros until the exponent reaches emin
	  ind = x_exp;
	  C1 = C1 * ten2k64[ind];
	  x_exp = EXP_MIN;
	}
      }
      if (!x_sign) {	// x > 0
	// add 1 ulp (add 1 to the significand)
	C1++;
	if (C1 == 0x002386f26fc10000ull) {	// if  C1 = 10^16
	  C1 = 0x00038d7ea4c68000ull;	// C1 = 10^15
	  x_exp++;
	}
	// Ok, because MAXFP = 999...99 * 10^emax was caught already
      } else {	// x < 0
	// subtract 1 ulp (subtract 1 from the significand)
	C1--;
	if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) {	// if  C1 = 10^15 - 1
	  C1 = 0x002386f26fc0ffffull;	// C1 = 10^16 - 1
	  x_exp--;
	}
      }
      // assemble the result
      // if significand has 54 bits
      if (C1 & MASK_BINARY_OR2) {
	res =
	  x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
							 MASK_BINARY_SIG2);
      } else {	// significand fits in 53 bits
	res = x_sign | (x_exp << 53) | C1;
      }
    }	// end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  }	// end x is not special and is not zero
  BID_RETURN (res);
}

/*****************************************************************************
 *  BID64 nextdown
 ****************************************************************************/

#if DECIMAL_CALL_BY_REFERENCE
void
bid64_nextdown (UINT64 * pres,
		UINT64 *
		px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  UINT64 x = *px;
#else
UINT64
bid64_nextdown (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
		_EXC_INFO_PARAM) {
#endif

  UINT64 res;
  UINT64 x_sign;
  UINT64 x_exp;
  BID_UI64DOUBLE tmp1;
  int x_nr_bits;
  int q1, ind;
  UINT64 C1;			// C1 represents x_signif (UINT64)

  // check for NaNs and infinities
  if ((x & MASK_NAN) == MASK_NAN) {	// check for NaN 
    if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
      x = x & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits 
    else
      x = x & 0xfe03ffffffffffffull;	// clear G6-G12 
    if ((x & MASK_SNAN) == MASK_SNAN) {	// SNaN 
      // set invalid flag
      *pfpsf |= INVALID_EXCEPTION;
      // return quiet (SNaN)
      res = x & 0xfdffffffffffffffull;
    } else {	// QNaN 
      res = x;
    }
    BID_RETURN (res);
  } else if ((x & MASK_INF) == MASK_INF) {	// check for Infinity
    if (x & 0x8000000000000000ull) {	// x is -inf
      res = 0xf800000000000000ull;
    } else {	// x is +inf
      res = 0x77fb86f26fc0ffffull;	// +MAXFP = +999...99 * 10^emax
    }
    BID_RETURN (res);
  }
  // unpack the argument
  x_sign = x & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
  // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
  if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
    x_exp = (x & MASK_BINARY_EXPONENT2) >> 51;	// biased
    C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
    if (C1 > 9999999999999999ull) {	// non-canonical
      x_exp = 0;
      C1 = 0;
    }
  } else {
    x_exp = (x & MASK_BINARY_EXPONENT1) >> 53;	// biased
    C1 = x & MASK_BINARY_SIG1;
  }

  // check for zeros (possibly from non-canonical values)
  if (C1 == 0x0ull) {
    // x is 0
    res = 0x8000000000000001ull;	// -MINFP = -1 * 10^emin
  } else {	// x is not special and is not zero
    if (x == 0xf7fb86f26fc0ffffull) {
      // x = -MAXFP = -999...99 * 10^emax
      res = 0xf800000000000000ull;	// -inf
    } else if (x == 0x0000000000000001ull) {
      // x = +MINFP = 1...99 * 10^emin
      res = 0x0000000000000000ull;	// -0
    } else {	// -MAXFP + 1ulp <= x <= -MINFP OR MINFP + 1 ulp <= x <= MAXFP
      // can add/subtract 1 ulp to the significand

      // Note: we could check here if x >= 10^16 to speed up the case q1 =16 
      // q1 = nr. of decimal digits in x (1 <= q1 <= 16)
      //  determine first the nr. of bits in x
      if (C1 >= 0x0020000000000000ull) {	// x >= 2^53
	// split the 64-bit value in two 32-bit halves to avoid 
	// rounding errors
	if (C1 >= 0x0000000100000000ull) {	// x >= 2^32
	  tmp1.d = (double) (C1 >> 32);	// exact conversion
	  x_nr_bits =
	    33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	} else {	// x < 2^32
	  tmp1.d = (double) C1;	// exact conversion
	  x_nr_bits =
	    1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	}
      } else {	// if x < 2^53
	tmp1.d = (double) C1;	// exact conversion
	x_nr_bits =
	  1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
      }
      q1 = nr_digits[x_nr_bits - 1].digits;
      if (q1 == 0) {
	q1 = nr_digits[x_nr_bits - 1].digits1;
	if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
	  q1++;
      }
      // if q1 < P16 then pad the significand with zeros
      if (q1 < P16) {
	if (x_exp > (UINT64) (P16 - q1)) {
	  ind = P16 - q1;	// 1 <= ind <= P16 - 1
	  // pad with P16 - q1 zeros, until exponent = emin
	  // C1 = C1 * 10^ind
	  C1 = C1 * ten2k64[ind];
	  x_exp = x_exp - ind;
	} else {	// pad with zeros until the exponent reaches emin
	  ind = x_exp;
	  C1 = C1 * ten2k64[ind];
	  x_exp = EXP_MIN;
	}
      }
      if (x_sign) {	// x < 0
	// add 1 ulp (add 1 to the significand)
	C1++;
	if (C1 == 0x002386f26fc10000ull) {	// if  C1 = 10^16
	  C1 = 0x00038d7ea4c68000ull;	// C1 = 10^15
	  x_exp++;
	  // Ok, because -MAXFP = -999...99 * 10^emax was caught already
	}
      } else {	// x > 0
	// subtract 1 ulp (subtract 1 from the significand)
	C1--;
	if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) {	// if  C1 = 10^15 - 1
	  C1 = 0x002386f26fc0ffffull;	// C1 = 10^16 - 1
	  x_exp--;
	}
      }
      // assemble the result
      // if significand has 54 bits
      if (C1 & MASK_BINARY_OR2) {
	res =
	  x_sign | (x_exp << 51) | MASK_STEERING_BITS | (C1 &
							 MASK_BINARY_SIG2);
      } else {	// significand fits in 53 bits
	res = x_sign | (x_exp << 53) | C1;
      }
    }	// end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
  }	// end x is not special and is not zero
  BID_RETURN (res);
}

/*****************************************************************************
 *  BID64 nextafter
 ****************************************************************************/

#if DECIMAL_CALL_BY_REFERENCE
void
bid64_nextafter (UINT64 * pres, UINT64 * px,
		 UINT64 *
		 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  UINT64 x = *px;
  UINT64 y = *py;
#else
UINT64
bid64_nextafter (UINT64 x,
		 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
		 _EXC_INFO_PARAM) {
#endif

  UINT64 res;
  UINT64 tmp1, tmp2;
  FPSC tmp_fpsf = 0;		// dummy fpsf for calls to comparison functions
  int res1, res2;

  // check for NaNs or infinities
  if (((x & MASK_SPECIAL) == MASK_SPECIAL) ||
      ((y & MASK_SPECIAL) == MASK_SPECIAL)) {
    // x is NaN or infinity or y is NaN or infinity

    if ((x & MASK_NAN) == MASK_NAN) {	// x is NAN
      if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
	x = x & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
      else
	x = x & 0xfe03ffffffffffffull;	// clear G6-G12
      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
	if ((y & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
	  // set invalid flag
	  *pfpsf |= INVALID_EXCEPTION;
	}
	// return x
	res = x;
      }
      BID_RETURN (res);
    } else if ((y & MASK_NAN) == MASK_NAN) {	// y is NAN
      if ((y & 0x0003ffffffffffffull) > 999999999999999ull)
	y = y & 0xfe00000000000000ull;	// clear G6-G12 and the payload bits
      else
	y = y & 0xfe03ffffffffffffull;	// clear G6-G12
      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;
      }
      BID_RETURN (res);
    } else {	// at least one is infinity
      if ((x & MASK_ANY_INF) == MASK_INF) {	// x = inf
	x = x & (MASK_SIGN | MASK_INF);
      }
      if ((y & MASK_ANY_INF) == MASK_INF) {	// y = inf
	y = y & (MASK_SIGN | MASK_INF);
      }
    }
  }
  // neither x nor y is NaN

  // if not infinity, check for non-canonical values x (treated as zero)
  if ((x & MASK_ANY_INF) != MASK_INF) {	// x != inf
    // unpack x
    if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
      // if the steering bits are 11 (condition will be 0), then
      // the exponent is G[0:w+1]
      if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
	  9999999999999999ull) {
	// non-canonical
	x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
      }
    } else {	// if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) x is unch.
      ;	// canonical
    }
  }
  // no need to check for non-canonical y

  // neither x nor y is NaN
  tmp_fpsf = *pfpsf;	// save fpsf
#if DECIMAL_CALL_BY_REFERENCE
  bid64_quiet_equal (&res1, px,
		     py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid64_quiet_greater (&res2, px,
		       py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#else
  res1 =
    bid64_quiet_equal (x,
		       y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res2 =
    bid64_quiet_greater (x,
			 y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#endif
  *pfpsf = tmp_fpsf;	// restore fpsf
  if (res1) {	// x = y
    // return x with the sign of y
    res = (y & 0x8000000000000000ull) | (x & 0x7fffffffffffffffull);
  } else if (res2) {	// x > y
#if DECIMAL_CALL_BY_REFERENCE
    bid64_nextdown (&res,
		    px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#else
    res =
      bid64_nextdown (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#endif
  } else {	// x < y
#if DECIMAL_CALL_BY_REFERENCE
    bid64_nextup (&res, px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#else
    res = bid64_nextup (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
#endif
  }
  // if the operand x is finite but the result is infinite, signal
  // overflow and inexact
  if (((x & MASK_INF) != MASK_INF) && ((res & MASK_INF) == MASK_INF)) {
    // set the inexact flag
    *pfpsf |= INEXACT_EXCEPTION;
    // set the overflow flag
    *pfpsf |= OVERFLOW_EXCEPTION;
  }
  // if the result is in (-10^emin, 10^emin), and is different from the
  // operand x, signal underflow and inexact 
  tmp1 = 0x00038d7ea4c68000ull;	// +100...0[16] * 10^emin
  tmp2 = res & 0x7fffffffffffffffull;
  tmp_fpsf = *pfpsf;	// save fpsf
#if DECIMAL_CALL_BY_REFERENCE
  bid64_quiet_greater (&res1, &tmp1,
		       &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
		       _EXC_INFO_ARG);
  bid64_quiet_not_equal (&res2, &x,
			 &res _EXC_FLAGS_ARG _EXC_MASKS_ARG
			 _EXC_INFO_ARG);
#else
  res1 =
    bid64_quiet_greater (tmp1,
			 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
			 _EXC_INFO_ARG);
  res2 =
    bid64_quiet_not_equal (x,
			   res _EXC_FLAGS_ARG _EXC_MASKS_ARG
			   _EXC_INFO_ARG);
#endif
  *pfpsf = tmp_fpsf;	// restore fpsf
  if (res1 && res2) {
    // if (bid64_quiet_greater (tmp1, tmp2, &tmp_fpsf) &&
    // bid64_quiet_not_equal (x, res, &tmp_fpsf)) {
    // set the inexact flag
    *pfpsf |= INEXACT_EXCEPTION;
    // set the underflow flag
    *pfpsf |= UNDERFLOW_EXCEPTION;
  }
  BID_RETURN (res);
}
Commit	Line	Data
8d9254fc	1	/* Copyright (C) 2007-2020 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	#include "bid_internal.h"
	25
	26	/*****************************************************************************
	27	* BID64 nextup
	28	****************************************************************************/
	29
	30	#if DECIMAL_CALL_BY_REFERENCE
	31	void
b2a00c89	32	bid64_nextup (UINT64 * pres,
200359e8 L	33	UINT64 *
	34	px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
	35	UINT64 x = *px;
	36	#else
	37	UINT64
b2a00c89	38	bid64_nextup (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
200359e8 L	39	_EXC_INFO_PARAM) {
	40	#endif
	41
	42	UINT64 res;
	43	UINT64 x_sign;
	44	UINT64 x_exp;
	45	BID_UI64DOUBLE tmp1;
	46	int x_nr_bits;
	47	int q1, ind;
b2a00c89	48	UINT64 C1; // C1 represents x_signif (UINT64)
200359e8 L	49
200359e8 L	50	// check for NaNs and infinities
b2a00c89 L	51	if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
	52	if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
	53	x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
	54	else
	55	x = x & 0xfe03ffffffffffffull; // clear G6-G12
	56	if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
	57	// set invalid flag
	58	*pfpsf \|= INVALID_EXCEPTION;
	59	// return quiet (SNaN)
	60	res = x & 0xfdffffffffffffffull;
	61	} else { // QNaN
	62	res = x;
	63	}
200359e8	64	BID_RETURN (res);
b2a00c89 L	65	} else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
b2a00c89 L	66	if (!(x & 0x8000000000000000ull)) { // x is +inf
200359e8	67	res = 0x7800000000000000ull;
b2a00c89 L	68	} else { // x is -inf
b2a00c89 L	69	res = 0xf7fb86f26fc0ffffull; // -MAXFP = -999...99 * 10^emax
200359e8 L	70	}
	71	BID_RETURN (res);
	72	}
	73	// unpack the argument
b2a00c89	74	x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8 L	75	// if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
200359e8 L	76	if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
b2a00c89	77	x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8	78	C1 = (x & MASK_BINARY_SIG2) \| MASK_BINARY_OR2;
b2a00c89	79	if (C1 > 9999999999999999ull) { // non-canonical
200359e8 L	80	x_exp = 0;
	81	C1 = 0;
	82	}
	83	} else {
b2a00c89	84	x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8 L	85	C1 = x & MASK_BINARY_SIG1;
	86	}
	87
	88	// check for zeros (possibly from non-canonical values)
	89	if (C1 == 0x0ull) {
	90	// x is 0
b2a00c89 L	91	res = 0x0000000000000001ull; // MINFP = 1 * 10^emin
b2a00c89 L	92	} else { // x is not special and is not zero
200359e8 L	93	if (x == 0x77fb86f26fc0ffffull) {
200359e8 L	94	// x = +MAXFP = 999...99 * 10^emax
b2a00c89	95	res = 0x7800000000000000ull; // +inf
200359e8 L	96	} else if (x == 0x8000000000000001ull) {
200359e8 L	97	// x = -MINFP = 1...99 * 10^emin
b2a00c89 L	98	res = 0x8000000000000000ull; // -0
b2a00c89 L	99	} else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
200359e8 L	100	// can add/subtract 1 ulp to the significand
	101
	102	// Note: we could check here if x >= 10^16 to speed up the case q1 =16
	103	// q1 = nr. of decimal digits in x (1 <= q1 <= 54)
	104	// determine first the nr. of bits in x
b2a00c89	105	if (C1 >= MASK_BINARY_OR2) { // x >= 2^53
200359e8	106	// split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89 L	107	if (C1 >= 0x0000000100000000ull) { // x >= 2^32
b2a00c89 L	108	tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8 L	109	x_nr_bits =
200359e8 L	110	33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89 L	111	} else { // x < 2^32
b2a00c89 L	112	tmp1.d = (double) C1; // exact conversion
200359e8 L	113	x_nr_bits =
	114	1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	115	}
b2a00c89 L	116	} else { // if x < 2^53
b2a00c89 L	117	tmp1.d = (double) C1; // exact conversion
200359e8 L	118	x_nr_bits =
	119	1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	120	}
b2a00c89	121	q1 = nr_digits[x_nr_bits - 1].digits;
200359e8	122	if (q1 == 0) {
b2a00c89 L	123	q1 = nr_digits[x_nr_bits - 1].digits1;
b2a00c89 L	124	if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8 L	125	q1++;
	126	}
	127	// if q1 < P16 then pad the significand with zeros
	128	if (q1 < P16) {
b2a00c89 L	129	if (x_exp > (UINT64) (P16 - q1)) {
b2a00c89 L	130	ind = P16 - q1; // 1 <= ind <= P16 - 1
200359e8 L	131	// pad with P16 - q1 zeros, until exponent = emin
200359e8 L	132	// C1 = C1 * 10^ind
b2a00c89	133	C1 = C1 * ten2k64[ind];
200359e8	134	x_exp = x_exp - ind;
b2a00c89	135	} else { // pad with zeros until the exponent reaches emin
200359e8	136	ind = x_exp;
b2a00c89	137	C1 = C1 * ten2k64[ind];
200359e8 L	138	x_exp = EXP_MIN;
	139	}
	140	}
b2a00c89	141	if (!x_sign) { // x > 0
200359e8 L	142	// add 1 ulp (add 1 to the significand)
200359e8 L	143	C1++;
b2a00c89 L	144	if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
b2a00c89 L	145	C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
200359e8 L	146	x_exp++;
	147	}
	148	// Ok, because MAXFP = 999...99 * 10^emax was caught already
b2a00c89	149	} else { // x < 0
200359e8 L	150	// subtract 1 ulp (subtract 1 from the significand)
200359e8 L	151	C1--;
b2a00c89 L	152	if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
b2a00c89 L	153	C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
200359e8 L	154	x_exp--;
	155	}
	156	}
	157	// assemble the result
	158	// if significand has 54 bits
	159	if (C1 & MASK_BINARY_OR2) {
	160	res =
	161	x_sign \| (x_exp << 51) \| MASK_STEERING_BITS \| (C1 &
	162	MASK_BINARY_SIG2);
b2a00c89	163	} else { // significand fits in 53 bits
200359e8 L	164	res = x_sign \| (x_exp << 53) \| C1;
200359e8 L	165	}
b2a00c89 L	166	} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
b2a00c89 L	167	} // end x is not special and is not zero
200359e8 L	168	BID_RETURN (res);
	169	}
	170
	171	/*****************************************************************************
	172	* BID64 nextdown
	173	****************************************************************************/
	174
	175	#if DECIMAL_CALL_BY_REFERENCE
	176	void
b2a00c89	177	bid64_nextdown (UINT64 * pres,
200359e8 L	178	UINT64 *
	179	px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
	180	UINT64 x = *px;
	181	#else
	182	UINT64
b2a00c89	183	bid64_nextdown (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
200359e8 L	184	_EXC_INFO_PARAM) {
	185	#endif
	186
	187	UINT64 res;
	188	UINT64 x_sign;
	189	UINT64 x_exp;
	190	BID_UI64DOUBLE tmp1;
	191	int x_nr_bits;
	192	int q1, ind;
b2a00c89	193	UINT64 C1; // C1 represents x_signif (UINT64)
200359e8 L	194
200359e8 L	195	// check for NaNs and infinities
b2a00c89 L	196	if ((x & MASK_NAN) == MASK_NAN) { // check for NaN
	197	if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
	198	x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
	199	else
	200	x = x & 0xfe03ffffffffffffull; // clear G6-G12
	201	if ((x & MASK_SNAN) == MASK_SNAN) { // SNaN
	202	// set invalid flag
	203	*pfpsf \|= INVALID_EXCEPTION;
	204	// return quiet (SNaN)
	205	res = x & 0xfdffffffffffffffull;
	206	} else { // QNaN
	207	res = x;
	208	}
200359e8	209	BID_RETURN (res);
b2a00c89 L	210	} else if ((x & MASK_INF) == MASK_INF) { // check for Infinity
b2a00c89 L	211	if (x & 0x8000000000000000ull) { // x is -inf
200359e8	212	res = 0xf800000000000000ull;
b2a00c89 L	213	} else { // x is +inf
b2a00c89 L	214	res = 0x77fb86f26fc0ffffull; // +MAXFP = +999...99 * 10^emax
200359e8 L	215	}
	216	BID_RETURN (res);
	217	}
	218	// unpack the argument
b2a00c89	219	x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8 L	220	// if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
200359e8 L	221	if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
b2a00c89	222	x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8	223	C1 = (x & MASK_BINARY_SIG2) \| MASK_BINARY_OR2;
b2a00c89	224	if (C1 > 9999999999999999ull) { // non-canonical
200359e8 L	225	x_exp = 0;
	226	C1 = 0;
	227	}
	228	} else {
b2a00c89	229	x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8 L	230	C1 = x & MASK_BINARY_SIG1;
	231	}
	232
	233	// check for zeros (possibly from non-canonical values)
	234	if (C1 == 0x0ull) {
	235	// x is 0
b2a00c89 L	236	res = 0x8000000000000001ull; // -MINFP = -1 * 10^emin
b2a00c89 L	237	} else { // x is not special and is not zero
200359e8 L	238	if (x == 0xf7fb86f26fc0ffffull) {
200359e8 L	239	// x = -MAXFP = -999...99 * 10^emax
b2a00c89	240	res = 0xf800000000000000ull; // -inf
200359e8 L	241	} else if (x == 0x0000000000000001ull) {
200359e8 L	242	// x = +MINFP = 1...99 * 10^emin
b2a00c89 L	243	res = 0x0000000000000000ull; // -0
b2a00c89 L	244	} else { // -MAXFP + 1ulp <= x <= -MINFP OR MINFP + 1 ulp <= x <= MAXFP
200359e8 L	245	// can add/subtract 1 ulp to the significand
	246
	247	// Note: we could check here if x >= 10^16 to speed up the case q1 =16
	248	// q1 = nr. of decimal digits in x (1 <= q1 <= 16)
	249	// determine first the nr. of bits in x
b2a00c89	250	if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 L	251	// split the 64-bit value in two 32-bit halves to avoid
200359e8 L	252	// rounding errors
b2a00c89 L	253	if (C1 >= 0x0000000100000000ull) { // x >= 2^32
b2a00c89 L	254	tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8 L	255	x_nr_bits =
200359e8 L	256	33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89 L	257	} else { // x < 2^32
b2a00c89 L	258	tmp1.d = (double) C1; // exact conversion
200359e8 L	259	x_nr_bits =
	260	1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	261	}
b2a00c89 L	262	} else { // if x < 2^53
b2a00c89 L	263	tmp1.d = (double) C1; // exact conversion
200359e8 L	264	x_nr_bits =
	265	1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
	266	}
b2a00c89	267	q1 = nr_digits[x_nr_bits - 1].digits;
200359e8	268	if (q1 == 0) {
b2a00c89 L	269	q1 = nr_digits[x_nr_bits - 1].digits1;
b2a00c89 L	270	if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8 L	271	q1++;
	272	}
	273	// if q1 < P16 then pad the significand with zeros
	274	if (q1 < P16) {
b2a00c89 L	275	if (x_exp > (UINT64) (P16 - q1)) {
b2a00c89 L	276	ind = P16 - q1; // 1 <= ind <= P16 - 1
200359e8 L	277	// pad with P16 - q1 zeros, until exponent = emin
200359e8 L	278	// C1 = C1 * 10^ind
b2a00c89	279	C1 = C1 * ten2k64[ind];
200359e8	280	x_exp = x_exp - ind;
b2a00c89	281	} else { // pad with zeros until the exponent reaches emin
200359e8	282	ind = x_exp;
b2a00c89	283	C1 = C1 * ten2k64[ind];
200359e8 L	284	x_exp = EXP_MIN;
	285	}
	286	}
b2a00c89	287	if (x_sign) { // x < 0
200359e8 L	288	// add 1 ulp (add 1 to the significand)
200359e8 L	289	C1++;
b2a00c89 L	290	if (C1 == 0x002386f26fc10000ull) { // if C1 = 10^16
b2a00c89 L	291	C1 = 0x00038d7ea4c68000ull; // C1 = 10^15
200359e8 L	292	x_exp++;
	293	// Ok, because -MAXFP = -999...99 * 10^emax was caught already
	294	}
b2a00c89	295	} else { // x > 0
200359e8 L	296	// subtract 1 ulp (subtract 1 from the significand)
200359e8 L	297	C1--;
b2a00c89 L	298	if (C1 == 0x00038d7ea4c67fffull && x_exp != 0) { // if C1 = 10^15 - 1
b2a00c89 L	299	C1 = 0x002386f26fc0ffffull; // C1 = 10^16 - 1
200359e8 L	300	x_exp--;
	301	}
	302	}
	303	// assemble the result
	304	// if significand has 54 bits
	305	if (C1 & MASK_BINARY_OR2) {
	306	res =
	307	x_sign \| (x_exp << 51) \| MASK_STEERING_BITS \| (C1 &
	308	MASK_BINARY_SIG2);
b2a00c89	309	} else { // significand fits in 53 bits
200359e8 L	310	res = x_sign \| (x_exp << 53) \| C1;
200359e8 L	311	}
b2a00c89 L	312	} // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
b2a00c89 L	313	} // end x is not special and is not zero
200359e8 L	314	BID_RETURN (res);
	315	}
	316
	317	/*****************************************************************************
	318	* BID64 nextafter
	319	****************************************************************************/
	320
	321	#if DECIMAL_CALL_BY_REFERENCE
	322	void
b2a00c89	323	bid64_nextafter (UINT64 * pres, UINT64 * px,
200359e8 L	324	UINT64 *
	325	py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
	326	UINT64 x = *px;
	327	UINT64 y = *py;
	328	#else
	329	UINT64
b2a00c89	330	bid64_nextafter (UINT64 x,
200359e8 L	331	UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	332	_EXC_INFO_PARAM) {
	333	#endif
	334
	335	UINT64 res;
	336	UINT64 tmp1, tmp2;
b2a00c89	337	FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
200359e8 L	338	int res1, res2;
200359e8 L	339
b2a00c89 L	340	// check for NaNs or infinities
	341	if (((x & MASK_SPECIAL) == MASK_SPECIAL) \|\|
	342	((y & MASK_SPECIAL) == MASK_SPECIAL)) {
	343	// x is NaN or infinity or y is NaN or infinity
200359e8	344
b2a00c89 L	345	if ((x & MASK_NAN) == MASK_NAN) { // x is NAN
	346	if ((x & 0x0003ffffffffffffull) > 999999999999999ull)
	347	x = x & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
	348	else
	349	x = x & 0xfe03ffffffffffffull; // clear G6-G12
	350	if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNAN
200359e8 L	351	// set invalid flag
	352	*pfpsf \|= INVALID_EXCEPTION;
	353	// return quiet (x)
	354	res = x & 0xfdffffffffffffffull;
b2a00c89 L	355	} else { // x is QNaN
b2a00c89 L	356	if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
200359e8 L	357	// set invalid flag
	358	*pfpsf \|= INVALID_EXCEPTION;
	359	}
	360	// return x
	361	res = x;
	362	}
	363	BID_RETURN (res);
b2a00c89 L	364	} else if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
	365	if ((y & 0x0003ffffffffffffull) > 999999999999999ull)
	366	y = y & 0xfe00000000000000ull; // clear G6-G12 and the payload bits
	367	else
	368	y = y & 0xfe03ffffffffffffull; // clear G6-G12
	369	if ((y & MASK_SNAN) == MASK_SNAN) { // y is SNAN
200359e8 L	370	// set invalid flag
	371	*pfpsf \|= INVALID_EXCEPTION;
	372	// return quiet (y)
	373	res = y & 0xfdffffffffffffffull;
b2a00c89	374	} else { // y is QNaN
200359e8 L	375	// return y
	376	res = y;
	377	}
	378	BID_RETURN (res);
b2a00c89 L	379	} else { // at least one is infinity
	380	if ((x & MASK_ANY_INF) == MASK_INF) { // x = inf
	381	x = x & (MASK_SIGN \| MASK_INF);
	382	}
	383	if ((y & MASK_ANY_INF) == MASK_INF) { // y = inf
	384	y = y & (MASK_SIGN \| MASK_INF);
	385	}
200359e8 L	386	}
	387	}
	388	// neither x nor y is NaN
b2a00c89 L	389
	390	// if not infinity, check for non-canonical values x (treated as zero)
	391	if ((x & MASK_ANY_INF) != MASK_INF) { // x != inf
	392	// unpack x
	393	if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
	394	// if the steering bits are 11 (condition will be 0), then
	395	// the exponent is G[0:w+1]
	396	if (((x & MASK_BINARY_SIG2) \| MASK_BINARY_OR2) >
	397	9999999999999999ull) {
	398	// non-canonical
	399	x = (x & MASK_SIGN) \| ((x & MASK_BINARY_EXPONENT2) << 2);
	400	}
	401	} else { // if ((x & MASK_STEERING_BITS) != MASK_STEERING_BITS) x is unch.
	402	; // canonical
	403	}
	404	}
	405	// no need to check for non-canonical y
	406
	407	// neither x nor y is NaN
	408	tmp_fpsf = *pfpsf; // save fpsf
200359e8	409	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	410	bid64_quiet_equal (&res1, px,
200359e8	411	py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
b2a00c89	412	bid64_quiet_greater (&res2, px,
200359e8 L	413	py _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
	414	#else
	415	res1 =
b2a00c89	416	bid64_quiet_equal (x,
200359e8 L	417	y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
200359e8 L	418	res2 =
b2a00c89	419	bid64_quiet_greater (x,
200359e8 L	420	y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
200359e8 L	421	#endif
b2a00c89 L	422	*pfpsf = tmp_fpsf; // restore fpsf
b2a00c89 L	423	if (res1) { // x = y
200359e8	424	// return x with the sign of y
b2a00c89 L	425	res = (y & 0x8000000000000000ull) \| (x & 0x7fffffffffffffffull);
b2a00c89 L	426	} else if (res2) { // x > y
200359e8	427	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	428	bid64_nextdown (&res,
200359e8 L	429	px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
	430	#else
	431	res =
b2a00c89	432	bid64_nextdown (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
200359e8	433	#endif
b2a00c89	434	} else { // x < y
200359e8	435	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	436	bid64_nextup (&res, px _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
200359e8	437	#else
b2a00c89	438	res = bid64_nextup (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
200359e8 L	439	#endif
	440	}
	441	// if the operand x is finite but the result is infinite, signal
	442	// overflow and inexact
	443	if (((x & MASK_INF) != MASK_INF) && ((res & MASK_INF) == MASK_INF)) {
	444	// set the inexact flag
	445	*pfpsf \|= INEXACT_EXCEPTION;
	446	// set the overflow flag
	447	*pfpsf \|= OVERFLOW_EXCEPTION;
	448	}
	449	// if the result is in (-10^emin, 10^emin), and is different from the
	450	// operand x, signal underflow and inexact
b2a00c89	451	tmp1 = 0x00038d7ea4c68000ull; // +100...0[16] * 10^emin
200359e8	452	tmp2 = res & 0x7fffffffffffffffull;
b2a00c89	453	tmp_fpsf = *pfpsf; // save fpsf
200359e8	454	#if DECIMAL_CALL_BY_REFERENCE
b2a00c89	455	bid64_quiet_greater (&res1, &tmp1,
200359e8 L	456	&tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
200359e8 L	457	_EXC_INFO_ARG);
b2a00c89	458	bid64_quiet_not_equal (&res2, &x,
200359e8 L	459	&res _EXC_FLAGS_ARG _EXC_MASKS_ARG
	460	_EXC_INFO_ARG);
	461	#else
	462	res1 =
b2a00c89	463	bid64_quiet_greater (tmp1,
200359e8 L	464	tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
	465	_EXC_INFO_ARG);
	466	res2 =
b2a00c89	467	bid64_quiet_not_equal (x,
200359e8 L	468	res _EXC_FLAGS_ARG _EXC_MASKS_ARG
	469	_EXC_INFO_ARG);
	470	#endif
b2a00c89	471	*pfpsf = tmp_fpsf; // restore fpsf
200359e8	472	if (res1 && res2) {
b2a00c89 L	473	// if (bid64_quiet_greater (tmp1, tmp2, &tmp_fpsf) &&
b2a00c89 L	474	// bid64_quiet_not_equal (x, res, &tmp_fpsf)) {
200359e8 L	475	// set the inexact flag
	476	*pfpsf \|= INEXACT_EXCEPTION;
	477	// set the underflow flag
	478	*pfpsf \|= UNDERFLOW_EXCEPTION;
	479	}
	480	BID_RETURN (res);
	481	}