1 /* Copyright (C) 2007-2021 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 #define BID_128RES 25 #include "bid_internal.h" 26 27 BID128_FUNCTION_ARG2 (bid128_quantize, x, y) 28 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); 37 38 // unpack arguments, check for NaN or Infinity 39 if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) { 40 // y is Inf. or NaN 41 #ifdef SET_STATUS_FLAGS 42 if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 43 __set_status_flags (pfpsf, INVALID_EXCEPTION); 44 #endif 45 46 // test if y is NaN 47 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 48 #ifdef SET_STATUS_FLAGS 49 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) { 50 // set status flags 51 __set_status_flags (pfpsf, INVALID_EXCEPTION); 52 } 53 #endif 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 } 63 // y is Infinity? 64 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 65 // check if x is not Inf. 66 if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) { 67 // return NaN 68 #ifdef SET_STATUS_FLAGS 69 // set status flags 70 __set_status_flags (pfpsf, INVALID_EXCEPTION); 71 #endif 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); 80 } 81 } 82 83 } 84 85 if (!valid_x) { 86 // test if x is NaN or Inf 87 if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) { 88 #ifdef SET_STATUS_FLAGS 89 // set status flags 90 __set_status_flags (pfpsf, INVALID_EXCEPTION); 91 #endif 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) { 97 #ifdef SET_STATUS_FLAGS 98 // set status flags 99 __set_status_flags (pfpsf, INVALID_EXCEPTION); 100 #endif 101 } 102 res.w[1] = CX.w[1] & QUIET_MASK64; 103 res.w[0] = CX.w[0]; 104 BID_RETURN (res); 105 } 106 if (!CX.w[1] && !CX.w[0]) { 107 get_BID128_very_fast (&res, sign_x, exponent_y, CX); 108 BID_RETURN (res); 109 } 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 } 119 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++; 125 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 } 136 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 137 #ifndef IEEE_ROUND_NEAREST 138 rmode = rnd_mode; 139 if (sign_x && (unsigned) (rmode - 1) < 2) 140 rmode = 3 - rmode; 141 #else 142 rmode = 0; 143 #endif 144 #else 145 rmode = 0; 146 #endif 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]); 150 151 // get P*(2^M[extra_digits])/10^extra_digits 152 __mul_128x128_to_256 (CT, CX, reciprocals10_128[extra_digits]); 153 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 } 164 165 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 166 #ifndef IEEE_ROUND_NEAREST 167 if (rnd_mode == 0) 168 #endif 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 173 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); 179 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]--; 186 } 187 } 188 #endif 189 190 #ifdef SET_STATUS_FLAGS 191 status = INEXACT_EXCEPTION; 192 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 } 201 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); 235 } 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 } 242 243 __set_status_flags (pfpsf, status); 244 245 #endif 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; 251 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY 252 #ifndef IEEE_ROUND_NEAREST 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; 258 #endif 259 #endif 260 #ifdef SET_STATUS_FLAGS 261 __set_status_flags (pfpsf, INEXACT_EXCEPTION); 262 #endif 263 get_BID128_very_fast (&res, sign_x, exponent_y, CR); 264 BID_RETURN (res); 265 } 266 // else more than 34 digits in coefficient 267 #ifdef SET_STATUS_FLAGS 268 __set_status_flags (pfpsf, INVALID_EXCEPTION); 269 #endif 270 res.w[1] = 0x7c00000000000000ull; 271 res.w[0] = 0; 272 BID_RETURN (res); 273 274 } 275