1 /* Copyright (C) 2007, 2009 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 /***************************************************************************** 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 47 void 48 bid64_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 55 UINT64 56 bid64_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