libdecnumber/bid/bid2dpd_dpd2bid.c

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