comparison libgcc/config/libbid/bid64_to_bid128.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 #define BID_128RES
25 #include "bid_internal.h"
26
27 /*
28 * Takes a BID64 as input and converts it to a BID128 and returns it.
29 */
30 TYPE0_FUNCTION_ARGTYPE1_NORND (UINT128, bid64_to_bid128, UINT64, x)
31
32 UINT128 new_coeff, res;
33 UINT64 sign_x;
34 int exponent_x;
35 UINT64 coefficient_x;
36
37 if (!unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x)) {
38 if (((x) << 1) >= 0xf000000000000000ull) {
39 #ifdef SET_STATUS_FLAGS
40 if (((x) & SNAN_MASK64) == SNAN_MASK64) // sNaN
41 __set_status_flags (pfpsf, INVALID_EXCEPTION);
42 #endif
43 res.w[0] = (coefficient_x & 0x0003ffffffffffffull);
44 __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
45 res.w[1] |= ((coefficient_x) & 0xfc00000000000000ull);
46 BID_RETURN (res);
47 }
48 }
49
50 new_coeff.w[0] = coefficient_x;
51 new_coeff.w[1] = 0;
52 get_BID128_very_fast (&res, sign_x,
53 exponent_x + DECIMAL_EXPONENT_BIAS_128 -
54 DECIMAL_EXPONENT_BIAS, new_coeff);
55 BID_RETURN (res);
56 } // convert_bid64_to_bid128
57
58
59
60 /*
61 * Takes a BID128 as input and converts it to a BID64 and returns it.
62 */
63 #if DECIMAL_CALL_BY_REFERENCE
64
65 void
66 bid128_to_bid64 (UINT64 * pres,
67 UINT128 *
68 px _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
69 _EXC_INFO_PARAM) {
70 UINT128 x = *px;
71 #else
72
73 UINT64
74 bid128_to_bid64 (UINT128 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
75 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
76 #endif
77 UINT128 CX, T128, TP128, Qh, Ql, Qh1, Stemp, Tmp, Tmp1, CX1;
78 UINT64 sign_x, carry, cy, res;
79 SINT64 D;
80 int_float f64, fx;
81 int exponent_x, extra_digits, amount, bin_expon_cx;
82 unsigned rmode, status, uf_check = 0;
83
84 #if DECIMAL_CALL_BY_REFERENCE
85 #if !DECIMAL_GLOBAL_ROUNDING
86 _IDEC_round rnd_mode = *prnd_mode;
87 #endif
88 #endif
89
90 BID_SWAP128 (x);
91 // unpack arguments, check for NaN or Infinity or 0
92 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
93 if ((x.w[1] << 1) >= 0xf000000000000000ull) {
94 Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull);
95 Tmp.w[0] = CX.w[0];
96 TP128 = reciprocals10_128[18];
97 __mul_128x128_full (Qh, Ql, Tmp, TP128);
98 amount = recip_scale[18];
99 __shr_128 (Tmp, Qh, amount);
100 res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
101 #ifdef SET_STATUS_FLAGS
102 if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // sNaN
103 __set_status_flags (pfpsf, INVALID_EXCEPTION);
104 #endif
105 BID_RETURN_VAL (res);
106 }
107 exponent_x =
108 exponent_x - DECIMAL_EXPONENT_BIAS_128 + DECIMAL_EXPONENT_BIAS;
109 if (exponent_x < 0) {
110 res = sign_x;
111 BID_RETURN_VAL (res);
112 }
113 if (exponent_x > DECIMAL_MAX_EXPON_64)
114 exponent_x = DECIMAL_MAX_EXPON_64;
115 res = sign_x | (((UINT64) exponent_x) << 53);
116 BID_RETURN_VAL (res);
117 }
118
119 if (CX.w[1] || (CX.w[0] >= 10000000000000000ull)) {
120 // find number of digits in coefficient
121 // 2^64
122 f64.i = 0x5f800000;
123 // fx ~ CX
124 fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
125 bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
126 extra_digits = estimate_decimal_digits[bin_expon_cx] - 16;
127 // scale = 38-estimate_decimal_digits[bin_expon_cx];
128 D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];
129 if (D > 0
130 || (!D
131 && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))
132 extra_digits++;
133
134 exponent_x += extra_digits;
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 if (exponent_x < DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS) {
148 uf_check = 1;
149 if (-extra_digits + exponent_x - DECIMAL_EXPONENT_BIAS_128 +
150 DECIMAL_EXPONENT_BIAS + 35 >= 0) {
151 if (exponent_x ==
152 DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS - 1) {
153 T128 = round_const_table_128[rmode][extra_digits];
154 __add_carry_out (CX1.w[0], carry, T128.w[0], CX.w[0]);
155 CX1.w[1] = CX.w[1] + T128.w[1] + carry;
156 if (__unsigned_compare_ge_128
157 (CX1, power10_table_128[extra_digits + 16]))
158 uf_check = 0;
159 }
160 extra_digits =
161 extra_digits + DECIMAL_EXPONENT_BIAS_128 -
162 DECIMAL_EXPONENT_BIAS - exponent_x;
163 exponent_x = DECIMAL_EXPONENT_BIAS_128 - DECIMAL_EXPONENT_BIAS;
164 //uf_check = 2;
165 } else
166 rmode = ROUNDING_TO_ZERO;
167 }
168
169 T128 = round_const_table_128[rmode][extra_digits];
170 __add_carry_out (CX.w[0], carry, T128.w[0], CX.w[0]);
171 CX.w[1] = CX.w[1] + T128.w[1] + carry;
172
173 TP128 = reciprocals10_128[extra_digits];
174 __mul_128x128_full (Qh, Ql, CX, TP128);
175 amount = recip_scale[extra_digits];
176
177 if (amount >= 64) {
178 CX.w[0] = Qh.w[1] >> (amount - 64);
179 CX.w[1] = 0;
180 } else {
181 __shr_128 (CX, Qh, amount);
182 }
183
184 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
185 #ifndef IEEE_ROUND_NEAREST
186 if (!(rmode))
187 #endif
188 if (CX.w[0] & 1) {
189 // check whether fractional part of initial_P/10^ed1 is exactly .5
190
191 // get remainder
192 __shl_128_long (Qh1, Qh, (128 - amount));
193
194 if (!Qh1.w[1] && !Qh1.w[0]
195 && (Ql.w[1] < reciprocals10_128[extra_digits].w[1]
196 || (Ql.w[1] == reciprocals10_128[extra_digits].w[1]
197 && Ql.w[0] < reciprocals10_128[extra_digits].w[0]))) {
198 CX.w[0]--;
199 }
200 }
201 #endif
202
203 {
204 status = INEXACT_EXCEPTION;
205 // get remainder
206 __shl_128_long (Qh1, Qh, (128 - amount));
207
208 switch (rmode) {
209 case ROUNDING_TO_NEAREST:
210 case ROUNDING_TIES_AWAY:
211 // test whether fractional part is 0
212 if (Qh1.w[1] == 0x8000000000000000ull && (!Qh1.w[0])
213 && (Ql.w[1] < reciprocals10_128[extra_digits].w[1]
214 || (Ql.w[1] == reciprocals10_128[extra_digits].w[1]
215 && Ql.w[0] < reciprocals10_128[extra_digits].w[0])))
216 status = EXACT_STATUS;
217 break;
218 case ROUNDING_DOWN:
219 case ROUNDING_TO_ZERO:
220 if ((!Qh1.w[1]) && (!Qh1.w[0])
221 && (Ql.w[1] < reciprocals10_128[extra_digits].w[1]
222 || (Ql.w[1] == reciprocals10_128[extra_digits].w[1]
223 && Ql.w[0] < reciprocals10_128[extra_digits].w[0])))
224 status = EXACT_STATUS;
225 break;
226 default:
227 // round up
228 __add_carry_out (Stemp.w[0], cy, Ql.w[0],
229 reciprocals10_128[extra_digits].w[0]);
230 __add_carry_in_out (Stemp.w[1], carry, Ql.w[1],
231 reciprocals10_128[extra_digits].w[1], cy);
232 __shr_128_long (Qh, Qh1, (128 - amount));
233 Tmp.w[0] = 1;
234 Tmp.w[1] = 0;
235 __shl_128_long (Tmp1, Tmp, amount);
236 Qh.w[0] += carry;
237 if (Qh.w[0] < carry)
238 Qh.w[1]++;
239 if (__unsigned_compare_ge_128 (Qh, Tmp1))
240 status = EXACT_STATUS;
241 }
242
243 if (status != EXACT_STATUS) {
244 if (uf_check)
245 status |= UNDERFLOW_EXCEPTION;
246 #ifdef SET_STATUS_FLAGS
247 __set_status_flags (pfpsf, status);
248 #endif
249 }
250
251
252 }
253
254 }
255
256 res =
257 get_BID64 (sign_x,
258 exponent_x - DECIMAL_EXPONENT_BIAS_128 +
259 DECIMAL_EXPONENT_BIAS, CX.w[0], rnd_mode, pfpsf);
260 BID_RETURN_VAL (res);
261
262 }