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

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

#define BID_128RES
#include "bid_internal.h"

BID128_FUNCTION_ARG2 (bid128_quantize, x, y)

     UINT256 CT;
     UINT128 CX, CY, T, CX2, CR, Stemp, res, REM_H, C2N;
     UINT64 sign_x, sign_y, remainder_h, carry, CY64, valid_x;
     int_float tempx;
     int exponent_x, exponent_y, digits_x, extra_digits, amount;
     int expon_diff, total_digits, bin_expon_cx, rmode, status;

valid_x = unpack_BID128_value (&sign_x, &exponent_x, &CX, x);

  // unpack arguments, check for NaN or Infinity
if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
    // y is Inf. or NaN
#ifdef SET_STATUS_FLAGS
if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64)	// y is sNaN
  __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif

    // test if y is NaN
if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
#ifdef SET_STATUS_FLAGS
  if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
    // set status flags
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
  }
#endif
  if ((x.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull) {
    res.w[1] = CY.w[1] & QUIET_MASK64;
    res.w[0] = CY.w[0];
  } else {
    res.w[1] = CX.w[1] & QUIET_MASK64;
    res.w[0] = CX.w[0];
  }
  BID_RETURN (res);
}
    // y is Infinity?
if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
  // check if x is not Inf.
  if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) {
    // return NaN 
#ifdef SET_STATUS_FLAGS
    // set status flags
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
    res.w[1] = 0x7c00000000000000ull;
    res.w[0] = 0;
    BID_RETURN (res);
  } else
    if (((x.w[1] & 0x7c00000000000000ull) <= 0x7800000000000000ull)) {
    res.w[1] = CX.w[1] & QUIET_MASK64;
    res.w[0] = CX.w[0];
    BID_RETURN (res);
  }
}

}

if (!valid_x) {
  // test if x is NaN or Inf
  if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) {
#ifdef SET_STATUS_FLAGS
    // set status flags
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
    res.w[1] = 0x7c00000000000000ull;
    res.w[0] = 0;
    BID_RETURN (res);
  } else if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
    if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
#ifdef SET_STATUS_FLAGS
      // set status flags
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
    }
    res.w[1] = CX.w[1] & QUIET_MASK64;
    res.w[0] = CX.w[0];
    BID_RETURN (res);
  }
  if (!CX.w[1] && !CX.w[0]) {
    get_BID128_very_fast (&res, sign_x, exponent_y, CX);
    BID_RETURN (res);
  }
}
  // get number of decimal digits in coefficient_x
if (CX.w[1]) {
  tempx.d = (float) CX.w[1];
  bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f + 64;
} else {
  tempx.d = (float) CX.w[0];
  bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
}

digits_x = estimate_decimal_digits[bin_expon_cx];
if (CX.w[1] > power10_table_128[digits_x].w[1]
    || (CX.w[1] == power10_table_128[digits_x].w[1]
	&& CX.w[0] >= power10_table_128[digits_x].w[0]))
  digits_x++;

expon_diff = exponent_x - exponent_y;
total_digits = digits_x + expon_diff;

if ((UINT32) total_digits <= 34) {
  if (expon_diff >= 0) {
    T = power10_table_128[expon_diff];
    __mul_128x128_low (CX2, T, CX);
    get_BID128_very_fast (&res, sign_x, exponent_y, CX2);
    BID_RETURN (res);
  }
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  rmode = rnd_mode;
  if (sign_x && (unsigned) (rmode - 1) < 2)
    rmode = 3 - rmode;
#else
  rmode = 0;
#endif
#else
  rmode = 0;
#endif
  // must round off -expon_diff digits
  extra_digits = -expon_diff;
  __add_128_128 (CX, CX, round_const_table_128[rmode][extra_digits]);

  // get P*(2^M[extra_digits])/10^extra_digits
  __mul_128x128_to_256 (CT, CX, reciprocals10_128[extra_digits]);

  // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
  amount = recip_scale[extra_digits];
  CX2.w[0] = CT.w[2];
  CX2.w[1] = CT.w[3];
  if (amount >= 64) {
    CR.w[1] = 0;
    CR.w[0] = CX2.w[1] >> (amount - 64);
  } else {
    __shr_128 (CR, CX2, amount);
  }

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  if (rnd_mode == 0)
#endif
    if (CR.w[0] & 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
      if (amount >= 64) {
	remainder_h = CX2.w[0] | (CX2.w[1] << (128 - amount));
      } else
	remainder_h = CX2.w[0] << (64 - amount);

      // test whether fractional part is 0
      if (!remainder_h
	  && (CT.w[1] < reciprocals10_128[extra_digits].w[1]
	      || (CT.w[1] == reciprocals10_128[extra_digits].w[1]
		  && CT.w[0] < reciprocals10_128[extra_digits].w[0]))) {
	CR.w[0]--;
      }
    }
#endif

#ifdef SET_STATUS_FLAGS
  status = INEXACT_EXCEPTION;

  // get remainder
  if (amount >= 64) {
    REM_H.w[1] = (CX2.w[1] << (128 - amount));
    REM_H.w[0] = CX2.w[0];
  } else {
    REM_H.w[1] = CX2.w[0] << (64 - amount);
    REM_H.w[0] = 0;
  }

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

  __set_status_flags (pfpsf, status);

#endif
  get_BID128_very_fast (&res, sign_x, exponent_y, CR);
  BID_RETURN (res);
}
if (total_digits < 0) {
  CR.w[1] = CR.w[0] = 0;
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
  rmode = rnd_mode;
  if (sign_x && (unsigned) (rmode - 1) < 2)
    rmode = 3 - rmode;
  if (rmode == ROUNDING_UP)
    CR.w[0] = 1;
#endif
#endif
#ifdef SET_STATUS_FLAGS
  __set_status_flags (pfpsf, INEXACT_EXCEPTION);
#endif
  get_BID128_very_fast (&res, sign_x, exponent_y, CR);
  BID_RETURN (res);
}
  // else  more than 34 digits in coefficient
