]>
Commit | Line | Data |
---|---|---|
a5544970 | 1 | /* Copyright (C) 2007-2019 Free Software Foundation, Inc. |
200359e8 L |
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 | |
748086b7 | 7 | Software Foundation; either version 3, or (at your option) any later |
200359e8 L |
8 | version. |
9 | ||
200359e8 L |
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 | ||
748086b7 JJ |
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/>. */ | |
200359e8 L |
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 | |
b2a00c89 | 48 | bid64_mul (UINT64 * pres, UINT64 * px, |
200359e8 L |
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 | |
b2a00c89 | 56 | bid64_mul (UINT64 x, |
200359e8 L |
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; | |
b2a00c89 | 65 | int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy, |
200359e8 L |
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 | |
b2a00c89 | 96 | BID_RETURN (coefficient_x & QUIET_MASK64); |
200359e8 L |
97 | } |
98 | // x is Infinity? | |
99 | if ((x & INFINITY_MASK64) == INFINITY_MASK64) { | |
100 | // check if y is 0 | |
b2a00c89 L |
101 | if (((y & INFINITY_MASK64) != INFINITY_MASK64) |
102 | && !coefficient_y) { | |
200359e8 L |
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 | |
b2a00c89 | 112 | BID_RETURN (coefficient_y & QUIET_MASK64); |
200359e8 L |
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 | |
b2a00c89 | 141 | BID_RETURN (coefficient_y & QUIET_MASK64); |
200359e8 L |
142 | } |
143 | // y is Infinity? | |
144 | if ((y & INFINITY_MASK64) == INFINITY_MASK64) { | |
145 | // check if x is 0 | |
b2a00c89 | 146 | if (!coefficient_x) { |
200359e8 L |
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 | |
b2a00c89 L |
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 | |
200359e8 L |
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) { | |
b2a00c89 | 233 | __add_128_64 (PU, P, round_const_table[rmode][extra_digits]); |
200359e8 | 234 | if (__unsigned_compare_ge_128 |
b2a00c89 | 235 | (PU, power10_table_128[extra_digits + 16])) |
200359e8 L |
236 | uf_status = 0; |
237 | } | |
238 | extra_digits -= final_exponent; | |
239 | final_exponent = 0; | |
240 | ||
241 | if (extra_digits > 17) { | |
b2a00c89 | 242 | __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]); |
200359e8 | 243 | |
b2a00c89 | 244 | amount = recip_scale[16]; |
200359e8 L |
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; | |
b2a00c89 | 255 | if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1] |
200359e8 | 256 | || (Q_low.w[1] == |
b2a00c89 | 257 | reciprocals10_128[16].w[1] |
200359e8 | 258 | && Q_low.w[0] >= |
b2a00c89 | 259 | reciprocals10_128[16].w[0]))) { |
200359e8 L |
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 | |
b2a00c89 | 284 | __add_128_64 (P, P, round_const_table[rmode][extra_digits]); |
200359e8 L |
285 | |
286 | // get P*(2^M[extra_digits])/10^extra_digits | |
287 | __mul_128x128_full (Q_high, Q_low, P, | |
b2a00c89 | 288 | reciprocals10_128[extra_digits]); |
200359e8 L |
289 | |
290 | // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128 | |
b2a00c89 | 291 | amount = recip_scale[extra_digits]; |
200359e8 L |
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 | |
b2a00c89 | 302 | // is exactly .5 |
200359e8 L |
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 | |
b2a00c89 L |
311 | && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
312 | || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
200359e8 | 313 | && Q_low.w[0] < |
b2a00c89 | 314 | reciprocals10_128[extra_digits].w[0]))) { |
200359e8 L |
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 | |
b2a00c89 L |
331 | && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
332 | || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
200359e8 | 333 | && Q_low.w[0] < |
b2a00c89 | 334 | reciprocals10_128[extra_digits].w[0]))) |
200359e8 L |
335 | status = EXACT_STATUS; |
336 | break; | |
337 | case ROUNDING_DOWN: | |
338 | case ROUNDING_TO_ZERO: | |
339 | if (!remainder_h | |
b2a00c89 L |
340 | && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1] |
341 | || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1] | |
200359e8 | 342 | && Q_low.w[0] < |
b2a00c89 | 343 | reciprocals10_128[extra_digits].w[0]))) |
200359e8 L |
344 | status = EXACT_STATUS; |
345 | break; | |
346 | default: | |
347 | // round up | |
348 | __add_carry_out (Stemp.w[0], CY, Q_low.w[0], | |
b2a00c89 | 349 | reciprocals10_128[extra_digits].w[0]); |
200359e8 | 350 | __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1], |
b2a00c89 | 351 | reciprocals10_128[extra_digits].w[1], CY); |
200359e8 L |
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 | } |