1 /* Copyright (C) 2007 Free Software Foundation, Inc.
3 This file is part of GCC.
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 2, or (at your option) any later
10 In addition to the permissions in the GNU General Public License, the
11 Free Software Foundation gives you unlimited permission to link the
12 compiled version of this file into combinations with other programs,
13 and to distribute those combinations without any restriction coming
14 from the use of this file. (The General Public License restrictions
15 do apply in other respects; for example, they cover modification of
16 the file, and distribution when not linked into a combine
19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24 You should have received a copy of the GNU General Public License
25 along with GCC; see the file COPYING. If not, write to the Free
26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
29 /*****************************************************************************
31 *****************************************************************************
33 * Algorithm description:
35 * if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed
37 * return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias,
38 * coefficient_x*coefficient_y)
40 * get long product: coefficient_x*coefficient_y
41 * determine number of digits to round off (extra_digits)
42 * rounding is performed as a 128x128-bit multiplication by
43 * 2^M[extra_digits]/10^extra_digits, followed by a shift
44 * M[extra_digits] is sufficiently large for required accuracy
46 ****************************************************************************/
48 #include "bid_internal.h"
50 #if DECIMAL_CALL_BY_REFERENCE
53 bid64_mul (UINT64
* pres
, UINT64
* px
,
55 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
62 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
63 _EXC_MASKS_PARAM _EXC_INFO_PARAM
) {
65 UINT128 P
, PU
, C128
, Q_high
, Q_low
, Stemp
;
66 UINT64 sign_x
, sign_y
, coefficient_x
, coefficient_y
;
67 UINT64 C64
, remainder_h
, carry
, CY
, res
;
68 UINT64 valid_x
, valid_y
;
69 int_double tempx
, tempy
;
70 int extra_digits
, exponent_x
, exponent_y
, bin_expon_cx
, bin_expon_cy
,
72 int rmode
, digits_p
, bp
, amount
, amount2
, final_exponent
, round_up
;
73 unsigned status
, uf_status
;
75 #if DECIMAL_CALL_BY_REFERENCE
76 #if !DECIMAL_GLOBAL_ROUNDING
77 _IDEC_round rnd_mode
= *prnd_mode
;
83 valid_x
= unpack_BID64 (&sign_x
, &exponent_x
, &coefficient_x
, x
);
84 valid_y
= unpack_BID64 (&sign_y
, &exponent_y
, &coefficient_y
, y
);
86 // unpack arguments, check for NaN or Infinity
89 #ifdef SET_STATUS_FLAGS
90 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
91 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
96 if ((x
& NAN_MASK64
) == NAN_MASK64
) {
97 #ifdef SET_STATUS_FLAGS
98 if ((x
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
99 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
101 BID_RETURN (coefficient_x
& QUIET_MASK64
);
104 if ((x
& INFINITY_MASK64
) == INFINITY_MASK64
) {
106 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
)
108 #ifdef SET_STATUS_FLAGS
109 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
112 BID_RETURN (NAN_MASK64
);
115 if ((y
& NAN_MASK64
) == NAN_MASK64
)
116 // y==NaN , return NaN
117 BID_RETURN (coefficient_y
& QUIET_MASK64
);
118 // otherwise return +/-Inf
119 BID_RETURN (((x
^ y
) & 0x8000000000000000ull
) | INFINITY_MASK64
);
122 if (((y
& INFINITY_MASK64
) != INFINITY_MASK64
)) {
123 if ((y
& SPECIAL_ENCODING_MASK64
) == SPECIAL_ENCODING_MASK64
)
124 exponent_y
= ((UINT32
) (y
>> 51)) & 0x3ff;
126 exponent_y
= ((UINT32
) (y
>> 53)) & 0x3ff;
127 sign_y
= y
& 0x8000000000000000ull
;
129 exponent_x
+= exponent_y
- DECIMAL_EXPONENT_BIAS
;
130 if (exponent_x
> DECIMAL_MAX_EXPON_64
)
131 exponent_x
= DECIMAL_MAX_EXPON_64
;
132 else if (exponent_x
< 0)
134 BID_RETURN ((sign_x
^ sign_y
) | (((UINT64
) exponent_x
) << 53));
141 if ((y
& NAN_MASK64
) == NAN_MASK64
) {
142 #ifdef SET_STATUS_FLAGS
143 if ((y
& SNAN_MASK64
) == SNAN_MASK64
) // sNaN
144 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
146 BID_RETURN (coefficient_y
& QUIET_MASK64
);
149 if ((y
& INFINITY_MASK64
) == INFINITY_MASK64
) {
151 if (!coefficient_x
) {
152 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
154 BID_RETURN (NAN_MASK64
);
156 // otherwise return +/-Inf
157 BID_RETURN (((x
^ y
) & 0x8000000000000000ull
) | INFINITY_MASK64
);
160 exponent_x
+= exponent_y
- DECIMAL_EXPONENT_BIAS
;
161 if (exponent_x
> DECIMAL_MAX_EXPON_64
)
162 exponent_x
= DECIMAL_MAX_EXPON_64
;
163 else if (exponent_x
< 0)
165 BID_RETURN ((sign_x
^ sign_y
) | (((UINT64
) exponent_x
) << 53));
167 //--- get number of bits in the coefficients of x and y ---
168 // version 2 (original)
169 tempx
.d
= (double) coefficient_x
;
170 bin_expon_cx
= ((tempx
.i
& MASK_BINARY_EXPONENT
) >> 52);
171 tempy
.d
= (double) coefficient_y
;
172 bin_expon_cy
= ((tempy
.i
& MASK_BINARY_EXPONENT
) >> 52);
174 // magnitude estimate for coefficient_x*coefficient_y is
175 // 2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)
176 bin_expon_product
= bin_expon_cx
+ bin_expon_cy
;
178 // check if coefficient_x*coefficient_y<2^(10*k+3)
179 // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1
180 if (bin_expon_product
< UPPER_EXPON_LIMIT
+ 2 * BINARY_EXPONENT_BIAS
) {
182 C64
= coefficient_x
* coefficient_y
;
185 get_BID64_small_mantissa (sign_x
^ sign_y
,
186 exponent_x
+ exponent_y
-
187 DECIMAL_EXPONENT_BIAS
, C64
, rnd_mode
,
192 // get 128-bit product: coefficient_x*coefficient_y
193 __mul_64x64_to_128 (P
, coefficient_x
, coefficient_y
);
195 // tighten binary range of P: leading bit is 2^bp
196 // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1
197 bin_expon_product
-= 2 * BINARY_EXPONENT_BIAS
;
199 __tight_bin_range_128 (bp
, P
, bin_expon_product
);
201 // get number of decimal digits in the product
202 digits_p
= estimate_decimal_digits
[bp
];
203 if (!(__unsigned_compare_gt_128 (power10_table_128
[digits_p
], P
)))
204 digits_p
++; // if power10_table_128[digits_p] <= P
206 // determine number of decimal digits to be rounded out
207 extra_digits
= digits_p
- MAX_FORMAT_DIGITS
;
209 exponent_x
+ exponent_y
+ extra_digits
- DECIMAL_EXPONENT_BIAS
;
211 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
212 #ifndef IEEE_ROUND_NEAREST
214 if (sign_x
^ sign_y
&& (unsigned) (rmode
- 1) < 2)
224 if (((unsigned) final_exponent
) >= 3 * 256) {
225 if (final_exponent
< 0) {
227 if (final_exponent
+ 16 < 0) {
228 res
= sign_x
^ sign_y
;
229 __set_status_flags (pfpsf
,
230 UNDERFLOW_EXCEPTION
| INEXACT_EXCEPTION
);
231 if (rmode
== ROUNDING_UP
)
236 uf_status
= UNDERFLOW_EXCEPTION
;
237 if (final_exponent
== -1) {
238 __add_128_64 (PU
, P
, round_const_table
[rmode
][extra_digits
]);
239 if (__unsigned_compare_ge_128
240 (PU
, power10_table_128
[extra_digits
+ 16]))
243 extra_digits
-= final_exponent
;
246 if (extra_digits
> 17) {
247 __mul_128x128_full (Q_high
, Q_low
, P
, reciprocals10_128
[16]);
249 amount
= recip_scale
[16];
250 __shr_128 (P
, Q_high
, amount
);
253 amount2
= 64 - amount
;
256 remainder_h
>>= amount2
;
257 remainder_h
= remainder_h
& Q_high
.w
[0];
260 if (remainder_h
|| (Q_low
.w
[1] > reciprocals10_128
[16].w
[1]
262 reciprocals10_128
[16].w
[1]
264 reciprocals10_128
[16].w
[0]))) {
266 __set_status_flags (pfpsf
,
267 UNDERFLOW_EXCEPTION
|
269 P
.w
[0] = (P
.w
[0] << 3) + (P
.w
[0] << 1);
276 fast_get_BID64_check_OF (sign_x
^ sign_y
, final_exponent
,
277 1000000000000000ull, rnd_mode
,
284 if (extra_digits
> 0) {
285 // will divide by 10^(digits_p - 16)
287 // add a constant to P, depending on rounding mode
288 // 0.5*10^(digits_p - 16) for round-to-nearest
289 __add_128_64 (P
, P
, round_const_table
[rmode
][extra_digits
]);
291 // get P*(2^M[extra_digits])/10^extra_digits
292 __mul_128x128_full (Q_high
, Q_low
, P
,
293 reciprocals10_128
[extra_digits
]);
295 // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
296 amount
= recip_scale
[extra_digits
];
297 __shr_128 (C128
, Q_high
, amount
);
299 C64
= __low_64 (C128
);
301 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
302 #ifndef IEEE_ROUND_NEAREST
303 if (rmode
== 0) //ROUNDING_TO_NEAREST
305 if ((C64
& 1) && !round_up
) {
306 // check whether fractional part of initial_P/10^extra_digits
308 // this is the same as fractional part of
309 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
312 remainder_h
= Q_high
.w
[0] << (64 - amount
);
314 // test whether fractional part is 0
316 && (Q_low
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
317 || (Q_low
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
319 reciprocals10_128
[extra_digits
].w
[0]))) {
325 #ifdef SET_STATUS_FLAGS
326 status
= INEXACT_EXCEPTION
| uf_status
;
329 remainder_h
= Q_high
.w
[0] << (64 - amount
);
332 case ROUNDING_TO_NEAREST
:
333 case ROUNDING_TIES_AWAY
:
334 // test whether fractional part is 0
335 if (remainder_h
== 0x8000000000000000ull
336 && (Q_low
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
337 || (Q_low
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
339 reciprocals10_128
[extra_digits
].w
[0])))
340 status
= EXACT_STATUS
;
343 case ROUNDING_TO_ZERO
:
345 && (Q_low
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
346 || (Q_low
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
348 reciprocals10_128
[extra_digits
].w
[0])))
349 status
= EXACT_STATUS
;
353 __add_carry_out (Stemp
.w
[0], CY
, Q_low
.w
[0],
354 reciprocals10_128
[extra_digits
].w
[0]);
355 __add_carry_in_out (Stemp
.w
[1], carry
, Q_low
.w
[1],
356 reciprocals10_128
[extra_digits
].w
[1], CY
);
357 if ((remainder_h
>> (64 - amount
)) + carry
>=
358 (((UINT64
) 1) << amount
))
359 status
= EXACT_STATUS
;
362 __set_status_flags (pfpsf
, status
);
365 // convert to BID and return
367 fast_get_BID64_check_OF (sign_x
^ sign_y
, final_exponent
, C64
,
371 // go to convert_format and exit
374 get_BID64 (sign_x
^ sign_y
,
375 exponent_x
+ exponent_y
- DECIMAL_EXPONENT_BIAS
, C64
,