1 /* Copyright (C) 2007-2024 Free Software Foundation, Inc. 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 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 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 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/>. */ 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 32 void 33 bid64_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 40 UINT64 41 bid64_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
Indexes created Tue Mar 03 05:31:39 UTC 2026