Mercurial > hg > CbC > CbC_gcc

comparison libgcc/config/libbid/bid128_quantize.c @ 0:a06113de4d67

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first commit
author kent <kent@cr.ie.u-ryukyu.ac.jp>
date Fri, 17 Jul 2009 14:47:48 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
-1:000000000000 0:a06113de4d67
1 /* Copyright (C) 2007, 2009 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 #define BID_128RES
25 #include "bid_internal.h"
26
27 BID128_FUNCTION_ARG2 (bid128_quantize, x, y)
28
29 UINT256 CT;
30 UINT128 CX, CY, T, CX2, CR, Stemp, res, REM_H, C2N;
31 UINT64 sign_x, sign_y, remainder_h, carry, CY64, valid_x;
32 int_float tempx;
33 int exponent_x, exponent_y, digits_x, extra_digits, amount;
34 int expon_diff, total_digits, bin_expon_cx, rmode, status;
35
36 valid_x = unpack_BID128_value (&sign_x, &exponent_x, &CX, x);
37
38 // unpack arguments, check for NaN or Infinity
39 if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
40 // y is Inf. or NaN
41 #ifdef SET_STATUS_FLAGS
42 if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
43 __set_status_flags (pfpsf, INVALID_EXCEPTION);
44 #endif
45
46 // test if y is NaN
47 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
48 #ifdef SET_STATUS_FLAGS
49 if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
50 // set status flags
51 __set_status_flags (pfpsf, INVALID_EXCEPTION);
52 }
53 #endif
54 if ((x.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull) {
55 res.w[1] = CY.w[1] & QUIET_MASK64;
56 res.w[0] = CY.w[0];
57 } else {
58 res.w[1] = CX.w[1] & QUIET_MASK64;
59 res.w[0] = CX.w[0];
60 }
61 BID_RETURN (res);
62 }
63 // y is Infinity?
64 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
65 // check if x is not Inf.
66 if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) {
67 // return NaN
68 #ifdef SET_STATUS_FLAGS
69 // set status flags
70 __set_status_flags (pfpsf, INVALID_EXCEPTION);
71 #endif
72 res.w[1] = 0x7c00000000000000ull;
73 res.w[0] = 0;
74 BID_RETURN (res);
75 } else
76 if (((x.w[1] & 0x7c00000000000000ull) <= 0x7800000000000000ull)) {
77 res.w[1] = CX.w[1] & QUIET_MASK64;
78 res.w[0] = CX.w[0];
79 BID_RETURN (res);
80 }
81 }
82
83 }
84
85 if (!valid_x) {
86 // test if x is NaN or Inf
87 if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) {
88 #ifdef SET_STATUS_FLAGS
89 // set status flags
90 __set_status_flags (pfpsf, INVALID_EXCEPTION);
91 #endif
92 res.w[1] = 0x7c00000000000000ull;
93 res.w[0] = 0;
94 BID_RETURN (res);
95 } else if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
96 if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
97 #ifdef SET_STATUS_FLAGS
98 // set status flags
99 __set_status_flags (pfpsf, INVALID_EXCEPTION);
100 #endif
101 }
102 res.w[1] = CX.w[1] & QUIET_MASK64;
103 res.w[0] = CX.w[0];
104 BID_RETURN (res);
105 }
106 if (!CX.w[1] && !CX.w[0]) {
107 get_BID128_very_fast (&res, sign_x, exponent_y, CX);
108 BID_RETURN (res);
109 }
110 }
111 // get number of decimal digits in coefficient_x
112 if (CX.w[1]) {
113 tempx.d = (float) CX.w[1];
114 bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f + 64;
115 } else {
116 tempx.d = (float) CX.w[0];
117 bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
118 }
119
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]))
124 digits_x++;
125
126 expon_diff = exponent_x - exponent_y;
127 total_digits = digits_x + expon_diff;
128
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);
134 BID_RETURN (res);
135 }
136 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
137 #ifndef IEEE_ROUND_NEAREST
138 rmode = rnd_mode;
139 if (sign_x && (unsigned) (rmode - 1) < 2)
140 rmode = 3 - rmode;
141 #else
142 rmode = 0;
143 #endif
144 #else
145 rmode = 0;
146 #endif
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]);
150
151 // get P*(2^M[extra_digits])/10^extra_digits
152 __mul_128x128_to_256 (CT, CX, reciprocals10_128[extra_digits]);
153
154 // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
155 amount = recip_scale[extra_digits];
156 CX2.w[0] = CT.w[2];
157 CX2.w[1] = CT.w[3];
158 if (amount >= 64) {
159 CR.w[1] = 0;
160 CR.w[0] = CX2.w[1] >> (amount - 64);
161 } else {
162 __shr_128 (CR, CX2, amount);
163 }
164
165 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
166 #ifndef IEEE_ROUND_NEAREST
167 if (rnd_mode == 0)
168 #endif
169 if (CR.w[0] & 1) {
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
173
174 // get remainder
175 if (amount >= 64) {
176 remainder_h = CX2.w[0] | (CX2.w[1] << (128 - amount));
177 } else
178 remainder_h = CX2.w[0] << (64 - amount);
179
180 // test whether fractional part is 0
181 if (!remainder_h
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]))) {
185 CR.w[0]--;
186 }
187 }
188 #endif
189
190 #ifdef SET_STATUS_FLAGS
191 status = INEXACT_EXCEPTION;
192
193 // get remainder
194 if (amount >= 64) {
195 REM_H.w[1] = (CX2.w[1] << (128 - amount));
196 REM_H.w[0] = CX2.w[0];
197 } else {
198 REM_H.w[1] = CX2.w[0] << (64 - amount);
199 REM_H.w[0] = 0;
200 }
201
202 switch (rmode) {
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;
211 break;
212 case ROUNDING_DOWN:
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;
219 break;
220 default:
221 // round up
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);
226 if (amount < 64) {
227 C2N.w[1] = 0;
228 C2N.w[0] = ((UINT64) 1) << amount;
229 REM_H.w[0] = REM_H.w[1] >> (64 - amount);
230 REM_H.w[1] = 0;
231 } else {
232 C2N.w[1] = ((UINT64) 1) << (amount - 64);
233 C2N.w[0] = 0;
234 REM_H.w[1] >>= (128 - amount);
235 }
236 REM_H.w[0] += carry;
237 if (REM_H.w[0] < carry)
238 REM_H.w[1]++;
239 if (__unsigned_compare_ge_128 (REM_H, C2N))
240 status = EXACT_STATUS;
241 }
242
243 __set_status_flags (pfpsf, status);
244
245 #endif
246 get_BID128_very_fast (&res, sign_x, exponent_y, CR);
247 BID_RETURN (res);
248 }
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
253 rmode = rnd_mode;
254 if (sign_x && (unsigned) (rmode - 1) < 2)
255 rmode = 3 - rmode;
256 if (rmode == ROUNDING_UP)
257 CR.w[0] = 1;
258 #endif
259 #endif
260 #ifdef SET_STATUS_FLAGS
261 __set_status_flags (pfpsf, INEXACT_EXCEPTION);
262 #endif
263 get_BID128_very_fast (&res, sign_x, exponent_y, CR);
264 BID_RETURN (res);
265 }
266 // else more than 34 digits in coefficient
267 #ifdef SET_STATUS_FLAGS
268 __set_status_flags (pfpsf, INVALID_EXCEPTION);
269 #endif
270 res.w[1] = 0x7c00000000000000ull;
271 res.w[0] = 0;
272 BID_RETURN (res);
273
274 }