libdecnumber/bid/bid2dpd_dpd2bid.c

   1 /* Copyright (C) 2007-2024 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 }