#ifdef SET_STATUS_FLAGS
__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
res.w[1] = 0x7c00000000000000ull;
res.w[0] = 0;
BID_RETURN (res);

}
Commit	Line	Data
fbd26352	1	/* Copyright (C) 2007-2019 Free Software Foundation, Inc.
9b6b0236	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
6bc9506f	7	Software Foundation; either version 3, or (at your option) any later
9b6b0236	8	version.
9b6b0236	9
9b6b0236	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
6bc9506f	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/>. */
9b6b0236	23
	24	#define BID_128RES
	25	#include "bid_internal.h"
	26
84d1fc49	27	BID128_FUNCTION_ARG2 (bid128_quantize, x, y)
9b6b0236	28
84d1fc49	29	UINT256 CT;
	30	UINT128 CX, CY, T, CX2, CR, Stemp, res, REM_H, C2N;
	31	UINT64 sign_x, sign_y, remainder_h, carry, CY64, valid_x;
	32	int_float tempx;
	33	int exponent_x, exponent_y, digits_x, extra_digits, amount;
	34	int expon_diff, total_digits, bin_expon_cx, rmode, status;
	35
	36	valid_x = unpack_BID128_value (&sign_x, &exponent_x, &CX, x);
9b6b0236	37
9b6b0236	38	// unpack arguments, check for NaN or Infinity
84d1fc49	39	if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
9b6b0236	40	// y is Inf. or NaN
9b6b0236	41	#ifdef SET_STATUS_FLAGS
84d1fc49	42	if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
84d1fc49	43	__set_status_flags (pfpsf, INVALID_EXCEPTION);
9b6b0236	44	#endif
	45
	46	// test if y is NaN
84d1fc49	47	if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
9b6b0236	48	#ifdef SET_STATUS_FLAGS
84d1fc49	49	if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
	50	// set status flags
	51	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	52	}
9b6b0236	53	#endif
84d1fc49	54	if ((x.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull) {
	55	res.w[1] = CY.w[1] & QUIET_MASK64;
	56	res.w[0] = CY.w[0];
	57	} else {
	58	res.w[1] = CX.w[1] & QUIET_MASK64;
	59	res.w[0] = CX.w[0];
	60	}
	61	BID_RETURN (res);
	62	}
9b6b0236	63	// y is Infinity?
84d1fc49	64	if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
	65	// check if x is not Inf.
	66	if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) {
	67	// return NaN
9b6b0236	68	#ifdef SET_STATUS_FLAGS
84d1fc49	69	// set status flags
84d1fc49	70	__set_status_flags (pfpsf, INVALID_EXCEPTION);
9b6b0236	71	#endif
84d1fc49	72	res.w[1] = 0x7c00000000000000ull;
	73	res.w[0] = 0;
	74	BID_RETURN (res);
	75	} else
	76	if (((x.w[1] & 0x7c00000000000000ull) <= 0x7800000000000000ull)) {
	77	res.w[1] = CX.w[1] & QUIET_MASK64;
	78	res.w[0] = CX.w[0];
	79	BID_RETURN (res);
9b6b0236	80	}
84d1fc49	81	}
9b6b0236	82
84d1fc49	83	}
	84
	85	if (!valid_x) {
	86	// test if x is NaN or Inf
	87	if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) {
9b6b0236	88	#ifdef SET_STATUS_FLAGS
84d1fc49	89	// set status flags
84d1fc49	90	__set_status_flags (pfpsf, INVALID_EXCEPTION);
9b6b0236	91	#endif
84d1fc49	92	res.w[1] = 0x7c00000000000000ull;
	93	res.w[0] = 0;
	94	BID_RETURN (res);
	95	} else if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
	96	if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
