Mercurial > hg > CbC > CbC_gcc
comparison libgcc/config/libbid/bid64_quantize.c @ 0:a06113de4d67
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 #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 } |