Subversion Repositories scarts

/* 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;

}