| 1 | /* Copyright (C) 2007-2018 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 | } else { |
| 193 | exp = (comb >> 2) & 0x3ff; |
| 194 | bcoeff = (x & 0x001fffffffffffffull); |
| 195 | } |
| 196 | D61 = 2305843009ull; /* Floor(2^61 / 10^9) */ |
| 197 | /* Multiply the binary coefficient by ceil(2^64 / 1000), and take the upper |
| 198 | 64-bits in order to compute a division by 1000. */ |
| 199 | yhi = (D61 * (UINT64)(bcoeff >> (UINT64)27)) >> (UINT64)34; |
| 200 | ylo = bcoeff - 1000000000ull * yhi; |
| 201 | if (ylo >= 1000000000) { |
| 202 | ylo = ylo - 1000000000; |
| 203 | yhi = yhi + 1; |
| 204 | } |
| 205 | d103 = 0x4189374c; |
| 206 | B34 = ((UINT64) ylo * d103) >> (32 + 8); |
| 207 | B01 = ((UINT64) yhi * d103) >> (32 + 8); |
| 208 | b5 = ylo - B34 * 1000; |
| 209 | b2 = yhi - B01 * 1000; |
| 210 | b3 = ((UINT64) B34 * d103) >> (32 + 8); |
| 211 | b0 = ((UINT64) B01 * d103) >> (32 + 8); |
| 212 | b4 = (unsigned int) B34 - (unsigned int) b3 *1000; |
| 213 | b1 = (unsigned int) B01 - (unsigned int) dm103[b0]; |
| 214 | dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1]; |
| 215 | if (b0 >= 8) /* is b0 8 or 9? */ |
| 216 | res = sign | ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) | |
| 217 | (exp & 0xff)) << 50) | dcoeff; |
| 218 | else /* else b0 is 0..7 */ |
| 219 | res = sign | ((((exp >> 8) << 11) | (b0 << 8) | |
| 220 | (exp & 0xff)) << 50) | dcoeff; |
| 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 | *pres = x; |
| 241 | return; |
| 242 | } |
| 243 | /* normal number */ |
| 244 | if ((comb & 0x1800) == 0x1800) { /* G0..G1 = 11 -> d0 = 8 + G4 */ |
| 245 | d0 = d2b6[((comb >> 8) & 1) | 8]; /* d0 = (comb & 0x0100 ? 9 : 8); */ |
| 246 | exp = (comb & 0x600) >> 1; /* exp = (comb & 0x0400 ? 1 : 0) * 0x200 + |
| 247 | (comb & 0x0200 ? 1 : 0) * 0x100; exp leading bits are G2..G3 */ |
| 248 | } else { |
| 249 | d0 = d2b6[(comb >> 8) & 0x7]; |
| 250 | exp = (comb & 0x1800) >> 3; /* exp = (comb & 0x1000 ? 1 : 0) * 0x200 + |
| 251 | (comb & 0x0800 ? 1 : 0) * 0x100; exp loading bits are G0..G1 */ |
| 252 | } |
| 253 | d1 = d2b5[(trailing >> 40) & 0x3ff]; |
| 254 | d2 = d2b4[(trailing >> 30) & 0x3ff]; |
| 255 | d3 = d2b3[(trailing >> 20) & 0x3ff]; |
| 256 | d4 = d2b2[(trailing >> 10) & 0x3ff]; |
| 257 | d5 = d2b[(trailing) & 0x3ff]; |
| 258 | bcoeff = (d5 + d4 + d3) + d2 + d1 + d0; |
| 259 | exp += (comb & 0xff); |
| 260 | mask = 1; |
| 261 | mask <<= 53; |
| 262 | if (bcoeff < mask) { /* check whether coefficient fits in 10*5+3 bits */ |
| 263 | res = exp; |
| 264 | res <<= 53; |
| 265 | res |= (bcoeff | sign); |
| 266 | *pres = res; |
| 267 | return; |
| 268 | } |
| 269 | /* special format */ |
| 270 | res = (exp << 51) | (sign | 0x6000000000000000ull); |
| 271 | /* add coeff, without leading bits */ |
| 272 | mask = (mask >> 2) - 1; |
| 273 | bcoeff &= mask; |
| 274 | res |= bcoeff; |
| 275 | *pres = res; |
| 276 | } |
| 277 | |
| 278 | void _bid_to_dpd128 (_Decimal128 *, _Decimal128 *); |
| 279 | |
| 280 | void |
| 281 | _bid_to_dpd128 (_Decimal128 *pres, _Decimal128 *px) { |
| 282 | UINT128 res; |
| 283 | UINT128 sign; |
| 284 | unsigned int comb; |
| 285 | UINT128 bcoeff; |
| 286 | UINT128 dcoeff; |
| 287 | UINT128 BH, d1018, BT2, BT1; |
| 288 | UINT64 exp, BL, d109; |
| 289 | UINT64 d106, d103; |
| 290 | UINT64 k1, k2, k4, k5, k7, k8, k10, k11; |
| 291 | unsigned int BHH32, BLL32, BHL32, BLH32, k0, k3, k6, k9, amount; |
| 292 | _Decimal128 x = *px; |
| 293 | |
| 294 | sign.w[1] = (x.w[1] & 0x8000000000000000ull); |
| 295 | sign.w[0] = 0; |
| 296 | comb = (x.w[1] /*& 0x7fffc00000000000ull */ ) >> 46; |
| 297 | exp = 0; |
| 298 | if ((comb & 0x1e000) == 0x1e000) { |
| 299 | res = x; |
| 300 | } else { /* normal number */ |
| 301 | exp = ((x.w[1] & 0x7fff000000000000ull) >> 49) & 0x3fff; |
| 302 | bcoeff.w[1] = (x.w[1] & 0x0001ffffffffffffull); |
| 303 | bcoeff.w[0] = x.w[0]; |
| 304 | d1018 = reciprocals10_128[18]; |
| 305 | __mul_128x128_high (BH, bcoeff, d1018); |
| 306 | amount = recip_scale[18]; |
| 307 | BH.w[0] = (BH.w[0] >> amount) | (BH.w[1] << (64 - amount)); |
| 308 | BL = bcoeff.w[0] - BH.w[0] * 1000000000000000000ull; |
| 309 | d109 = reciprocals10_64[9]; |
| 310 | __mul_64x64_to_128 (BT1, BH.w[0], d109); |
| 311 | BHH32 = (unsigned int) (BT1.w[1] >> short_recip_scale[9]); |
| 312 | BHL32 = (unsigned int) BH.w[0] - BHH32 * 1000000000; |
| 313 | __mul_64x64_to_128 (BT2, BL, d109); |
| 314 | BLH32 = (unsigned int) (BT2.w[1] >> short_recip_scale[9]); |
| 315 | BLL32 = (unsigned int) BL - BLH32 * 1000000000; |
| 316 | d106 = 0x431BDE83; |
| 317 | d103 = 0x4189374c; |
| 318 | k0 = ((UINT64) BHH32 * d106) >> (32 + 18); |
| 319 | BHH32 -= (unsigned int) k0 *1000000; |
| 320 | k1 = ((UINT64) BHH32 * d103) >> (32 + 8); |
| 321 | k2 = BHH32 - (unsigned int) k1 *1000; |
| 322 | k3 = ((UINT64) BHL32 * d106) >> (32 + 18); |
| 323 | BHL32 -= (unsigned int) k3 *1000000; |
| 324 | k4 = ((UINT64) BHL32 * d103) >> (32 + 8); |
| 325 | k5 = BHL32 - (unsigned int) k4 *1000; |
| 326 | k6 = ((UINT64) BLH32 * d106) >> (32 + 18); |
| 327 | BLH32 -= (unsigned int) k6 *1000000; |
| 328 | k7 = ((UINT64) BLH32 * d103) >> (32 + 8); |
| 329 | k8 = BLH32 - (unsigned int) k7 *1000; |
| 330 | k9 = ((UINT64) BLL32 * d106) >> (32 + 18); |
| 331 | BLL32 -= (unsigned int) k9 *1000000; |
| 332 | k10 = ((UINT64) BLL32 * d103) >> (32 + 8); |
| 333 | k11 = BLL32 - (unsigned int) k10 *1000; |
| 334 | dcoeff.w[1] = (b2d[k5] >> 4) | (b2d[k4] << 6) | (b2d[k3] << 16) | |
| 335 | (b2d[k2] << 26) | (b2d[k1] << 36); |
| 336 | dcoeff.w[0] = b2d[k11] | (b2d[k10] << 10) | (b2d[k9] << 20) | |
| 337 | (b2d[k8] << 30) | (b2d[k7] << 40) | (b2d[k6] << 50) | (b2d[k5] << 60); |
| 338 | res.w[0] = dcoeff.w[0]; |
| 339 | if (k0 >= 8) { |
| 340 | res.w[1] = sign.w[1] | ((0x18000 | ((exp >> 12) << 13) | |
| 341 | ((k0 & 1) << 12) | (exp & 0xfff)) << 46) | dcoeff.w[1]; |
| 342 | } else { |
| 343 | res.w[1] = sign.w[1] | ((((exp >> 12) << 15) | (k0 << 12) | |
| 344 | (exp & 0xfff)) << 46) | dcoeff.w[1]; |
| 345 | } |
| 346 | } |
| 347 | *pres = res; |
| 348 | } |
| 349 | |
| 350 | void _dpd_to_bid128 (_Decimal128 *, _Decimal128 *); |
| 351 | |
| 352 | void |
| 353 | _dpd_to_bid128 (_Decimal128 *pres, _Decimal128 *px) { |
| 354 | UINT128 res; |
| 355 | UINT128 sign; |
| 356 | UINT64 exp, comb; |
| 357 | UINT128 trailing; |
| 358 | UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11; |
| 359 | UINT128 bcoeff; |
| 360 | UINT64 tl, th; |
| 361 | _Decimal128 x = *px; |
| 362 | |
| 363 | sign.w[1] = (x.w[1] & 0x8000000000000000ull); |
| 364 | sign.w[0] = 0; |
| 365 | comb = (x.w[1] & 0x7fffc00000000000ull) >> 46; |
| 366 | trailing.w[1] = x.w[1]; |
| 367 | trailing.w[0] = x.w[0]; |
| 368 | if ((comb & 0x1e000) == 0x1e000) { |
| 369 | *pres = x; |
| 370 | return; |
| 371 | } |
| 372 | if ((comb & 0x18000) == 0x18000) { /* G0..G1 = 11 -> d0 = 8 + G4 */ |
| 373 | d0 = d2b6[8 + ((comb & 0x01000) >> 12)]; |
| 374 | exp = (comb & 0x06000) >> 1; /* exp leading bits are G2..G3 */ |
| 375 | } else { |
| 376 | d0 = d2b6[((comb & 0x07000) >> 12)]; |
| 377 | exp = (comb & 0x18000) >> 3; /* exp loading bits are G0..G1 */ |
| 378 | } |
| 379 | d11 = d2b[(trailing.w[0]) & 0x3ff]; |
| 380 | d10 = d2b2[(trailing.w[0] >> 10) & 0x3ff]; |
| 381 | d9 = d2b3[(trailing.w[0] >> 20) & 0x3ff]; |
| 382 | d8 = d2b4[(trailing.w[0] >> 30) & 0x3ff]; |
| 383 | d7 = d2b5[(trailing.w[0] >> 40) & 0x3ff]; |
| 384 | d6 = d2b6[(trailing.w[0] >> 50) & 0x3ff]; |
| 385 | d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)]; |
| 386 | d4 = d2b2[(trailing.w[1] >> 6) & 0x3ff]; |
| 387 | d3 = d2b3[(trailing.w[1] >> 16) & 0x3ff]; |
| 388 | d2 = d2b4[(trailing.w[1] >> 26) & 0x3ff]; |
| 389 | d1 = d2b5[(trailing.w[1] >> 36) & 0x3ff]; |
| 390 | tl = d11 + d10 + d9 + d8 + d7 + d6; |
| 391 | th = d5 + d4 + d3 + d2 + d1 + d0; |
| 392 | __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull); |
| 393 | __add_128_64 (bcoeff, bcoeff, tl); |
| 394 | exp += (comb & 0xfff); |
| 395 | res.w[0] = bcoeff.w[0]; |
| 396 | res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1]; |
| 397 | *pres = res; |
| 398 | } |