1 /* Copyright (C) 2007-2023 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 }
142 /* normal number */
143 if ((x & 0x60000000) == 0x60000000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
144 d0 = d2b3[((x >> 26) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
145 exp = (x >> 27) & 3; /* exp leading bits are G2..G3 */
146 } else {
147 d0 = d2b3[(x >> 26) & 0x7];
148 exp = (x >> 29) & 3; /* exp loading bits are G0..G1 */
149 }
150 d1 = d2b2[(trailing >> 10) & 0x3ff];
151 d2 = d2b[(trailing) & 0x3ff];
152 bcoeff = d2 + d1 + d0;
153 exp = (exp << 6) + ((x >> 20) & 0x3f);
154 if (bcoeff < (1 << 23)) {
155 r = exp;
156 r <<= 23;
157 r |= (bcoeff | sign);
158 } else {
159 r = exp;
160 r <<= 21;
161 r |= (sign | 0x60000000ul);
162 /* add coeff, without leading bits */
163 r |= (((unsigned int) bcoeff) & 0x1fffff);
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 }
188 /* Normal number */
189 if ((comb & 0xc00) == 0xc00) { /* G0..G1 = 11 -> exp is G2..G11 */
190 exp = (comb) & 0x3ff;
191 bcoeff = (x & 0x0007ffffffffffffull) | 0x0020000000000000ull;
192 if (bcoeff >= 10000000000000000ull)
193 bcoeff = 0;
194 } else {
195 exp = (comb >> 2) & 0x3ff;
196 bcoeff = (x & 0x001fffffffffffffull);
197 }
198 D61 = 2305843009ull; /* Floor(2^61 / 10^9) */
199 /* Multiply the binary coefficient by ceil(2^64 / 1000), and take the upper
200 64-bits in order to compute a division by 1000. */
201 yhi = (D61 * (UINT64)(bcoeff >> (UINT64)27)) >> (UINT64)34;
202 ylo = bcoeff - 1000000000ull * yhi;
203 if (ylo >= 1000000000) {
204 ylo = ylo - 1000000000;
205 yhi = yhi + 1;
206 }
207 d103 = 0x4189374c;
208 B34 = ((UINT64) ylo * d103) >> (32 + 8);
209 B01 = ((UINT64) yhi * d103) >> (32 + 8);
210 b5 = ylo - B34 * 1000;
211 b2 = yhi - B01 * 1000;
212 b3 = ((UINT64) B34 * d103) >> (32 + 8);
213 b0 = ((UINT64) B01 * d103) >> (32 + 8);
214 b4 = (unsigned int) B34 - (unsigned int) b3 *1000;
215 b1 = (unsigned int) B01 - (unsigned int) dm103[b0];
216 dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1];
217 if (b0 >= 8) /* is b0 8 or 9? */
218 res = sign | ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) |
219 (exp & 0xff)) << 50) | dcoeff;
220 else /* else b0 is 0..7 */
221 res = sign | ((((exp >> 8) << 11) | (b0 << 8) |
222 (exp & 0xff)) << 50) | dcoeff;
223 *pres = res;
224 }
225
226 void _dpd_to_bid64 (_Decimal64 *, _Decimal64 *);
227
228 void
229 _dpd_to_bid64 (_Decimal64 *pres, _Decimal64 *px) {
230 UINT64 res;
231 UINT64 sign, comb, exp;
232 UINT64 trailing;
233 UINT64 d0, d1, d2;
234 unsigned int d3, d4, d5;
235 UINT64 bcoeff, mask;
236 _Decimal64 x = *px;
237
238 sign = (x & 0x8000000000000000ull);
239 comb = (x & 0x7ffc000000000000ull) >> 50;
240 trailing = (x & 0x0003ffffffffffffull);
241 if ((comb & 0x1e00) == 0x1e00) {
242 *pres = x;
243 return;
244 }
245 /* normal number */
246 if ((comb & 0x1800) == 0x1800) { /* G0..G1 = 11 -> d0 = 8 + G4 */
247 d0 = d2b6[((comb >> 8) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */
248 exp = (comb & 0x600) >> 1; /* exp = (comb & 0x0400 ? 1 : 0) * 0x200 +
249 (comb & 0x0200 ? 1 : 0) * 0x100; exp leading bits are G2..G3 */
250 } else {
251 d0 = d2b6[(comb >> 8) & 0x7];
252 exp = (comb & 0x1800) >> 3; /* exp = (comb & 0x1000 ? 1 : 0) * 0x200 +
253 (comb & 0x0800 ? 1 : 0) * 0x100; exp loading bits are G0..G1 */
254 }
255 d1 = d2b5[(trailing >> 40) & 0x3ff];
256 d2 = d2b4[(trailing >> 30) & 0x3ff];
257 d3 = d2b3[(trailing >> 20) & 0x3ff];
258 d4 = d2b2[(trailing >> 10) & 0x3ff];
259 d5 = d2b[(trailing) & 0x3ff];
260 bcoeff = (d5 + d4 + d3) + d2 + d1 + d0;
261 exp += (comb & 0xff);
262 mask = 1;
263 mask <<= 53;
264 if (bcoeff < mask) { /* check whether coefficient fits in 10*5+3 bits */
265 res = exp;
266 res <<= 53;
267 res |= (bcoeff | sign);
268 *pres = res;
269 return;
270 }
271 /* special format */
272 res = (exp << 51) | (sign | 0x6000000000000000ull);
273 /* add coeff, without leading bits */
274 mask = (mask >> 2) - 1;
275 bcoeff &= mask;
276 res |= bcoeff;
277 *pres = res;
278 }
279
280 void _bid_to_dpd128 (_Decimal128 *, _Decimal128 *);
281
282 void
283 _bid_to_dpd128 (_Decimal128 *pres, _Decimal128 *px) {
284 UINT128 res;
285 UINT128 sign;
286 unsigned int comb;
287 UINT128 bcoeff;
288 UINT128 dcoeff;
289 UINT128 BH, d1018, BT2, BT1;
290 UINT64 exp, BL, d109;
291 UINT64 d106, d103;
292 UINT64 k1, k2, k4, k5, k7, k8, k10, k11;
293 unsigned int BHH32, BLL32, BHL32, BLH32, k0, k3, k6, k9, amount;
294 _Decimal128 x = *px;
295
296 sign.w[1] = (x.w[1] & 0x8000000000000000ull);
297 sign.w[0] = 0;
298 comb = (x.w[1] /*& 0x7fffc00000000000ull */ ) >> 46;
299 exp = 0;
300 if ((comb & 0x1e000) == 0x1e000) {
301 res = x;
302 } else { /* normal number */
303 if ((comb & 0x18000) == 0x18000) {
304 /* Noncanonical significand (prepending 8 or 9 to any 110-bit
305 trailing significand field produces a value above 10^34). */
306 exp = (comb & 0x7fff) >> 1;
307 bcoeff.w[1] = 0;
308 bcoeff.w[0] = 0;
309 } else {
310 exp = ((x.w[1] & 0x7fff000000000000ull) >> 49) & 0x3fff;
311 bcoeff.w[1] = (x.w[1] & 0x0001ffffffffffffull);
312 bcoeff.w[0] = x.w[0];
313 if (bcoeff.w[1] > 0x1ed09bead87c0ull
314 || (bcoeff.w[1] == 0x1ed09bead87c0ull
315 && bcoeff.w[0] >= 0x378d8e6400000000ull)) {
316 bcoeff.w[1] = 0;
317 bcoeff.w[0] = 0;
318 }
319 }
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 *pres = x;
386 return;
387 }
388 if ((comb & 0x18000) == 0x18000) { /* G0..G1 = 11 -> d0 = 8 + G4 */
389 d0 = d2b6[8 + ((comb & 0x01000) >> 12)];
390 exp = (comb & 0x06000) >> 1; /* exp leading bits are G2..G3 */
391 } else {
392 d0 = d2b6[((comb & 0x07000) >> 12)];
393 exp = (comb & 0x18000) >> 3; /* exp loading bits are G0..G1 */
394 }
395 d11 = d2b[(trailing.w[0]) & 0x3ff];
396 d10 = d2b2[(trailing.w[0] >> 10) & 0x3ff];
397 d9 = d2b3[(trailing.w[0] >> 20) & 0x3ff];
398 d8 = d2b4[(trailing.w[0] >> 30) & 0x3ff];
399 d7 = d2b5[(trailing.w[0] >> 40) & 0x3ff];
400 d6 = d2b6[(trailing.w[0] >> 50) & 0x3ff];
401 d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)];
402 d4 = d2b2[(trailing.w[1] >> 6) & 0x3ff];
403 d3 = d2b3[(trailing.w[1] >> 16) & 0x3ff];
404 d2 = d2b4[(trailing.w[1] >> 26) & 0x3ff];
405 d1 = d2b5[(trailing.w[1] >> 36) & 0x3ff];
406 tl = d11 + d10 + d9 + d8 + d7 + d6;
407 th = d5 + d4 + d3 + d2 + d1 + d0;
408 __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull);
409 __add_128_64 (bcoeff, bcoeff, tl);
410 exp += (comb & 0xfff);
411 res.w[0] = bcoeff.w[0];
412 res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1];
413 *pres = res;
414 }