1 /* Copyright (C) 2007-2023 Free Software Foundation, Inc.
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
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
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
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/>. */
23
24 #include "bid_internal.h"
25
26 #define MAX_FORMAT_DIGITS 16
27 #define DECIMAL_EXPONENT_BIAS 398
28 #define MAX_DECIMAL_EXPONENT 767
29
30 #if DECIMAL_CALL_BY_REFERENCE
31
32 void
33 bid64_quantize (UINT64 * pres, UINT64 * px,
34 UINT64 *
35 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
36 _EXC_INFO_PARAM) {
37 UINT64 x, y;
38 #else
39
40 UINT64
41 bid64_quantize (UINT64 x,
42 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
43 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
44 #endif
45 UINT128 CT;
46 UINT64 sign_x, sign_y, coefficient_x, coefficient_y, remainder_h, C64,
47 valid_x;
48 UINT64 tmp, carry, res;
49 int_float tempx;
50 int exponent_x, exponent_y, digits_x, extra_digits, amount, amount2;
51 int expon_diff, total_digits, bin_expon_cx;
52 unsigned rmode, status;
53
54 #if DECIMAL_CALL_BY_REFERENCE
55 #if !DECIMAL_GLOBAL_ROUNDING
56 _IDEC_round rnd_mode = *prnd_mode;
57 #endif
58 x = *px;
59 y = *py;
60 #endif
61
62 valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
63 // unpack arguments, check for NaN or Infinity
64 if (!unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y)) {
65 // Inf. or NaN or 0
66 #ifdef SET_STATUS_FLAGS
67 if ((x & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
68 __set_status_flags (pfpsf, INVALID_EXCEPTION);
69 #endif
70
71 // x=Inf, y=Inf?
72 if (((coefficient_x << 1) == 0xf000000000000000ull)
73 && ((coefficient_y << 1) == 0xf000000000000000ull)) {
74 res = coefficient_x;
75 BID_RETURN (res);
76 }
77 // Inf or NaN?
78 if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
79 #ifdef SET_STATUS_FLAGS
80 if (((y & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
81 || (((y & 0x7c00000000000000ull) == 0x7800000000000000ull) && //Inf
82 ((x & 0x7c00000000000000ull) < 0x7800000000000000ull)))
83 __set_status_flags (pfpsf, INVALID_EXCEPTION);
84 #endif
85 if ((y & NAN_MASK64) != NAN_MASK64)
86 coefficient_y = 0;
87 if ((x & NAN_MASK64) != NAN_MASK64) {
88 res = 0x7c00000000000000ull | (coefficient_y & QUIET_MASK64);
89 if (((y & NAN_MASK64) != NAN_MASK64) && ((x & NAN_MASK64) == 0x7800000000000000ull))
90 res = x;
91 BID_RETURN (res);
92 }
93 }
94 }
95 // unpack arguments, check for NaN or Infinity
96 if (!valid_x) {
97 // x is Inf. or NaN or 0
98
99 // Inf or NaN?
100 if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
101 #ifdef SET_STATUS_FLAGS
102 if (((x & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
103 || ((x & 0x7c00000000000000ull) == 0x7800000000000000ull)) //Inf
104 __set_status_flags (pfpsf, INVALID_EXCEPTION);
105 #endif
106 if ((x & NAN_MASK64) != NAN_MASK64)
107 coefficient_x = 0;
108 res = 0x7c00000000000000ull | (coefficient_x & QUIET_MASK64);
109 BID_RETURN (res);
110 }
111
112 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0);
113 BID_RETURN (res);
114 }
115 // get number of decimal digits in coefficient_x
116 tempx.d = (float) coefficient_x;
117 bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
118 digits_x = estimate_decimal_digits[bin_expon_cx];
119 if (coefficient_x >= power10_table_128[digits_x].w[0])
120 digits_x++;
121
122 expon_diff = exponent_x - exponent_y;
123 total_digits = digits_x + expon_diff;
124
125 // check range of scaled coefficient
126 if ((UINT32) (total_digits + 1) <= 17) {
127 if (expon_diff >= 0) {
128 coefficient_x *= power10_table_128[expon_diff].w[0];
129 res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x);
130 BID_RETURN (res);
131 }
132 // must round off -expon_diff digits
133 extra_digits = -expon_diff;
134 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
135 #ifndef IEEE_ROUND_NEAREST
136 rmode = rnd_mode;
137 if (sign_x && (unsigned) (rmode - 1) < 2)
138 rmode = 3 - rmode;
139 #else
140 rmode = 0;
141 #endif
142 #else
143 rmode = 0;
144 #endif
145 coefficient_x += round_const_table[rmode][extra_digits];
146
147 // get P*(2^M[extra_digits])/10^extra_digits
148 __mul_64x64_to_128 (CT, coefficient_x,
149 reciprocals10_64[extra_digits]);
150
151 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
152 amount = short_recip_scale[extra_digits];
153 C64 = CT.w[1] >> amount;
154 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
155 #ifndef IEEE_ROUND_NEAREST
156 if (rnd_mode == 0)
157 #endif
158 if (C64 & 1) {
159 // check whether fractional part of initial_P/10^extra_digits
160 // is exactly .5
161 // this is the same as fractional part of
162 // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
163
164 // get remainder
165 amount2 = 64 - amount;
166 remainder_h = 0;
167 remainder_h--;
168 remainder_h >>= amount2;
169 remainder_h = remainder_h & CT.w[1];
170
171 // test whether fractional part is 0
172 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
173 C64--;
174 }
175 }
176 #endif
177
178 #ifdef SET_STATUS_FLAGS
179 status = INEXACT_EXCEPTION;
180 // get remainder
181 remainder_h = CT.w[1] << (64 - amount);
182 switch (rmode) {
183 case ROUNDING_TO_NEAREST:
184 case ROUNDING_TIES_AWAY:
185 // test whether fractional part is 0
186 if ((remainder_h == 0x8000000000000000ull)
187 && (CT.w[0] < reciprocals10_64[extra_digits]))
188 status = EXACT_STATUS;
189 break;
190 case ROUNDING_DOWN:
191 case ROUNDING_TO_ZERO:
192 if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
193 status = EXACT_STATUS;
194 //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
195 break;
196 default:
197 // round up
198 __add_carry_out (tmp, carry, CT.w[0],
199 reciprocals10_64[extra_digits]);
200 if ((remainder_h >> (64 - amount)) + carry >=
201 (((UINT64) 1) << amount))
202 status = EXACT_STATUS;
203 break;
204 }
205 __set_status_flags (pfpsf, status);
206 #endif
207
208 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64);
209 BID_RETURN (res);
210 }
211
212 if (total_digits < 0) {
213 #ifdef SET_STATUS_FLAGS
214 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
215 #endif
216 C64 = 0;
217 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
218 #ifndef IEEE_ROUND_NEAREST
219 rmode = rnd_mode;
220 if (sign_x && (unsigned) (rmode - 1) < 2)
221 rmode = 3 - rmode;
222 if (rmode == ROUNDING_UP)
223 C64 = 1;
224 #endif
225 #endif
226 res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64);
227 BID_RETURN (res);
228 }
229 // else more than 16 digits in coefficient
230 #ifdef SET_STATUS_FLAGS
231 __set_status_flags (pfpsf, INVALID_EXCEPTION);
232 #endif
233 res = 0x7c00000000000000ull;
234 BID_RETURN (res);
235
236 }