]>
git.ipfire.org Git - thirdparty/gcc.git/blob - libgcc/config/libbid/bid128_quantize.c
1 /* Copyright (C) 2007, 2009 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 3, or (at your option) any later
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
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.
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/>. */
25 #include "bid_internal.h"
27 BID128_FUNCTION_ARG2 (bid128_quantize
, x
, y
)
30 UINT128 CX
, CY
, T
, CX2
, CR
, Stemp
, res
, REM_H
, C2N
;
31 UINT64 sign_x
, sign_y
, remainder_h
, carry
, CY64
, valid_x
;
33 int exponent_x
, exponent_y
, digits_x
, extra_digits
, amount
;
34 int expon_diff
, total_digits
, bin_expon_cx
, rmode
, status
;
36 valid_x
= unpack_BID128_value (&sign_x
, &exponent_x
, &CX
, x
);
38 // unpack arguments, check for NaN or Infinity
39 if (!unpack_BID128_value (&sign_y
, &exponent_y
, &CY
, y
)) {
41 #ifdef SET_STATUS_FLAGS
42 if ((x
.w
[1] & SNAN_MASK64
) == SNAN_MASK64
) // y is sNaN
43 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
47 if ((y
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
48 #ifdef SET_STATUS_FLAGS
49 if ((y
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) {
51 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
54 if ((x
.w
[1] & 0x7c00000000000000ull
) != 0x7c00000000000000ull
) {
55 res
.w
[1] = CY
.w
[1] & QUIET_MASK64
;
58 res
.w
[1] = CX
.w
[1] & QUIET_MASK64
;
64 if ((y
.w
[1] & 0x7800000000000000ull
) == 0x7800000000000000ull
) {
65 // check if x is not Inf.
66 if (((x
.w
[1] & 0x7c00000000000000ull
) < 0x7800000000000000ull
)) {
68 #ifdef SET_STATUS_FLAGS
70 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
72 res
.w
[1] = 0x7c00000000000000ull
;
76 if (((x
.w
[1] & 0x7c00000000000000ull
) <= 0x7800000000000000ull
)) {
77 res
.w
[1] = CX
.w
[1] & QUIET_MASK64
;
86 // test if x is NaN or Inf
87 if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7800000000000000ull
) {
88 #ifdef SET_STATUS_FLAGS
90 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
92 res
.w
[1] = 0x7c00000000000000ull
;
95 } else if ((x
.w
[1] & 0x7c00000000000000ull
) == 0x7c00000000000000ull
) {
96 if ((x
.w
[1] & 0x7e00000000000000ull
) == 0x7e00000000000000ull
) {
97 #ifdef SET_STATUS_FLAGS
99 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
102 res
.w
[1] = CX
.w
[1] & QUIET_MASK64
;
106 if (!CX
.w
[1] && !CX
.w
[0]) {
107 get_BID128_very_fast (&res
, sign_x
, exponent_y
, CX
);
111 // get number of decimal digits in coefficient_x
113 tempx
.d
= (float) CX
.w
[1];
114 bin_expon_cx
= ((tempx
.i
>> 23) & 0xff) - 0x7f + 64;
116 tempx
.d
= (float) CX
.w
[0];
117 bin_expon_cx
= ((tempx
.i
>> 23) & 0xff) - 0x7f;
120 digits_x
= estimate_decimal_digits
[bin_expon_cx
];
121 if (CX
.w
[1] > power10_table_128
[digits_x
].w
[1]
122 || (CX
.w
[1] == power10_table_128
[digits_x
].w
[1]
123 && CX
.w
[0] >= power10_table_128
[digits_x
].w
[0]))
126 expon_diff
= exponent_x
- exponent_y
;
127 total_digits
= digits_x
+ expon_diff
;
129 if ((UINT32
) total_digits
<= 34) {
130 if (expon_diff
>= 0) {
131 T
= power10_table_128
[expon_diff
];
132 __mul_128x128_low (CX2
, T
, CX
);
133 get_BID128_very_fast (&res
, sign_x
, exponent_y
, CX2
);
136 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
137 #ifndef IEEE_ROUND_NEAREST
139 if (sign_x
&& (unsigned) (rmode
- 1) < 2)
147 // must round off -expon_diff digits
148 extra_digits
= -expon_diff
;
149 __add_128_128 (CX
, CX
, round_const_table_128
[rmode
][extra_digits
]);
151 // get P*(2^M[extra_digits])/10^extra_digits
152 __mul_128x128_to_256 (CT
, CX
, reciprocals10_128
[extra_digits
]);
154 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
155 amount
= recip_scale
[extra_digits
];
160 CR
.w
[0] = CX2
.w
[1] >> (amount
- 64);
162 __shr_128 (CR
, CX2
, amount
);
165 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
166 #ifndef IEEE_ROUND_NEAREST
170 // check whether fractional part of initial_P/10^extra_digits is
171 // exactly .5 this is the same as fractional part of
172 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
176 remainder_h
= CX2
.w
[0] | (CX2
.w
[1] << (128 - amount
));
178 remainder_h
= CX2
.w
[0] << (64 - amount
);
180 // test whether fractional part is 0
182 && (CT
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
183 || (CT
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
184 && CT
.w
[0] < reciprocals10_128
[extra_digits
].w
[0]))) {
190 #ifdef SET_STATUS_FLAGS
191 status
= INEXACT_EXCEPTION
;
195 REM_H
.w
[1] = (CX2
.w
[1] << (128 - amount
));
196 REM_H
.w
[0] = CX2
.w
[0];
198 REM_H
.w
[1] = CX2
.w
[0] << (64 - amount
);
203 case ROUNDING_TO_NEAREST
:
204 case ROUNDING_TIES_AWAY
:
205 // test whether fractional part is 0
206 if (REM_H
.w
[1] == 0x8000000000000000ull
&& !REM_H
.w
[0]
207 && (CT
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
208 || (CT
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
209 && CT
.w
[0] < reciprocals10_128
[extra_digits
].w
[0])))
210 status
= EXACT_STATUS
;
213 case ROUNDING_TO_ZERO
:
214 if (!(REM_H
.w
[1] | REM_H
.w
[0])
215 && (CT
.w
[1] < reciprocals10_128
[extra_digits
].w
[1]
216 || (CT
.w
[1] == reciprocals10_128
[extra_digits
].w
[1]
217 && CT
.w
[0] < reciprocals10_128
[extra_digits
].w
[0])))
218 status
= EXACT_STATUS
;
222 __add_carry_out (Stemp
.w
[0], CY64
, CT
.w
[0],
223 reciprocals10_128
[extra_digits
].w
[0]);
224 __add_carry_in_out (Stemp
.w
[1], carry
, CT
.w
[1],
225 reciprocals10_128
[extra_digits
].w
[1], CY64
);
228 C2N
.w
[0] = ((UINT64
) 1) << amount
;
229 REM_H
.w
[0] = REM_H
.w
[1] >> (64 - amount
);
232 C2N
.w
[1] = ((UINT64
) 1) << (amount
- 64);
234 REM_H
.w
[1] >>= (128 - amount
);
237 if (REM_H
.w
[0] < carry
)
239 if (__unsigned_compare_ge_128 (REM_H
, C2N
))
240 status
= EXACT_STATUS
;
243 __set_status_flags (pfpsf
, status
);
246 get_BID128_very_fast (&res
, sign_x
, exponent_y
, CR
);
249 if (total_digits
< 0) {
250 CR
.w
[1] = CR
.w
[0] = 0;
251 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
252 #ifndef IEEE_ROUND_NEAREST
254 if (sign_x
&& (unsigned) (rmode
- 1) < 2)
256 if (rmode
== ROUNDING_UP
)
260 #ifdef SET_STATUS_FLAGS
261 __set_status_flags (pfpsf
, INEXACT_EXCEPTION
);
263 get_BID128_very_fast (&res
, sign_x
, exponent_y
, CR
);
266 // else more than 34 digits in coefficient
267 #ifdef SET_STATUS_FLAGS
268 __set_status_flags (pfpsf
, INVALID_EXCEPTION
);
270 res
.w
[1] = 0x7c00000000000000ull
;