Mercurial > hg > CbC > CbC_gcc

comparison libdecnumber/bid/bid2dpd_dpd2bid.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 #undef IN_LIBGCC2
25 #include "bid-dpd.h"
26
27 /* get full 64x64bit product */
28 #define __mul_64x64_to_128(P, CX, CY) \
29 { \
30 UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2; \
31 CXH = (CX) >> 32; \
32 CXL = (UINT32)(CX); \
33 CYH = (CY) >> 32; \
34 CYL = (UINT32)(CY); \
35 \
36 PM = CXH*CYL; \
37 PH = CXH*CYH; \
38 PL = CXL*CYL; \
39 PM2 = CXL*CYH; \
40 PH += (PM>>32); \
41 PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \
42 \
43 (P).w[1] = PH + (PM>>32); \
44 (P).w[0] = (PM<<32)+(UINT32)PL; \
45 }
46
47 /* add 64-bit value to 128-bit */
48 #define __add_128_64(R128, A128, B64) \
49 { \
50 UINT64 R64H; \
51 R64H = (A128).w[1]; \
52 (R128).w[0] = (B64) + (A128).w[0]; \
53 if((R128).w[0] < (B64)) R64H ++; \
54 (R128).w[1] = R64H; \
55 }
56
57 /* add 128-bit value to 128-bit (assume no carry-out) */
58 #define __add_128_128(R128, A128, B128) \
59 { \
60 UINT128 Q128; \
61 Q128.w[1] = (A128).w[1]+(B128).w[1]; \
62 Q128.w[0] = (B128).w[0] + (A128).w[0]; \
63 if(Q128.w[0] < (B128).w[0]) Q128.w[1] ++; \
64 (R128).w[1] = Q128.w[1]; \
65 (R128).w[0] = Q128.w[0]; \
66 }
67
68 #define __mul_128x128_high(Q, A, B) \
69 { \
70 UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2; \
71 \
72 __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]); \
73 __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]); \
74 __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]); \
75 __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]); \
76 \
77 __add_128_128(QM, ALBH, AHBL); \
78 __add_128_64(QM2, QM, ALBL.w[1]); \
79 __add_128_64((Q), AHBH, QM2.w[1]); \
80 }
81
82 #include "bid2dpd_dpd2bid.h"
83
84 static const unsigned int dm103[] =
85 { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000 };
86
87 void _bid_to_dpd32 (_Decimal32 *, _Decimal32 *);
88
89 void
90 _bid_to_dpd32 (_Decimal32 *pres, _Decimal32 *px) {
91 unsigned int sign, coefficient_x, exp, dcoeff;
92 unsigned int b2, b1, b0, b01, res;
93 _Decimal32 x = *px;
94
95 sign = (x & 0x80000000);
96 if ((x & 0x60000000ul) == 0x60000000ul) {
97 /* special encodings */
98 if ((x & 0x78000000ul) == 0x78000000ul) {
99 *pres = x; /* NaN or Infinity */
100 return;
101 }
102 /* coefficient */
103 coefficient_x = (x & 0x001ffffful) | 0x00800000ul;
104 if (coefficient_x >= 10000000) coefficient_x = 0;
105 /* get exponent */
106 exp = (x >> 21) & 0xff;
107 } else {
108 exp = (x >> 23) & 0xff;
109 coefficient_x = (x & 0x007ffffful);
110 }
111 b01 = coefficient_x / 1000;
112 b2 = coefficient_x - 1000 * b01;
113 b0 = b01 / 1000;
114 b1 = b01 - 1000 * b0;
115 dcoeff = b2d[b2] | b2d2[b1];
116 if (b0 >= 8) { /* is b0 8 or 9? */
117 res = sign | ((0x600 | ((exp >> 6) << 7) |
118 ((b0 & 1) << 6) | (exp & 0x3f)) << 20) | dcoeff;
119 } else { /* else b0 is 0..7 */
120 res = sign | ((((exp >> 6) << 9) | (b0 << 6) |
121 (exp & 0x3f)) << 20) | dcoeff;
122 }
123 *pres = res;
124 }
125
126 void _dpd_to_bid32 (_Decimal32 *, _Decimal32 *);
127
128 void
129 _dpd_to_bid32 (_Decimal32 *pres, _Decimal32 *px) {
130 unsigned int r;
131 unsigned int sign, exp, bcoeff;
132 UINT64 trailing;
133 unsigned int d0, d1, d2;
134 _Decimal32 x = *px;
135
136 sign = (x & 0x80000000);
137 trailing = (x & 0x000fffff);
138 if ((x & 0x78000000) == 0x78000000) {
139 *pres = x;
140 return;
141 } else { /* normal number */
142 if ((x & 0x60000000) == 0x60000000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
143 d0 = d2b3[((x >> 26) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
144 exp = (x >> 27) & 3; /* exp leading bits are G2..G3 */
145 } else {
146 d0 = d2b3[(x >> 26) & 0x7];
147 exp = (x >> 29) & 3; /* exp loading bits are G0..G1 */
148 }
149 d1 = d2b2[(trailing >> 10) & 0x3ff];
150 d2 = d2b[(trailing) & 0x3ff];
151 bcoeff = d2 + d1 + d0;
152 exp = (exp << 6) + ((x >> 20) & 0x3f);
153 if (bcoeff < (1 << 23)) {
154 r = exp;
155 r <<= 23;
156 r |= (bcoeff | sign);
157 } else {
158 r = exp;
159 r <<= 21;
160 r |= (sign | 0x60000000ul);
161 /* add coeff, without leading bits */
162 r |= (((unsigned int) bcoeff) & 0x1fffff);
163 }
164 }
165 *pres = r;
166 }
167
168 void _bid_to_dpd64 (_Decimal64 *, _Decimal64 *);
169
170 void
171 _bid_to_dpd64 (_Decimal64 *pres, _Decimal64 *px) {
172 UINT64 res;
173 UINT64 sign, comb, exp, B34, B01;
174 UINT64 d103, D61;
175 UINT64 b0, b2, b3, b5;
176 unsigned int b1, b4;
177 UINT64 bcoeff;
178 UINT64 dcoeff;
179 unsigned int yhi, ylo;
180 _Decimal64 x = *px;
181
182 sign = (x & 0x8000000000000000ull);
183 comb = (x & 0x7ffc000000000000ull) >> 51;
184 if ((comb & 0xf00) == 0xf00) {
185 *pres = x;
186 return;
187 } else { /* Normal number */
188 if ((comb & 0xc00) == 0xc00) { /* G0..G1 = 11 -> exp is G2..G11 */
189 exp = (comb) & 0x3ff;
190 bcoeff = (x & 0x0007ffffffffffffull) | 0x0020000000000000ull;
191 } else {
192 exp = (comb >> 2) & 0x3ff;
193 bcoeff = (x & 0x001fffffffffffffull);
194 }
195 D61 = 2305843009ull; /* Floor(2^61 / 10^9) */
196 /* Multiply the binary coefficient by ceil(2^64 / 1000), and take the upper
197 64-bits in order to compute a division by 1000. */
198 yhi = (D61 * (UINT64)(bcoeff >> (UINT64)27)) >> (UINT64)34;
199 ylo = bcoeff - 1000000000ull * yhi;
200 if (ylo >= 1000000000) {
201 ylo = ylo - 1000000000;
202 yhi = yhi + 1;
203 }
204 d103 = 0x4189374c;
205 B34 = ((UINT64) ylo * d103) >> (32 + 8);
206 B01 = ((UINT64) yhi * d103) >> (32 + 8);
207 b5 = ylo - B34 * 1000;
208 b2 = yhi - B01 * 1000;
209 b3 = ((UINT64) B34 * d103) >> (32 + 8);
210 b0 = ((UINT64) B01 * d103) >> (32 + 8);
211 b4 = (unsigned int) B34 - (unsigned int) b3 *1000;
212 b1 = (unsigned int) B01 - (unsigned int) dm103[b0];
213 dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1];
214 if (b0 >= 8) /* is b0 8 or 9? */
215 res = sign | ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) |
216 (exp & 0xff)) << 50) | dcoeff;
217 else /* else b0 is 0..7 */
218 res = sign | ((((exp >> 8) << 11) | (b0 << 8) |
219 (exp & 0xff)) << 50) | dcoeff;
220 }
221 *pres = res;
222 }
223
224 void _dpd_to_bid64 (_Decimal64 *, _Decimal64 *);
225
226 void
227 _dpd_to_bid64 (_Decimal64 *pres, _Decimal64 *px) {
228 UINT64 res;
229 UINT64 sign, comb, exp;
230 UINT64 trailing;
231 UINT64 d0, d1, d2;
232 unsigned int d3, d4, d5;
233 UINT64 bcoeff, mask;
234 _Decimal64 x = *px;
235
236 sign = (x & 0x8000000000000000ull);
237 comb = (x & 0x7ffc000000000000ull) >> 50;
238 trailing = (x & 0x0003ffffffffffffull);
239 if ((comb & 0x1e00) == 0x1e00) {
240 if ((comb & 0x1f00) == 0x1f00) { /* G0..G4 = 11111 -> NaN */
241 if (comb & 0x0100) { /* G5 = 1 -> sNaN */
242 *pres = x;
243 } else { /* G5 = 0 -> qNaN */
244 *pres = x;
245 }
246 } else { /*if ((comb & 0x1e00) == 0x1e00); G0..G4 = 11110 -> INF */
247 *pres = x;
248 }
249 return;
250 } else { /* normal number */
251 if ((comb & 0x1800) == 0x1800) { /* G0..G1 = 11 -> d0 = 8 + G4 */
252 d0 = d2b6[((comb >> 8) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
253 exp = (comb & 0x600) >> 1; /* exp = (comb & 0x0400 ? 1 : 0) * 0x200 +
254 (comb & 0x0200 ? 1 : 0) * 0x100; exp leading bits are G2..G3 */
255 } else {
256 d0 = d2b6[(comb >> 8) & 0x7];
257 exp = (comb & 0x1800) >> 3; /* exp = (comb & 0x1000 ? 1 : 0) * 0x200 +
258 (comb & 0x0800 ? 1 : 0) * 0x100; exp loading bits are G0..G1 */
259 }
260 d1 = d2b5[(trailing >> 40) & 0x3ff];
261 d2 = d2b4[(trailing >> 30) & 0x3ff];
262 d3 = d2b3[(trailing >> 20) & 0x3ff];
263 d4 = d2b2[(trailing >> 10) & 0x3ff];
264 d5 = d2b[(trailing) & 0x3ff];
265 bcoeff = (d5 + d4 + d3) + d2 + d1 + d0;
266 exp += (comb & 0xff);
267 mask = 1;
268 mask <<= 53;
269 if (bcoeff < mask) { /* check whether coefficient fits in 10*5+3 bits */
270 res = exp;
271 res <<= 53;
272 res |= (bcoeff | sign);
273 *pres = res;
274 return;
275 }
276 /* special format */
277 res = (exp << 51) | (sign | 0x6000000000000000ull);
278 /* add coeff, without leading bits */
279 mask = (mask >> 2) - 1;
280 bcoeff &= mask;
281 res |= bcoeff;
282 }
283 *pres = res;
284 }
285
286 void _bid_to_dpd128 (_Decimal128 *, _Decimal128 *);
287
288 void
289 _bid_to_dpd128 (_Decimal128 *pres, _Decimal128 *px) {
290 UINT128 res;
291 UINT128 sign;
292 unsigned int comb;
293 UINT128 bcoeff;
294 UINT128 dcoeff;
295 UINT128 BH, d1018, BT2, BT1;
296 UINT64 exp, BL, d109;
297 UINT64 d106, d103;
298 UINT64 k1, k2, k4, k5, k7, k8, k10, k11;
299 unsigned int BHH32, BLL32, BHL32, BLH32, k0, k3, k6, k9, amount;
300 _Decimal128 x = *px;
301
302 sign.w[1] = (x.w[1] & 0x8000000000000000ull);
303 sign.w[0] = 0;
304 comb = (x.w[1] /*& 0x7fffc00000000000ull */ ) >> 46;
305 exp = 0;
306 if ((comb & 0x1e000) == 0x1e000) {
307 if ((comb & 0x1f000) == 0x1f000) { /* G0..G4 = 11111 -> NaN */
308 if (comb & 0x01000) { /* G5 = 1 -> sNaN */
309 res = x;
310 } else { /* G5 = 0 -> qNaN */
311 res = x;
312 }
313 } else { /* G0..G4 = 11110 -> INF */
314 res = x;
315 }
316 } else { /* normal number */
317 exp = ((x.w[1] & 0x7fff000000000000ull) >> 49) & 0x3fff;
318 bcoeff.w[1] = (x.w[1] & 0x0001ffffffffffffull);
319 bcoeff.w[0] = x.w[0];
320 d1018 = reciprocals10_128[18];
321 __mul_128x128_high (BH, bcoeff, d1018);
322 amount = recip_scale[18];
323 BH.w[0] = (BH.w[0] >> amount) | (BH.w[1] << (64 - amount));
324 BL = bcoeff.w[0] - BH.w[0] * 1000000000000000000ull;
325 d109 = reciprocals10_64[9];
326 __mul_64x64_to_128 (BT1, BH.w[0], d109);
327 BHH32 = (unsigned int) (BT1.w[1] >> short_recip_scale[9]);
328 BHL32 = (unsigned int) BH.w[0] - BHH32 * 1000000000;
329 __mul_64x64_to_128 (BT2, BL, d109);
330 BLH32 = (unsigned int) (BT2.w[1] >> short_recip_scale[9]);
331 BLL32 = (unsigned int) BL - BLH32 * 1000000000;
332 d106 = 0x431BDE83;
333 d103 = 0x4189374c;
334 k0 = ((UINT64) BHH32 * d106) >> (32 + 18);
335 BHH32 -= (unsigned int) k0 *1000000;
336 k1 = ((UINT64) BHH32 * d103) >> (32 + 8);
337 k2 = BHH32 - (unsigned int) k1 *1000;
338 k3 = ((UINT64) BHL32 * d106) >> (32 + 18);
339 BHL32 -= (unsigned int) k3 *1000000;
340 k4 = ((UINT64) BHL32 * d103) >> (32 + 8);
341 k5 = BHL32 - (unsigned int) k4 *1000;
342 k6 = ((UINT64) BLH32 * d106) >> (32 + 18);
343 BLH32 -= (unsigned int) k6 *1000000;
344 k7 = ((UINT64) BLH32 * d103) >> (32 + 8);
345 k8 = BLH32 - (unsigned int) k7 *1000;
346 k9 = ((UINT64) BLL32 * d106) >> (32 + 18);
347 BLL32 -= (unsigned int) k9 *1000000;
348 k10 = ((UINT64) BLL32 * d103) >> (32 + 8);
349 k11 = BLL32 - (unsigned int) k10 *1000;
350 dcoeff.w[1] = (b2d[k5] >> 4) | (b2d[k4] << 6) | (b2d[k3] << 16) |
351 (b2d[k2] << 26) | (b2d[k1] << 36);
352 dcoeff.w[0] = b2d[k11] | (b2d[k10] << 10) | (b2d[k9] << 20) |
353 (b2d[k8] << 30) | (b2d[k7] << 40) | (b2d[k6] << 50) | (b2d[k5] << 60);
354 res.w[0] = dcoeff.w[0];
355 if (k0 >= 8) {
356 res.w[1] = sign.w[1] | ((0x18000 | ((exp >> 12) << 13) |
357 ((k0 & 1) << 12) | (exp & 0xfff)) << 46) | dcoeff.w[1];
358 } else {
359 res.w[1] = sign.w[1] | ((((exp >> 12) << 15) | (k0 << 12) |
360 (exp & 0xfff)) << 46) | dcoeff.w[1];
361 }
362 }
363 *pres = res;
364 }
365
366 void _dpd_to_bid128 (_Decimal128 *, _Decimal128 *);
367
368 void
369 _dpd_to_bid128 (_Decimal128 *pres, _Decimal128 *px) {
370 UINT128 res;
371 UINT128 sign;
372 UINT64 exp, comb;
373 UINT128 trailing;
374 UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
375 UINT128 bcoeff;
376 UINT64 tl, th;
377 _Decimal128 x = *px;
378
379 sign.w[1] = (x.w[1] & 0x8000000000000000ull);
380 sign.w[0] = 0;
381 comb = (x.w[1] & 0x7fffc00000000000ull) >> 46;
382 trailing.w[1] = x.w[1];
383 trailing.w[0] = x.w[0];
384 if ((comb & 0x1e000) == 0x1e000) {
385 if ((comb & 0x1f000) == 0x1f000) { /* G0..G4 = 11111 -> NaN */
386 if (comb & 0x01000) { /* G5 = 1 -> sNaN */
387 *pres = x;
388 } else { /* G5 = 0 -> qNaN */
389 *pres = x;
390 }
391 } else { /* G0..G4 = 11110 -> INF */
392 *pres = x;
393 }
394 return;
395 } else { /* Normal number */
396 if ((comb & 0x18000) == 0x18000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
397 d0 = d2b6[8 + ((comb & 0x01000) >> 12)];
398 exp = (comb & 0x06000) >> 1; /* exp leading bits are G2..G3 */
399 } else {
400 d0 = d2b6[((comb & 0x07000) >> 12)];
401 exp = (comb & 0x18000) >> 3; /* exp loading bits are G0..G1 */
402 }
403 d11 = d2b[(trailing.w[0]) & 0x3ff];
404 d10 = d2b2[(trailing.w[0] >> 10) & 0x3ff];
405 d9 = d2b3[(trailing.w[0] >> 20) & 0x3ff];
406 d8 = d2b4[(trailing.w[0] >> 30) & 0x3ff];
407 d7 = d2b5[(trailing.w[0] >> 40) & 0x3ff];
408 d6 = d2b6[(trailing.w[0] >> 50) & 0x3ff];
409 d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)];
410 d4 = d2b2[(trailing.w[1] >> 6) & 0x3ff];
411 d3 = d2b3[(trailing.w[1] >> 16) & 0x3ff];
412 d2 = d2b4[(trailing.w[1] >> 26) & 0x3ff];
413 d1 = d2b5[(trailing.w[1] >> 36) & 0x3ff];
414 tl = d11 + d10 + d9 + d8 + d7 + d6;
415 th = d5 + d4 + d3 + d2 + d1 + d0;
416 __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull);
417 __add_128_64 (bcoeff, bcoeff, tl);
418 exp += (comb & 0xfff);
419 res.w[0] = bcoeff.w[0];
420 res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1];
421 }
422 *pres = res;
423 }