1 /****************************************************************************** 2 Copyright (c) 2007-2011, Intel Corp. 3 All rights reserved. 4 5 Redistribution and use in source and binary forms, with or without 6 modification, are permitted provided that the following conditions are met: 7 8 * Redistributions of source code must retain the above copyright notice, 9 this list of conditions and the following disclaimer. 10 * Redistributions in binary form must reproduce the above copyright 11 notice, this list of conditions and the following disclaimer in the 12 documentation and/or other materials provided with the distribution. 13 * Neither the name of Intel Corporation nor the names of its contributors 14 may be used to endorse or promote products derived from this software 15 without specific prior written permission. 16 17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF 27 THE POSSIBILITY OF SUCH DAMAGE. 28 ******************************************************************************/ 29 30 #define BID_FUNCTION_SETS_BINARY_FLAGS 31 32 #define BID_128RES 33 #include "bid_div_macros.h" 34 35 36 BID128_FUNCTION_ARG2_NORND ( bid128_rem, x, y) 37 38 BID_UINT256 P256; 39 BID_UINT128 CX, CY, CX2, CQ, CR, T, CXS, P128, res; 40 BID_UINT64 sign_x, sign_y, valid_y; 41 BID_SINT64 D; 42 int_float f64, fx; 43 int exponent_x, exponent_y, diff_expon, bin_expon_cx, scale, 44 scale0; 45 46 BID_OPT_SAVE_BINARY_FLAGS() 47 48 // unpack arguments, check for NaN or Infinity 49 50 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y); 51 52 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) { 53 #ifdef BID_SET_STATUS_FLAGS 54 if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 55 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 56 #endif 57 // test if x is NaN 58 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 59 #ifdef BID_SET_STATUS_FLAGS 60 if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 61 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 62 #endif 63 res.w[1] = CX.w[1] & QUIET_MASK64; 64 res.w[0] = CX.w[0]; 65 BID_RETURN (res); 66 } 67 // x is Infinity? 68 if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 69 // check if y is Inf. 70 if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull)) 71 // return NaN 72 { 73 #ifdef BID_SET_STATUS_FLAGS 74 // set status flags 75 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 76 #endif 77 res.w[1] = 0x7c00000000000000ull; 78 res.w[0] = 0; 79 BID_RETURN (res); 80 } 81 82 } 83 // x is 0 84 if ((!CY.w[1]) && (!CY.w[0])) { 85 #ifdef BID_SET_STATUS_FLAGS 86 // set status flags 87 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 88 #endif 89 // x=y=0, return NaN 90 res.w[1] = 0x7c00000000000000ull; 91 res.w[0] = 0; 92 BID_RETURN (res); 93 } 94 if (valid_y || ((y.w[1] & NAN_MASK64) == INFINITY_MASK64)) { 95 // return 0 96 if ((exponent_x > exponent_y) 97 && ((y.w[1] & NAN_MASK64) != INFINITY_MASK64)) 98 exponent_x = exponent_y; 99 100 res.w[1] = sign_x | (((BID_UINT64) exponent_x) << 49); 101 res.w[0] = 0; 102 BID_RETURN (res); 103 } 104 } 105 if (!valid_y) { 106 // y is Inf. or NaN 107 108 // test if y is NaN 109 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) { 110 #ifdef BID_SET_STATUS_FLAGS 111 if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN 112 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 113 #endif 114 res.w[1] = CY.w[1] & QUIET_MASK64; 115 res.w[0] = CY.w[0]; 116 BID_RETURN (res); 117 } 118 // y is Infinity? 119 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) { 120 // return x 121 res.w[1] = x.w[1]; 122 res.w[0] = x.w[0]; 123 BID_RETURN (res); 124 } 125 // y is 0 126 #ifdef BID_SET_STATUS_FLAGS 127 // set status flags 128 __set_status_flags (pfpsf, BID_INVALID_EXCEPTION); 129 #endif 130 res.w[1] = 0x7c00000000000000ull; 131 res.w[0] = 0; 132 BID_RETURN (res); 133 } 134 135 diff_expon = exponent_x - exponent_y; 136 137 if (diff_expon <= 0) { 138 diff_expon = -diff_expon; 139 140 if (diff_expon > 34) { 141 // |x|<|y| in this case 142 res = x; 143 BID_RETURN (res); 144 } 145 // set exponent of y to exponent_x, scale coefficient_y 146 T = bid_power10_table_128[diff_expon]; 147 __mul_128x128_to_256 (P256, CY, T); 148 149 if (P256.w[2] || P256.w[3]) { 150 // |x|<|y| in this case 151 res = x; 152 BID_RETURN (res); 153 } 154 155 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63); 156 CX2.w[0] = CX.w[0] << 1; 157 if (__unsigned_compare_ge_128 (P256, CX2)) { 158 // |x|<|y| in this case 159 res = x; 160 BID_RETURN (res); 161 } 162 163 P128.w[0] = P256.w[0]; 164 P128.w[1] = P256.w[1]; 165 bid___div_128_by_128 (&CQ, &CR, CX, P128); 166 167 CX2.w[1] = (CR.w[1] << 1) | (CR.w[0] >> 63); 168 CX2.w[0] = CR.w[0] << 1; 169 if ((__unsigned_compare_gt_128 (CX2, P256)) 170 || (CX2.w[1] == P256.w[1] && CX2.w[0] == P256.w[0] 171 && (CQ.w[0] & 1))) { 172 __sub_128_128 (CR, P256, CR); 173 sign_x ^= 0x8000000000000000ull; 174 } 175 176 bid_get_BID128_very_fast (&res, sign_x, exponent_x, CR); 177 BID_RETURN (res); 178 } 179 // 2^64 180 f64.i = 0x5f800000; 181 182 scale0 = 38; 183 if (!CY.w[1]) 184 scale0 = 34; 185 186 while (diff_expon > 0) { 187 // get number of digits in CX and scale=38-digits 188 // fx ~ CX 189 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0]; 190 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f; 191 scale = scale0 - bid_estimate_decimal_digits[bin_expon_cx]; 192 // scale = 38-estimate_decimal_digits[bin_expon_cx]; 193 D = CX.w[1] - bid_power10_index_binexp_128[bin_expon_cx].w[1]; 194 if (D > 0 195 || (!D && CX.w[0] >= bid_power10_index_binexp_128[bin_expon_cx].w[0])) 196 scale--; 197 198 if (diff_expon >= scale) 199 diff_expon -= scale; 200 else { 201 scale = diff_expon; 202 diff_expon = 0; 203 } 204 205 T = bid_power10_table_128[scale]; 206 __mul_128x128_low (CXS, CX, T); 207 208 bid___div_128_by_128 (&CQ, &CX, CXS, CY); 209 210 // check for remainder == 0 211 if (!CX.w[1] && !CX.w[0]) { 212 bid_get_BID128_very_fast (&res, sign_x, exponent_y, CX); 213 BID_RETURN (res); 214 } 215 } 216 217 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63); 218 CX2.w[0] = CX.w[0] << 1; 219 if ((__unsigned_compare_gt_128 (CX2, CY)) 220 || (CX2.w[1] == CY.w[1] && CX2.w[0] == CY.w[0] && (CQ.w[0] & 1))) { 221 __sub_128_128 (CX, CY, CX); 222 sign_x ^= 0x8000000000000000ull; 223 } 224 225 bid_get_BID128_very_fast (&res, sign_x, exponent_y, CX); 226 BID_RETURN (res); 227 } 228