9b6b0236	97	#ifdef SET_STATUS_FLAGS
84d1fc49	98	// set status flags
84d1fc49	99	__set_status_flags (pfpsf, INVALID_EXCEPTION);
9b6b0236	100	#endif
9b6b0236	101	}
84d1fc49	102	res.w[1] = CX.w[1] & QUIET_MASK64;
	103	res.w[0] = CX.w[0];
	104	BID_RETURN (res);
9b6b0236	105	}
84d1fc49	106	if (!CX.w[1] && !CX.w[0]) {
	107	get_BID128_very_fast (&res, sign_x, exponent_y, CX);
	108	BID_RETURN (res);
9b6b0236	109	}
84d1fc49	110	}
	111	// get number of decimal digits in coefficient_x
	112	if (CX.w[1]) {
	113	tempx.d = (float) CX.w[1];
	114	bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f + 64;
	115	} else {
	116	tempx.d = (float) CX.w[0];
	117	bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
	118	}
9b6b0236	119
84d1fc49	120	digits_x = estimate_decimal_digits[bin_expon_cx];
	121	if (CX.w[1] > power10_table_128[digits_x].w[1]
	122	\|\| (CX.w[1] == power10_table_128[digits_x].w[1]
	123	&& CX.w[0] >= power10_table_128[digits_x].w[0]))
	124	digits_x++;
