bid64_rem.c revision 1.3
1/* Copyright (C) 2007-2013 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 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 53void 54bid64_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 60UINT64 61bid64_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