0
|
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 }
|