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 remainder 26 ***************************************************************************** 27 * 28 * Algorithm description: 29 * 30 * if(exponent_x < exponent_y) 31 * scale coefficient_y so exponents are aligned 32 * perform coefficient divide (64-bit integer divide), unless 33 * coefficient_y is longer than 64 bits (clearly larger 34 * than coefficient_x) 35 * else // exponent_x > exponent_y 36 * use a loop to scale coefficient_x to 18_digits, divide by 37 * coefficient_y (64-bit integer divide), calculate remainder 38 * as new_coefficient_x and repeat until final remainder is obtained 39 * (when new_exponent_x < exponent_y) 40 * 41 ****************************************************************************/ 42 43 #include "bid_internal.h" 44 45 #define MAX_FORMAT_DIGITS 16 46 #define DECIMAL_EXPONENT_BIAS 398 47 #define MASK_BINARY_EXPONENT 0x7ff0000000000000ull 48 #define BINARY_EXPONENT_BIAS 0x3ff 49 #define UPPER_EXPON_LIMIT 51 50 51 #if DECIMAL_CALL_BY_REFERENCE 52 53 void 54 bid64_rem (UINT64 * pres, UINT64 * px, 55 UINT64 * 56 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 57 UINT64 x, y; 58 #else 59 60 UINT64 61 bid64_rem (UINT64 x, 62 UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) { 63 #endif 64 UINT128 CY; 65 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, res; 66 UINT64 Q, R, R2, T, valid_y, valid_x; 67 int_float tempx; 68 int exponent_x, exponent_y, bin_expon, e_scale; 69 int digits_x, diff_expon; 70 71 #if DECIMAL_CALL_BY_REFERENCE 72 x = *px; 73 y = *py; 74 #endif 75 76 valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y); 77 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x); 78 79 // unpack arguments, check for NaN or Infinity 80 if (!valid_x) { 81 // x is Inf. or NaN or 0 82 #ifdef SET_STATUS_FLAGS 83 if ((y & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 84 __set_status_flags (pfpsf, INVALID_EXCEPTION); 85 #endif 86 87 // test if x is NaN 88 if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 89 #ifdef SET_STATUS_FLAGS 90 if (((x & SNAN_MASK64) == SNAN_MASK64)) 91 __set_status_flags (pfpsf, INVALID_EXCEPTION); 92 #endif 93 res = coefficient_x & QUIET_MASK64;; 94 BID_RETURN (res); 95 } 96 // x is Infinity? 97 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) { 98 if (((y & NAN_MASK64) != NAN_MASK64)) { 99 #ifdef SET_STATUS_FLAGS 100 __set_status_flags (pfpsf, INVALID_EXCEPTION); 101 #endif 102 // return NaN 103 res = 0x7c00000000000000ull; 104 BID_RETURN (res); 105 } 106 } 107 // x is 0 108 // return x if y != 0 109 if (((y & 0x7800000000000000ull) < 0x7800000000000000ull) && 110 coefficient_y) { 111 if ((y & 0x6000000000000000ull) == 0x6000000000000000ull) 112 exponent_y = (y >> 51) & 0x3ff; 113 else 114 exponent_y = (y >> 53) & 0x3ff; 115 116 if (exponent_y < exponent_x) 117 exponent_x = exponent_y; 118 119 x = exponent_x; 120 x <<= 53; 121 122 res = x | sign_x; 123 BID_RETURN (res); 124 } 125 126 } 127 if (!valid_y) { 128 // y is Inf. or NaN 129 130 // test if y is NaN 131 if ((y & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 132 #ifdef SET_STATUS_FLAGS 133 if (((y & SNAN_MASK64) == SNAN_MASK64)) 134 __set_status_flags (pfpsf, INVALID_EXCEPTION); 135 #endif 136 res = coefficient_y & QUIET_MASK64;; 137 BID_RETURN (res); 138 } 139 // y is Infinity? 140 if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) { 141 res = very_fast_get_BID64 (sign_x, exponent_x, coefficient_x); 142 BID_RETURN (res); 143 } 144 // y is 0, return NaN 145 { 146 #ifdef SET_STATUS_FLAGS 147 __set_status_flags (pfpsf, INVALID_EXCEPTION); 148 #endif 149 res = 0x7c00000000000000ull; 150 BID_RETURN (res); 151 } 152 } 153 154 155 diff_expon = exponent_x - exponent_y; 156 if (diff_expon <= 0) { 157 diff_expon = -diff_expon; 158 159 if (diff_expon > 16) { 160 // |x|<|y| in this case 161 res = x; 162 BID_RETURN (res); 163 } 164 // set exponent of y to exponent_x, scale coefficient_y 165 T = power10_table_128[diff_expon].w[0]; 166 __mul_64x64_to_128 (CY, coefficient_y, T); 167 168 if (CY.w[1] || CY.w[0] > (coefficient_x << 1)) { 169 res = x; 170 BID_RETURN (res); 171 } 172 173 Q = coefficient_x / CY.w[0]; 174 R = coefficient_x - Q * CY.w[0]; 175 176 R2 = R + R; 177 if (R2 > CY.w[0] || (R2 == CY.w[0] && (Q & 1))) { 178 R = CY.w[0] - R; 179 sign_x ^= 0x8000000000000000ull; 180 } 181 182 res = very_fast_get_BID64 (sign_x, exponent_x, R); 183 BID_RETURN (res); 184 } 185 186 187 while (diff_expon > 0) { 188 // get number of digits in coeff_x 189 tempx.d = (float) coefficient_x; 190 bin_expon = ((tempx.i >> 23) & 0xff) - 0x7f; 191 digits_x = estimate_decimal_digits[bin_expon]; 192 // will not use this test, dividend will have 18 or 19 digits 193 //if(coefficient_x >= power10_table_128[digits_x].w[0]) 194 // digits_x++; 195 196 e_scale = 18 - digits_x; 197 if (diff_expon >= e_scale) { 198 diff_expon -= e_scale; 199 } else { 200 e_scale = diff_expon; 201 diff_expon = 0; 202 } 203 204 // scale dividend to 18 or 19 digits 205 coefficient_x *= power10_table_128[e_scale].w[0]; 206 207 // quotient 208 Q = coefficient_x / coefficient_y; 209 // remainder 210 coefficient_x -= Q * coefficient_y; 211 212 // check for remainder == 0 213 if (!coefficient_x) { 214 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0); 215 BID_RETURN (res); 216 } 217 } 218 219 R2 = coefficient_x + coefficient_x; 220 if (R2 > coefficient_y || (R2 == coefficient_y && (Q & 1))) { 221 coefficient_x = coefficient_y - coefficient_x; 222 sign_x ^= 0x8000000000000000ull; 223 } 224 225 res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x); 226 BID_RETURN (res); 227 228 } 229