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/***************************************************************************** 25 * BID64 multiply 26 ***************************************************************************** 27 * 28 * Algorithm description: 29 * 30 * if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed 31 * below 16) 32 * return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias, 33 * coefficient_x*coefficient_y) 34 * else 35 * get long product: coefficient_x*coefficient_y 36 * determine number of digits to round off (extra_digits) 37 * rounding is performed as a 128x128-bit multiplication by 38 * 2^M[extra_digits]/10^extra_digits, followed by a shift 39 * M[extra_digits] is sufficiently large for required accuracy 40 * 41 ****************************************************************************/ 42 43#include "bid_internal.h" 44 45#if DECIMAL_CALL_BY_REFERENCE 46 47void 48bid64_mul (UINT64 * pres, UINT64 * px, 49 UINT64 * 50 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM 51 _EXC_INFO_PARAM) { 52 UINT64 x, y; 53#else 54 55UINT64 56bid64_mul (UINT64 x, 57 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM 58 _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 59#endif 60 UINT128 P, PU, C128, Q_high, Q_low, Stemp; 61 UINT64 sign_x, sign_y, coefficient_x, coefficient_y; 62 UINT64 C64, remainder_h, carry, CY, res; 63 UINT64 valid_x, valid_y; 64 int_double tempx, tempy; 65 int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy, 66 bin_expon_product; 67 int rmode, digits_p, bp, amount, amount2, final_exponent, round_up; 68 unsigned status, uf_status; 69 70#if DECIMAL_CALL_BY_REFERENCE 71#if !DECIMAL_GLOBAL_ROUNDING 72 _IDEC_round rnd_mode = *prnd_mode; 73#endif 74 x = *px; 75 y = *py; 76#endif 77 78 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); 79 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); 80 81 // unpack arguments, check for NaN or Infinity 82 if (!valid_x) { 83 84#ifdef SET_STATUS_FLAGS 85 if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 86 __set_status_flags (pfpsf, INVALID_EXCEPTION); 87#endif 88 // x is Inf. or NaN 89 90 // test if x is NaN 91 if ((x & NAN_MASK64) == NAN_MASK64) { 92#ifdef SET_STATUS_FLAGS 93 if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN 94 __set_status_flags (pfpsf, INVALID_EXCEPTION); 95#endif 96 BID_RETURN (coefficient_x & QUIET_MASK64); 97 } 98 // x is Infinity? 99 if ((x & INFINITY_MASK64) == INFINITY_MASK64) { 100 // check if y is 0 101 if (((y & INFINITY_MASK64) != INFINITY_MASK64) 102 && !coefficient_y) { 103#ifdef SET_STATUS_FLAGS 104 __set_status_flags (pfpsf, INVALID_EXCEPTION); 105#endif 106 // y==0 , return NaN 107 BID_RETURN (NAN_MASK64); 108 } 109 // check if y is NaN 110 if ((y & NAN_MASK64) == NAN_MASK64) 111 // y==NaN , return NaN 112 BID_RETURN (coefficient_y & QUIET_MASK64); 113 // otherwise return +/-Inf 114 BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64); 115 } 116 // x is 0 117 if (((y & INFINITY_MASK64) != INFINITY_MASK64)) { 118 if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64) 119 exponent_y = ((UINT32) (y >> 51)) & 0x3ff; 120 else 121 exponent_y = ((UINT32) (y >> 53)) & 0x3ff; 122 sign_y = y & 0x8000000000000000ull; 123 124 exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS; 125 if (exponent_x > DECIMAL_MAX_EXPON_64) 126 exponent_x = DECIMAL_MAX_EXPON_64; 127 else if (exponent_x < 0) 128 exponent_x = 0; 129 BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); 130 } 131 } 132 if (!valid_y) { 133 // y is Inf. or NaN 134 135 // test if y is NaN 136 if ((y & NAN_MASK64) == NAN_MASK64) { 137#ifdef SET_STATUS_FLAGS 138 if ((y & SNAN_MASK64) == SNAN_MASK64) // sNaN 139 __set_status_flags (pfpsf, INVALID_EXCEPTION); 140#endif 141 BID_RETURN (coefficient_y & QUIET_MASK64); 142 } 143 // y is Infinity? 144 if ((y & INFINITY_MASK64) == INFINITY_MASK64) { 145 // check if x is 0 146 if (!coefficient_x) { 147 __set_status_flags (pfpsf, INVALID_EXCEPTION); 148 // x==0, return NaN 149 BID_RETURN (NAN_MASK64); 150 } 151 // otherwise return +/-Inf 152 BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64); 153 } 154 // y is 0 155 exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS; 156 if (exponent_x > DECIMAL_MAX_EXPON_64) 157 exponent_x = DECIMAL_MAX_EXPON_64; 158 else if (exponent_x < 0) 159 exponent_x = 0; 160 BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53)); 161 } 162 //--- get number of bits in the coefficients of x and y --- 163 // version 2 (original) 164 tempx.d = (double) coefficient_x; 165 bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52); 166 tempy.d = (double) coefficient_y; 167 bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52); 168 169 // magnitude estimate for coefficient_x*coefficient_y is 170 // 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx) 171 bin_expon_product = bin_expon_cx + bin_expon_cy; 172 173 // check if coefficient_x*coefficient_y<2^(10*k+3) 174 // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1 175 if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) { 176 // easy multiply 177 C64 = coefficient_x * coefficient_y; 178 179 res = 180 get_BID64_small_mantissa (sign_x ^ sign_y, 181 exponent_x + exponent_y - 182 DECIMAL_EXPONENT_BIAS, C64, rnd_mode, 183 pfpsf); 184 BID_RETURN (res); 185 } else { 186 uf_status = 0; 187 // get 128-bit product: coefficient_x*coefficient_y 188 __mul_64x64_to_128 (P, coefficient_x, coefficient_y); 189 190 // tighten binary range of P: leading bit is 2^bp 191 // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1 192 bin_expon_product -= 2 * BINARY_EXPONENT_BIAS; 193 194 __tight_bin_range_128 (bp, P, bin_expon_product); 195 196 // get number of decimal digits in the product 197 digits_p = estimate_decimal_digits[bp]; 198 if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P))) 199 digits_p++; // if power10_table_128[digits_p] <= P 200 201 // determine number of decimal digits to be rounded out 202 extra_digits = digits_p - MAX_FORMAT_DIGITS; 203 final_exponent = 204 exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS; 205 206#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 207#ifndef IEEE_ROUND_NEAREST 208 rmode = rnd_mode; 209 if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2) 210 rmode = 3 - rmode; 211#else 212 rmode = 0; 213#endif 214#else 215 rmode = 0; 216#endif 217 218 round_up = 0; 219 if (((unsigned) final_exponent) >= 3 * 256) { 220 if (final_exponent < 0) { 221 // underflow 222 if (final_exponent + 16 < 0) { 223 res = sign_x ^ sign_y; 224 __set_status_flags (pfpsf, 225 UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION); 226 if (rmode == ROUNDING_UP) 227 res |= 1; 228 BID_RETURN (res); 229 } 230 231 uf_status = UNDERFLOW_EXCEPTION; 232 if (final_exponent == -1) { 233 __add_128_64 (PU, P, round_const_table[rmode][extra_digits]); 234 if (__unsigned_compare_ge_128 235 (PU, power10_table_128[extra_digits + 16])) 236 uf_status = 0; 237 } 238 extra_digits -= final_exponent; 239 final_exponent = 0; 240 241 if (extra_digits > 17) { 242 __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]); 243 244 amount = recip_scale[16]; 245 __shr_128 (P, Q_high, amount); 246 247 // get sticky bits 248 amount2 = 64 - amount; 249 remainder_h = 0; 250 remainder_h--; 251 remainder_h >>= amount2; 252 remainder_h = remainder_h & Q_high.w[0]; 253 254 extra_digits -= 16; 255 if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1] 256 || (Q_low.w[1] == 257 reciprocals10_128[16].w[1] 258 && Q_low.w[0] >= 259 reciprocals10_128[16].w[0]))) { 260 round_up = 1; 261 __set_status_flags (pfpsf, 262 UNDERFLOW_EXCEPTION | 263 INEXACT_EXCEPTION); 264 P.w[0] = (P.w[0] << 3) + (P.w[0] << 1); 265 P.w[0] |= 1; 266 extra_digits++; 267 } 268 } 269 } else { 270 res = 271 fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, 272 1000000000000000ull, rnd_mode, 273 pfpsf); 274 BID_RETURN (res); 275 } 276 } 277 278 279 if (extra_digits > 0) { 280 // will divide by 10^(digits_p - 16) 281 282 // add a constant to P, depending on rounding mode 283 // 0.5*10^(digits_p - 16) for round-to-nearest 284 __add_128_64 (P, P, round_const_table[rmode][extra_digits]); 285 286 // get P*(2^M[extra_digits])/10^extra_digits 287 __mul_128x128_full (Q_high, Q_low, P, 288 reciprocals10_128[extra_digits]); 289 290 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 291 amount = recip_scale[extra_digits]; 292 __shr_128 (C128, Q_high, amount); 293 294 C64 = __low_64 (C128); 295 296#ifndef IEEE_ROUND_NEAREST_TIES_AWAY 297#ifndef IEEE_ROUND_NEAREST 298 if (rmode == 0) //ROUNDING_TO_NEAREST 299#endif 300 if ((C64 & 1) && !round_up) { 301 // check whether fractional part of initial_P/10^extra_digits 302 // is exactly .5 303 // this is the same as fractional part of 304 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero 305 306 // get remainder 307 remainder_h = Q_high.w[0] << (64 - amount); 308 309 // test whether fractional part is 0 310 if (!remainder_h 311 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 312 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 313 && Q_low.w[0] < 314 reciprocals10_128[extra_digits].w[0]))) { 315 C64--; 316 } 317 } 318#endif 319 320#ifdef SET_STATUS_FLAGS 321 status = INEXACT_EXCEPTION | uf_status; 322 323 // get remainder 324 remainder_h = Q_high.w[0] << (64 - amount); 325 326 switch (rmode) { 327 case ROUNDING_TO_NEAREST: 328 case ROUNDING_TIES_AWAY: 329 // test whether fractional part is 0 330 if (remainder_h == 0x8000000000000000ull 331 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 332 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 333 && Q_low.w[0] < 334 reciprocals10_128[extra_digits].w[0]))) 335 status = EXACT_STATUS; 336 break; 337 case ROUNDING_DOWN: 338 case ROUNDING_TO_ZERO: 339 if (!remainder_h 340 && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] 341 || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] 342 && Q_low.w[0] < 343 reciprocals10_128[extra_digits].w[0]))) 344 status = EXACT_STATUS; 345 break; 346 default: 347 // round up 348 __add_carry_out (Stemp.w[0], CY, Q_low.w[0], 349 reciprocals10_128[extra_digits].w[0]); 350 __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], 351 reciprocals10_128[extra_digits].w[1], CY); 352 if ((remainder_h >> (64 - amount)) + carry >= 353 (((UINT64) 1) << amount)) 354 status = EXACT_STATUS; 355 } 356 357 __set_status_flags (pfpsf, status); 358#endif 359 360 // convert to BID and return 361 res = 362 fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64, 363 rmode, pfpsf); 364 BID_RETURN (res); 365 } 366 // go to convert_format and exit 367 C64 = __low_64 (P); 368 res = 369 get_BID64 (sign_x ^ sign_y, 370 exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64, 371 rmode, pfpsf); 372 BID_RETURN (res); 373 } 374} 375