diff options
Diffstat (limited to 'libdecnumber/bid/bid2dpd_dpd2bid.c')
| -rw-r--r-- | libdecnumber/bid/bid2dpd_dpd2bid.c | 428 |
1 files changed, 428 insertions, 0 deletions
diff --git a/libdecnumber/bid/bid2dpd_dpd2bid.c b/libdecnumber/bid/bid2dpd_dpd2bid.c new file mode 100644 index 0000000..01e83bc --- /dev/null +++ b/libdecnumber/bid/bid2dpd_dpd2bid.c @@ -0,0 +1,428 @@ +/* Copyright (C) 2007 Free Software Foundation, Inc. + +This file is part of GCC. + +GCC is free software; you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free +Software Foundation; either version 2, or (at your option) any later +version. + +In addition to the permissions in the GNU General Public License, the +Free Software Foundation gives you unlimited permission to link the +compiled version of this file into combinations with other programs, +and to distribute those combinations without any restriction coming +from the use of this file. (The General Public License restrictions +do apply in other respects; for example, they cover modification of +the file, and distribution when not linked into a combine +executable.) + +GCC is distributed in the hope that it will be useful, but WITHOUT ANY +WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with GCC; see the file COPYING. If not, write to the Free +Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. */ + +#undef IN_LIBGCC2 +#include "bid-dpd.h" + +/* get full 64x64bit product */ +#define __mul_64x64_to_128(P, CX, CY) \ +{ \ + UINT64 CXH, CXL, CYH,CYL,PL,PH,PM,PM2; \ + CXH = (CX) >> 32; \ + CXL = (UINT32)(CX); \ + CYH = (CY) >> 32; \ + CYL = (UINT32)(CY); \ + \ + PM = CXH*CYL; \ + PH = CXH*CYH; \ + PL = CXL*CYL; \ + PM2 = CXL*CYH; \ + PH += (PM>>32); \ + PM = (UINT64)((UINT32)PM)+PM2+(PL>>32); \ + \ + (P).w[1] = PH + (PM>>32); \ + (P).w[0] = (PM<<32)+(UINT32)PL; \ +} + +/* add 64-bit value to 128-bit */ +#define __add_128_64(R128, A128, B64) \ +{ \ + UINT64 R64H; \ + R64H = (A128).w[1]; \ + (R128).w[0] = (B64) + (A128).w[0]; \ + if((R128).w[0] < (B64)) R64H ++; \ + (R128).w[1] = R64H; \ +} + +/* add 128-bit value to 128-bit (assume no carry-out) */ +#define __add_128_128(R128, A128, B128) \ +{ \ + UINT128 Q128; \ + Q128.w[1] = (A128).w[1]+(B128).w[1]; \ + Q128.w[0] = (B128).w[0] + (A128).w[0]; \ + if(Q128.w[0] < (B128).w[0]) Q128.w[1] ++; \ + (R128).w[1] = Q128.w[1]; \ + (R128).w[0] = Q128.w[0]; \ +} + +#define __mul_128x128_high(Q, A, B) \ +{ \ + UINT128 ALBL, ALBH, AHBL, AHBH, QM, QM2; \ + \ + __mul_64x64_to_128(ALBH, (A).w[0], (B).w[1]); \ + __mul_64x64_to_128(AHBL, (B).w[0], (A).w[1]); \ + __mul_64x64_to_128(ALBL, (A).w[0], (B).w[0]); \ + __mul_64x64_to_128(AHBH, (A).w[1],(B).w[1]); \ + \ + __add_128_128(QM, ALBH, AHBL); \ + __add_128_64(QM2, QM, ALBL.w[1]); \ + __add_128_64((Q), AHBH, QM2.w[1]); \ +} + +#include "bid2dpd_dpd2bid.h" + +static const unsigned int dm103[] = + { 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 11000 }; + +void _bid_to_dpd32 (_Decimal32 *, _Decimal32 *); + +void +_bid_to_dpd32 (_Decimal32 *pres, _Decimal32 *px) { + unsigned int sign, coefficient_x, exp, dcoeff; + unsigned int b2, b1, b0, b01, res; + _Decimal32 x = *px; + + sign = (x & 0x80000000); + if ((x & 0x60000000ul) == 0x60000000ul) { + /* special encodings */ + if ((x & 0x78000000ul) == 0x78000000ul) { + *pres = x; /* NaN or Infinity */ + return; + } + /* coefficient */ + coefficient_x = (x & 0x001ffffful) | 0x00800000ul; + if (coefficient_x >= 10000000) coefficient_x = 0; + /* get exponent */ + exp = (x >> 21) & 0xff; + } else { + exp = (x >> 23) & 0xff; + coefficient_x = (x & 0x007ffffful); + } + b01 = coefficient_x / 1000; + b2 = coefficient_x - 1000 * b01; + b0 = b01 / 1000; + b1 = b01 - 1000 * b0; + dcoeff = b2d[b2] | b2d2[b1]; + if (b0 >= 8) { /* is b0 8 or 9? */ + res = sign | ((0x600 | ((exp >> 6) << 7) | + ((b0 & 1) << 6) | (exp & 0x3f)) << 20) | dcoeff; + } else { /* else b0 is 0..7 */ + res = sign | ((((exp >> 6) << 9) | (b0 << 6) | + (exp & 0x3f)) << 20) | dcoeff; + } + *pres = res; +} + +void _dpd_to_bid32 (_Decimal32 *, _Decimal32 *); + +void +_dpd_to_bid32 (_Decimal32 *pres, _Decimal32 *px) { + unsigned int r; + unsigned int sign, exp, bcoeff; + UINT64 trailing; + unsigned int d0, d1, d2; + _Decimal32 x = *px; + + sign = (x & 0x80000000); + trailing = (x & 0x000fffff); + if ((x & 0x78000000) == 0x78000000) { + *pres = x; + return; + } else { /* normal number */ + if ((x & 0x60000000) == 0x60000000) { /* G0..G1 = 11 -> d0 = 8 + G4 */ + d0 = d2b3[((x >> 26) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */ + exp = (x >> 27) & 3; /* exp leading bits are G2..G3 */ + } else { + d0 = d2b3[(x >> 26) & 0x7]; + exp = (x >> 29) & 3; /* exp loading bits are G0..G1 */ + } + d1 = d2b2[(trailing >> 10) & 0x3ff]; + d2 = d2b[(trailing) & 0x3ff]; + bcoeff = d2 + d1 + d0; + exp = (exp << 6) + ((x >> 20) & 0x3f); + if (bcoeff < (1 << 23)) { + r = exp; + r <<= 23; + r |= (bcoeff | sign); + } else { + r = exp; + r <<= 21; + r |= (sign | 0x60000000ul); + /* add coeff, without leading bits */ + r |= (((unsigned int) bcoeff) & 0x1fffff); + } + } + *pres = r; +} + +void _bid_to_dpd64 (_Decimal64 *, _Decimal64 *); + +void +_bid_to_dpd64 (_Decimal64 *pres, _Decimal64 *px) { + UINT64 res; + UINT64 sign, comb, exp, B34, B01; + UINT64 d103, D61; + UINT64 b0, b2, b3, b5; + unsigned int b1, b4; + UINT64 bcoeff; + UINT64 dcoeff; + unsigned int yhi, ylo; + _Decimal64 x = *px; + + sign = (x & 0x8000000000000000ull); + comb = (x & 0x7ffc000000000000ull) >> 51; + if ((comb & 0xf00) == 0xf00) { + *pres = x; + return; + } else { /* Normal number */ + if ((comb & 0xc00) == 0xc00) { /* G0..G1 = 11 -> exp is G2..G11 */ + exp = (comb) & 0x3ff; + bcoeff = (x & 0x0007ffffffffffffull) | 0x0020000000000000ull; + } else { + exp = (comb >> 2) & 0x3ff; + bcoeff = (x & 0x001fffffffffffffull); + } + D61 = 2305843009ull; /* Floor(2^61 / 10^9) */ + /* Multiply the binary coefficient by ceil(2^64 / 1000), and take the upper + 64-bits in order to compute a division by 1000. */ + yhi = (D61 * (UINT64)(bcoeff >> (UINT64)27)) >> (UINT64)34; + ylo = bcoeff - 1000000000ull * yhi; + if (ylo >= 1000000000) { + ylo = ylo - 1000000000; + yhi = yhi + 1; + } + d103 = 0x4189374c; + B34 = ((UINT64) ylo * d103) >> (32 + 8); + B01 = ((UINT64) yhi * d103) >> (32 + 8); + b5 = ylo - B34 * 1000; + b2 = yhi - B01 * 1000; + b3 = ((UINT64) B34 * d103) >> (32 + 8); + b0 = ((UINT64) B01 * d103) >> (32 + 8); + b4 = (unsigned int) B34 - (unsigned int) b3 *1000; + b1 = (unsigned int) B01 - (unsigned int) dm103[b0]; + dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1]; + if (b0 >= 8) /* is b0 8 or 9? */ + res = sign | ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) | + (exp & 0xff)) << 50) | dcoeff; + else /* else b0 is 0..7 */ + res = sign | ((((exp >> 8) << 11) | (b0 << 8) | + (exp & 0xff)) << 50) | dcoeff; + } + *pres = res; +} + +void _dpd_to_bid64 (_Decimal64 *, _Decimal64 *); + +void +_dpd_to_bid64 (_Decimal64 *pres, _Decimal64 *px) { + UINT64 res; + UINT64 sign, comb, exp; + UINT64 trailing; + UINT64 d0, d1, d2; + unsigned int d3, d4, d5; + UINT64 bcoeff, mask; + _Decimal64 x = *px; + + sign = (x & 0x8000000000000000ull); + comb = (x & 0x7ffc000000000000ull) >> 50; + trailing = (x & 0x0003ffffffffffffull); + if ((comb & 0x1e00) == 0x1e00) { + if ((comb & 0x1f00) == 0x1f00) { /* G0..G4 = 11111 -> NaN */ + if (comb & 0x0100) { /* G5 = 1 -> sNaN */ + *pres = x; + } else { /* G5 = 0 -> qNaN */ + *pres = x; + } + } else { /*if ((comb & 0x1e00) == 0x1e00); G0..G4 = 11110 -> INF */ + *pres = x; + } + return; + } else { /* normal number */ + if ((comb & 0x1800) == 0x1800) { /* G0..G1 = 11 -> d0 = 8 + G4 */ + d0 = d2b6[((comb >> 8) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */ + exp = (comb & 0x600) >> 1; /* exp = (comb & 0x0400 ? 1 : 0) * 0x200 + + (comb & 0x0200 ? 1 : 0) * 0x100; exp leading bits are G2..G3 */ + } else { + d0 = d2b6[(comb >> 8) & 0x7]; + exp = (comb & 0x1800) >> 3; /* exp = (comb & 0x1000 ? 1 : 0) * 0x200 + + (comb & 0x0800 ? 1 : 0) * 0x100; exp loading bits are G0..G1 */ + } + d1 = d2b5[(trailing >> 40) & 0x3ff]; + d2 = d2b4[(trailing >> 30) & 0x3ff]; + d3 = d2b3[(trailing >> 20) & 0x3ff]; + d4 = d2b2[(trailing >> 10) & 0x3ff]; + d5 = d2b[(trailing) & 0x3ff]; + bcoeff = (d5 + d4 + d3) + d2 + d1 + d0; + exp += (comb & 0xff); + mask = 1; + mask <<= 53; + if (bcoeff < mask) { /* check whether coefficient fits in 10*5+3 bits */ + res = exp; + res <<= 53; + res |= (bcoeff | sign); + *pres = res; + return; + } + /* special format */ + res = (exp << 51) | (sign | 0x6000000000000000ull); + /* add coeff, without leading bits */ + mask = (mask >> 2) - 1; + bcoeff &= mask; + res |= bcoeff; + } + *pres = res; +} + +void _bid_to_dpd128 (_Decimal128 *, _Decimal128 *); + +void +_bid_to_dpd128 (_Decimal128 *pres, _Decimal128 *px) { + UINT128 res; + UINT128 sign; + unsigned int comb; + UINT128 bcoeff; + UINT128 dcoeff; + UINT128 BH, d1018, BT2, BT1; + UINT64 exp, BL, d109; + UINT64 d106, d103; + UINT64 k1, k2, k4, k5, k7, k8, k10, k11; + unsigned int BHH32, BLL32, BHL32, BLH32, k0, k3, k6, k9, amount; + _Decimal128 x = *px; + + sign.w[1] = (x.w[1] & 0x8000000000000000ull); + sign.w[0] = 0; + comb = (x.w[1] /*& 0x7fffc00000000000ull */ ) >> 46; + exp = 0; + if ((comb & 0x1e000) == 0x1e000) { + if ((comb & 0x1f000) == 0x1f000) { /* G0..G4 = 11111 -> NaN */ + if (comb & 0x01000) { /* G5 = 1 -> sNaN */ + res = x; + } else { /* G5 = 0 -> qNaN */ + res = x; + } + } else { /* G0..G4 = 11110 -> INF */ + res = x; + } + } else { /* normal number */ + exp = ((x.w[1] & 0x7fff000000000000ull) >> 49) & 0x3fff; + bcoeff.w[1] = (x.w[1] & 0x0001ffffffffffffull); + bcoeff.w[0] = x.w[0]; + d1018 = reciprocals10_128[18]; + __mul_128x128_high (BH, bcoeff, d1018); + amount = recip_scale[18]; + BH.w[0] = (BH.w[0] >> amount) | (BH.w[1] << (64 - amount)); + BL = bcoeff.w[0] - BH.w[0] * 1000000000000000000ull; + d109 = reciprocals10_64[9]; + __mul_64x64_to_128 (BT1, BH.w[0], d109); + BHH32 = (unsigned int) (BT1.w[1] >> short_recip_scale[9]); + BHL32 = (unsigned int) BH.w[0] - BHH32 * 1000000000; + __mul_64x64_to_128 (BT2, BL, d109); + BLH32 = (unsigned int) (BT2.w[1] >> short_recip_scale[9]); + BLL32 = (unsigned int) BL - BLH32 * 1000000000; + d106 = 0x431BDE83; + d103 = 0x4189374c; + k0 = ((UINT64) BHH32 * d106) >> (32 + 18); + BHH32 -= (unsigned int) k0 *1000000; + k1 = ((UINT64) BHH32 * d103) >> (32 + 8); + k2 = BHH32 - (unsigned int) k1 *1000; + k3 = ((UINT64) BHL32 * d106) >> (32 + 18); + BHL32 -= (unsigned int) k3 *1000000; + k4 = ((UINT64) BHL32 * d103) >> (32 + 8); + k5 = BHL32 - (unsigned int) k4 *1000; + k6 = ((UINT64) BLH32 * d106) >> (32 + 18); + BLH32 -= (unsigned int) k6 *1000000; + k7 = ((UINT64) BLH32 * d103) >> (32 + 8); + k8 = BLH32 - (unsigned int) k7 *1000; + k9 = ((UINT64) BLL32 * d106) >> (32 + 18); + BLL32 -= (unsigned int) k9 *1000000; + k10 = ((UINT64) BLL32 * d103) >> (32 + 8); + k11 = BLL32 - (unsigned int) k10 *1000; + dcoeff.w[1] = (b2d[k5] >> 4) | (b2d[k4] << 6) | (b2d[k3] << 16) | + (b2d[k2] << 26) | (b2d[k1] << 36); + dcoeff.w[0] = b2d[k11] | (b2d[k10] << 10) | (b2d[k9] << 20) | + (b2d[k8] << 30) | (b2d[k7] << 40) | (b2d[k6] << 50) | (b2d[k5] << 60); + res.w[0] = dcoeff.w[0]; + if (k0 >= 8) { + res.w[1] = sign.w[1] | ((0x18000 | ((exp >> 12) << 13) | + ((k0 & 1) << 12) | (exp & 0xfff)) << 46) | dcoeff.w[1]; + } else { + res.w[1] = sign.w[1] | ((((exp >> 12) << 15) | (k0 << 12) | + (exp & 0xfff)) << 46) | dcoeff.w[1]; + } + } + *pres = res; +} + +void _dpd_to_bid128 (_Decimal128 *, _Decimal128 *); + +void +_dpd_to_bid128 (_Decimal128 *pres, _Decimal128 *px) { + UINT128 res; + UINT128 sign; + UINT64 exp, comb; + UINT128 trailing; + UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11; + UINT128 bcoeff; + UINT64 tl, th; + _Decimal128 x = *px; + + sign.w[1] = (x.w[1] & 0x8000000000000000ull); + sign.w[0] = 0; + comb = (x.w[1] & 0x7fffc00000000000ull) >> 46; + trailing.w[1] = x.w[1]; + trailing.w[0] = x.w[0]; + if ((comb & 0x1e000) == 0x1e000) { + if ((comb & 0x1f000) == 0x1f000) { /* G0..G4 = 11111 -> NaN */ + if (comb & 0x01000) { /* G5 = 1 -> sNaN */ + *pres = x; + } else { /* G5 = 0 -> qNaN */ + *pres = x; + } + } else { /* G0..G4 = 11110 -> INF */ + *pres = x; + } + return; + } else { /* Normal number */ + if ((comb & 0x18000) == 0x18000) { /* G0..G1 = 11 -> d0 = 8 + G4 */ + d0 = d2b6[8 + ((comb & 0x01000) >> 12)]; + exp = (comb & 0x06000) >> 1; /* exp leading bits are G2..G3 */ + } else { + d0 = d2b6[((comb & 0x07000) >> 12)]; + exp = (comb & 0x18000) >> 3; /* exp loading bits are G0..G1 */ + } + d11 = d2b[(trailing.w[0]) & 0x3ff]; + d10 = d2b2[(trailing.w[0] >> 10) & 0x3ff]; + d9 = d2b3[(trailing.w[0] >> 20) & 0x3ff]; + d8 = d2b4[(trailing.w[0] >> 30) & 0x3ff]; + d7 = d2b5[(trailing.w[0] >> 40) & 0x3ff]; + d6 = d2b6[(trailing.w[0] >> 50) & 0x3ff]; + d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)]; + d4 = d2b2[(trailing.w[1] >> 6) & 0x3ff]; + d3 = d2b3[(trailing.w[1] >> 16) & 0x3ff]; + d2 = d2b4[(trailing.w[1] >> 26) & 0x3ff]; + d1 = d2b5[(trailing.w[1] >> 36) & 0x3ff]; + tl = d11 + d10 + d9 + d8 + d7 + d6; + th = d5 + d4 + d3 + d2 + d1 + d0; + __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull); + __add_128_64 (bcoeff, bcoeff, tl); + exp += (comb & 0xfff); + res.w[0] = bcoeff.w[0]; + res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1]; + } + *pres = res; +} |