1/* Copyright (C) 2007-2022 Free Software Foundation, Inc. 2 3This file is part of GCC. 4 5GCC is free software; you can redistribute it and/or modify it under 6the terms of the GNU General Public License as published by the Free 7Software Foundation; either version 3, or (at your option) any later 8version. 9 10GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11WARRANTY; without even the implied warranty of MERCHANTABILITY or 12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13for more details. 14 15Under Section 7 of GPL version 3, you are granted additional 16permissions described in the GCC Runtime Library Exception, version 173.1, as published by the Free Software Foundation. 18 19You should have received a copy of the GNU General Public License and 20a copy of the GCC Runtime Library Exception along with this program; 21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22<http://www.gnu.org/licenses/>. */ 23 24#include "bid_internal.h" 25 26#define MAX_FORMAT_DIGITS 16 27#define DECIMAL_EXPONENT_BIAS 398 28#define MAX_DECIMAL_EXPONENT 767 29 30#if DECIMAL_CALL_BY_REFERENCE 31 32void 33bid64_quantize (UINT64 * pres, UINT64 * px, 34 UINT64 * 35 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 36 _EXC_INFO_PARAM) { 37 UINT64 x, y; 38#else 39 40UINT64 41bid64_quantize (UINT64 x, 42 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM 43 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 44#endif 45 UINT128 CT; 46 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, remainder_h, C64, 47 valid_x; 48 UINT64 tmp, carry, res; 49 int_float tempx; 50 int exponent_x, exponent_y, digits_x, extra_digits, amount, amount2; 51 int expon_diff, total_digits, bin_expon_cx; 52 unsigned rmode, status; 53 54#if DECIMAL_CALL_BY_REFERENCE 55#if !DECIMAL_GLOBAL_ROUNDING 56 _IDEC_round rnd_mode = *prnd_mode; 57#endif 58 x = *px; 59 y = *py; 60#endif 61 62 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); 63 // unpack arguments, check for NaN or Infinity 64 if (!unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y)) { 65 // Inf. or NaN or 0 66#ifdef SET_STATUS_FLAGS 67 if ((x & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 68 __set_status_flags (pfpsf, INVALID_EXCEPTION); 69#endif 70 71 // x=Inf, y=Inf? 72 if (((coefficient_x << 1) == 0xf000000000000000ull) 73 && ((coefficient_y << 1) == 0xf000000000000000ull)) { 74 res = coefficient_x; 75 BID_RETURN (res); 76 } 77 // Inf or NaN? 78 if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) { 79#ifdef SET_STATUS_FLAGS 80 if (((y & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN 81 || (((y & 0x7c00000000000000ull) == 0x7800000000000000ull) && //Inf 82 ((x & 0x7c00000000000000ull) < 0x7800000000000000ull))) 83 __set_status_flags (pfpsf, INVALID_EXCEPTION); 84#endif 85 if ((y & NAN_MASK64) != NAN_MASK64) 86 coefficient_y = 0; 87 if ((x & NAN_MASK64) != NAN_MASK64) { 88 res = 0x7c00000000000000ull | (coefficient_y & QUIET_MASK64); 89 if (((y & NAN_MASK64) != NAN_MASK64) && ((x & NAN_MASK64) == 0x7800000000000000ull)) 90 res = x; 91 BID_RETURN (res); 92 } 93 } 94 } 95 // unpack arguments, check for NaN or Infinity 96 if (!valid_x) { 97 // x is Inf. or NaN or 0 98 99 // Inf or NaN? 100 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) { 101#ifdef SET_STATUS_FLAGS 102 if (((x & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN 103 || ((x & 0x7c00000000000000ull) == 0x7800000000000000ull)) //Inf 104 __set_status_flags (pfpsf, INVALID_EXCEPTION); 105#endif 106 if ((x & NAN_MASK64) != NAN_MASK64) 107 coefficient_x = 0; 108 res = 0x7c00000000000000ull | (coefficient_x & QUIET_MASK64); 109 BID_RETURN (res); 110 } 111 112 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0); 113 BID_RETURN (res); 114 } 115 // get number of decimal digits in coefficient_x 116 tempx.d = (float) coefficient_x; 117 bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f; 118 digits_x = estimate_decimal_digits[bin_expon_cx]; 119 if (coefficient_x >= power10_table_128[digits_x].w[0]) 120 digits_x++; 121 122 expon_diff = exponent_x - exponent_y; 123 total_digits = digits_x + expon_diff; 124 125 // check range of scaled coefficient 126 if ((UINT32) (total_digits + 1) <= 17) { 127 if (expon_diff >= 0) { 128 coefficient_x *= power10_table_128[expon_diff].w[0]; 129 res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x); 130 BID_RETURN (res); 131 } 132 // must round off -expon_diff digits 133 extra_digits = -expon_diff; 134#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 135#ifndef IEEE_ROUND_NEAREST 136 rmode = rnd_mode; 137 if (sign_x && (unsigned) (rmode - 1) < 2) 138 rmode = 3 - rmode; 139#else 140 rmode = 0; 141#endif 142#else 143 rmode = 0; 144#endif 145 coefficient_x += round_const_table[rmode][extra_digits]; 146 147 // get P*(2^M[extra_digits])/10^extra_digits 148 __mul_64x64_to_128 (CT, coefficient_x, 149 reciprocals10_64[extra_digits]); 150 151 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128 152 amount = short_recip_scale[extra_digits]; 153 C64 = CT.w[1] >> amount; 154#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 155#ifndef IEEE_ROUND_NEAREST 156 if (rnd_mode == 0) 157#endif 158 if (C64 & 1) { 159 // check whether fractional part of initial_P/10^extra_digits 160 // is exactly .5 161 // this is the same as fractional part of 162 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero 163 164 // get remainder 165 amount2 = 64 - amount; 166 remainder_h = 0; 167 remainder_h--; 168 remainder_h >>= amount2; 169 remainder_h = remainder_h & CT.w[1]; 170 171 // test whether fractional part is 0 172 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) { 173 C64--; 174 } 175 } 176#endif 177 178#ifdef SET_STATUS_FLAGS 179 status = INEXACT_EXCEPTION; 180 // get remainder 181 remainder_h = CT.w[1] << (64 - amount); 182 switch (rmode) { 183 case ROUNDING_TO_NEAREST: 184 case ROUNDING_TIES_AWAY: 185 // test whether fractional part is 0 186 if ((remainder_h == 0x8000000000000000ull) 187 && (CT.w[0] < reciprocals10_64[extra_digits])) 188 status = EXACT_STATUS; 189 break; 190 case ROUNDING_DOWN: 191 case ROUNDING_TO_ZERO: 192 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) 193 status = EXACT_STATUS; 194 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y; 195 break; 196 default: 197 // round up 198 __add_carry_out (tmp, carry, CT.w[0], 199 reciprocals10_64[extra_digits]); 200 if ((remainder_h >> (64 - amount)) + carry >= 201 (((UINT64) 1) << amount)) 202 status = EXACT_STATUS; 203 break; 204 } 205 __set_status_flags (pfpsf, status); 206#endif 207 208 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64); 209 BID_RETURN (res); 210 } 211 212 if (total_digits < 0) { 213#ifdef SET_STATUS_FLAGS 214 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 215#endif 216 C64 = 0; 217#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 218#ifndef IEEE_ROUND_NEAREST 219 rmode = rnd_mode; 220 if (sign_x && (unsigned) (rmode - 1) < 2) 221 rmode = 3 - rmode; 222 if (rmode == ROUNDING_UP) 223 C64 = 1; 224#endif 225#endif 226 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64); 227 BID_RETURN (res); 228 } 229 // else more than 16 digits in coefficient 230#ifdef SET_STATUS_FLAGS 231 __set_status_flags (pfpsf, INVALID_EXCEPTION); 232#endif 233 res = 0x7c00000000000000ull; 234 BID_RETURN (res); 235 236} 237