9b6b0236	125
84d1fc49	126	expon_diff = exponent_x - exponent_y;
	127	total_digits = digits_x + expon_diff;
	128
	129	if ((UINT32) total_digits <= 34) {
	130	if (expon_diff >= 0) {
	131	T = power10_table_128[expon_diff];
	132	__mul_128x128_low (CX2, T, CX);
	133	get_BID128_very_fast (&res, sign_x, exponent_y, CX2);
	134	BID_RETURN (res);
	135	}
9b6b0236	136	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
9b6b0236	137	#ifndef IEEE_ROUND_NEAREST
84d1fc49	138	rmode = rnd_mode;
	139	if (sign_x && (unsigned) (rmode - 1) < 2)
	140	rmode = 3 - rmode;
9b6b0236	141	#else
84d1fc49	142	rmode = 0;
9b6b0236	143	#endif
9b6b0236	144	#else
84d1fc49	145	rmode = 0;
9b6b0236	146	#endif
84d1fc49	147	// must round off -expon_diff digits
	148	extra_digits = -expon_diff;
	149	__add_128_128 (CX, CX, round_const_table_128[rmode][extra_digits]);
9b6b0236	150
84d1fc49	151	// get P*(2^M[extra_digits])/10^extra_digits
84d1fc49	152	__mul_128x128_to_256 (CT, CX, reciprocals10_128[extra_digits]);
9b6b0236	153
84d1fc49	154	// now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
	155	amount = recip_scale[extra_digits];
	156	CX2.w[0] = CT.w[2];
	157	CX2.w[1] = CT.w[3];
	158	if (amount >= 64) {
	159	CR.w[1] = 0;
	160	CR.w[0] = CX2.w[1] >> (amount - 64);
	161	} else {
	162	__shr_128 (CR, CX2, amount);
	163	}
9b6b0236	164
	165	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
	166	#ifndef IEEE_ROUND_NEAREST
84d1fc49	167	if (rnd_mode == 0)
9b6b0236	168	#endif
84d1fc49	169	if (CR.w[0] & 1) {
	170	// check whether fractional part of initial_P/10^extra_digits is
	171	// exactly .5 this is the same as fractional part of
	172	// (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
9b6b0236	173
84d1fc49	174	// get remainder
	175	if (amount >= 64) {
	176	remainder_h = CX2.w[0] \| (CX2.w[1] << (128 - amount));
	177	} else
	178	remainder_h = CX2.w[0] << (64 - amount);
9b6b0236	179
84d1fc49	180	// test whether fractional part is 0
	181	if (!remainder_h
	182	&& (CT.w[1] < reciprocals10_128[extra_digits].w[1]
	183	\|\| (CT.w[1] == reciprocals10_128[extra_digits].w[1]
	184	&& CT.w[0] < reciprocals10_128[extra_digits].w[0]))) {
	185	CR.w[0]--;
9b6b0236	186	}
84d1fc49	187	}
9b6b0236	188	#endif
	189
	190	#ifdef SET_STATUS_FLAGS
84d1fc49	191	status = INEXACT_EXCEPTION;
9b6b0236	192
84d1fc49	193	// get remainder
	194	if (amount >= 64) {
	195	REM_H.w[1] = (CX2.w[1] << (128 - amount));
	196	REM_H.w[0] = CX2.w[0];
	197	} else {
	198	REM_H.w[1] = CX2.w[0] << (64 - amount);
	199	REM_H.w[0] = 0;
	200	}
9b6b0236	201
84d1fc49	202	switch (rmode) {
	203	case ROUNDING_TO_NEAREST:
	204	case ROUNDING_TIES_AWAY:
	205	// test whether fractional part is 0
	206	if (REM_H.w[1] == 0x8000000000000000ull && !REM_H.w[0]
	207	&& (CT.w[1] < reciprocals10_128[extra_digits].w[1]
	208	\|\| (CT.w[1] == reciprocals10_128[extra_digits].w[1]
	209	&& CT.w[0] < reciprocals10_128[extra_digits].w[0])))
	210	status = EXACT_STATUS;
	211	break;
	212	case ROUNDING_DOWN:
	213	case ROUNDING_TO_ZERO:
	214	if (!(REM_H.w[1] \| REM_H.w[0])
	215	&& (CT.w[1] < reciprocals10_128[extra_digits].w[1]
	216	\|\| (CT.w[1] == reciprocals10_128[extra_digits].w[1]
	217	&& CT.w[0] < reciprocals10_128[extra_digits].w[0])))
	218	status = EXACT_STATUS;
	219	break;
	220	default:
	221	// round up
	222	__add_carry_out (Stemp.w[0], CY64, CT.w[0],
	223	reciprocals10_128[extra_digits].w[0]);
	224	__add_carry_in_out (Stemp.w[1], carry, CT.w[1],
	225	reciprocals10_128[extra_digits].w[1], CY64);
	226	if (amount < 64) {
	227	C2N.w[1] = 0;
	228	C2N.w[0] = ((UINT64) 1) << amount;
	229	REM_H.w[0] = REM_H.w[1] >> (64 - amount);
	230	REM_H.w[1] = 0;
	231	} else {
	232	C2N.w[1] = ((UINT64) 1) << (amount - 64);
	233	C2N.w[0] = 0;
	234	REM_H.w[1] >>= (128 - amount);
9b6b0236	235	}
84d1fc49	236	REM_H.w[0] += carry;
	237	if (REM_H.w[0] < carry)
	238	REM_H.w[1]++;
	239	if (__unsigned_compare_ge_128 (REM_H, C2N))
	240	status = EXACT_STATUS;
	241	}
9b6b0236	242
84d1fc49	243	__set_status_flags (pfpsf, status);
9b6b0236	244
9b6b0236	245	#endif
84d1fc49	246	get_BID128_very_fast (&res, sign_x, exponent_y, CR);
	247	BID_RETURN (res);
	248	}
	249	if (total_digits < 0) {
	250	CR.w[1] = CR.w[0] = 0;
9b6b0236	251	#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
9b6b0236	252	#ifndef IEEE_ROUND_NEAREST
84d1fc49	253	rmode = rnd_mode;
	254	if (sign_x && (unsigned) (rmode - 1) < 2)
	255	rmode = 3 - rmode;
	256	if (rmode == ROUNDING_UP)
	257	CR.w[0] = 1;
9b6b0236	258	#endif
	259	#endif
	260	#ifdef SET_STATUS_FLAGS
84d1fc49	261	__set_status_flags (pfpsf, INEXACT_EXCEPTION);
9b6b0236	262	#endif
84d1fc49	263	get_BID128_very_fast (&res, sign_x, exponent_y, CR);
	264	BID_RETURN (res);
	265	}
9b6b0236	266	// else more than 34 digits in coefficient
9b6b0236	267	#ifdef SET_STATUS_FLAGS
84d1fc49	268	__set_status_flags (pfpsf, INVALID_EXCEPTION);
9b6b0236	269	#endif
84d1fc49	270	res.w[1] = 0x7c00000000000000ull;
	271	res.w[0] = 0;
	272	BID_RETURN (res);
9b6b0236	273
9b6b0236	274	}