libgcc/config/libbid/bid_dpd.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 #define  DECNUMDIGITS 34        // work with up to 34 digits
  30
  31 #include "bid_internal.h"
  32 #include "bid_b2d.h"
  33
  34 UINT32
  35 bid_to_bid32 (UINT32 ba) {
  36   UINT32 res;
  37   UINT32 sign, comb, exp;
  38   UINT32 trailing;
  39   UINT32 bcoeff;
  40
  41   sign = (ba & 0x80000000ul);
  42   comb = (ba & 0x7ff00000ul) >> 20;
  43   trailing = (ba & 0x000ffffful);
  44
  45   if ((comb & 0x780) == 0x780) {
  46     ba &= 0xfff0000ul;
  47     return ba;
  48   } else {
  49     if ((comb & 0x600) == 0x600) {      // G0..G1 = 11 -> exp is G2..G11
  50       exp = (comb >> 1) & 0xff;
  51       bcoeff = ((8 + (comb & 1)) << 20) | trailing;
  52     } else {
  53       exp = (comb >> 3) & 0xff;
  54       bcoeff = ((comb & 7) << 20) | trailing;
  55     }
  56     if (bcoeff >= 10000000)
  57       bcoeff = 0;
  58     res = very_fast_get_BID32 (sign, exp, bcoeff);
  59   }
  60   return res;
  61 }
  62
  63 UINT64
  64 bid_to_bid64 (UINT64 ba) {
  65   UINT64 res;
  66   UINT64 sign, comb, exp;
  67   UINT64 trailing;
  68   UINT64 bcoeff;
  69
  70   sign = (ba & 0x8000000000000000ull);
  71   comb = (ba & 0x7ffc000000000000ull) >> 50;
  72   trailing = (ba & 0x0003ffffffffffffull);
  73
  74   if ((comb & 0x1e00) == 0x1e00) {
  75     ba &= 0xfff000000000000ULL;
  76     return ba;
  77   } else {
  78     if ((comb & 0x1800) == 0x1800) {    // G0..G1 = 11 -> exp is G2..G11
  79       exp = (comb >> 1) & 0x3ff;
  80       bcoeff = ((8 + (comb & 1)) << 50) | trailing;
  81     } else {
  82       exp = (comb >> 3) & 0x3ff;
  83       bcoeff = ((comb & 7) << 50) | trailing;
  84     }
  85     if (bcoeff >= 10000000000000000ull)
  86       bcoeff = 0ull;
  87     res = very_fast_get_BID64 (sign, exp, bcoeff);
  88   }
  89   return res;
  90 }
  91
  92 #if DECIMAL_CALL_BY_REFERENCE
  93 void
  94 bid_to_dpd32 (UINT32 * pres, UINT32 * pba) {
  95   UINT32 ba = *pba;
  96 #else
  97 UINT32
  98 bid_to_dpd32 (UINT32 ba) {
  99 #endif
 100   UINT32 res;
 101
 102   UINT32 sign, comb, exp, trailing;
 103   UINT32 b0, b1, b2;
 104   UINT32 bcoeff, dcoeff;
 105   UINT32 nanb = 0;
 106
 107   sign = (ba & 0x80000000);
 108   comb = (ba & 0x7ff00000) >> 20;
 109   trailing = (ba & 0xfffff);
 110
 111   // Detect infinity, and return canonical infinity
 112   if ((comb & 0x7c0) == 0x780) {
 113     res = sign | 0x78000000;
 114     BID_RETURN (res);
 115     // Detect NaN, and canonicalize trailing
 116   } else if ((comb & 0x7c0) == 0x7c0) {
 117     if (trailing > 999999)
 118       trailing = 0;
 119     nanb = ba & 0xfe000000;
 120     exp = 0;
 121     bcoeff = trailing;
 122   } else {      // Normal number
 123     if ((comb & 0x600) == 0x600) {      // G0..G1 = 11 -> exp is G2..G11
 124       exp = (comb >> 1) & 0xff;
 125       bcoeff = ((8 + (comb & 1)) << 20) | trailing;
 126     } else {
 127       exp = (comb >> 3) & 0xff;
 128       bcoeff = ((comb & 7) << 20) | trailing;
 129     }
 130     // Zero the coefficient if non-canonical (>= 10^7)
 131     if (bcoeff >= 10000000)
 132       bcoeff = 0;
 133   }
 134
 135   b0 = bcoeff / 1000000;
 136   b1 = (bcoeff / 1000) % 1000;
 137   b2 = bcoeff % 1000;
 138   dcoeff = (b2d[b1] << 10) | b2d[b2];
 139
 140   if (b0 >= 8)  // is b0 8 or 9?
 141     res =
 142       sign |
 143       ((0x600 | ((exp >> 6) << 7) | ((b0 & 1) << 6) | (exp & 0x3f)) <<
 144        20) | dcoeff;
 145   else  // else b0 is 0..7
 146     res =
 147       sign | ((((exp >> 6) << 9) | (b0 << 6) | (exp & 0x3f)) << 20) |
 148       dcoeff;
 149
 150   res |= nanb;
 151
 152   BID_RETURN (res);
 153 }
 154
 155 #if DECIMAL_CALL_BY_REFERENCE
 156 void
 157 bid_to_dpd64 (UINT64 * pres, UINT64 * pba) {
 158   UINT64 ba = *pba;
 159 #else
 160 UINT64
 161 bid_to_dpd64 (UINT64 ba) {
 162 #endif
 163   UINT64 res;
 164
 165   UINT64 sign, comb, exp;
 166   UINT64 trailing;
 167   UINT32 b0, b1, b2, b3, b4, b5;
 168   UINT64 bcoeff;
 169   UINT64 dcoeff;
 170   UINT32 yhi, ylo;
 171   UINT64 nanb = 0;
 172
 173 //printf("arg bid "FMT_LLX16" \n", ba);
 174   sign = (ba & 0x8000000000000000ull);
 175   comb = (ba & 0x7ffc000000000000ull) >> 50;
 176   trailing = (ba & 0x0003ffffffffffffull);
 177
 178   // Detect infinity, and return canonical infinity
 179   if ((comb & 0x1f00) == 0x1e00) {
 180     res = sign | 0x7800000000000000ull;
 181     BID_RETURN (res);
 182     // Detect NaN, and canonicalize trailing
 183   } else if ((comb & 0x1e00) == 0x1e00) {
 184     if (trailing > 999999999999999ull)
 185       trailing = 0;
 186     nanb = ba & 0xfe00000000000000ull;
 187     exp = 0;
 188     bcoeff = trailing;
 189   } else {      // Normal number
 190     if ((comb & 0x1800) == 0x1800) {    // G0..G1 = 11 -> exp is G2..G11
 191       exp = (comb >> 1) & 0x3ff;
 192       bcoeff = ((8 + (comb & 1)) << 50) | trailing;
 193     } else {
 194       exp = (comb >> 3) & 0x3ff;
 195       bcoeff = ((comb & 7) << 50) | trailing;
 196     }
 197
 198     // Zero the coefficient if it is non-canonical (>= 10^16)
 199     if (bcoeff >= 10000000000000000ull)
 200       bcoeff = 0;
 201   }
 202
 203 // Floor(2^61 / 10^9)
 204 #define D61 (2305843009ull)
 205
 206 // Multipy the binary coefficient by ceil(2^64 / 1000), and take the upper
 207 // 64-bits in order to compute a division by 1000.
 208
 209 #if 1
 210   yhi =
 211     ((UINT64) D61 *
 212      (UINT64) (UINT32) (bcoeff >> (UINT64) 27)) >> (UINT64) 34;
 213   ylo = bcoeff - 1000000000ull * yhi;
 214   if (ylo >= 1000000000) {
 215     ylo = ylo - 1000000000;
 216     yhi = yhi + 1;
 217   }
 218 #else
 219   yhi = bcoeff / 1000000000ull;
 220   ylo = bcoeff % 1000000000ull;
 221 #endif
 222
 223   // yhi = ABBBCCC ylo = DDDEEEFFF
 224   b5 = ylo % 1000;      // b5 = FFF
 225   b3 = ylo / 1000000;   // b3 = DDD
 226   b4 = (ylo / 1000) - (1000 * b3);      // b4 = EEE
 227   b2 = yhi % 1000;      // b2 = CCC
 228   b0 = yhi / 1000000;   // b0 = A
 229   b1 = (yhi / 1000) - (1000 * b0);      // b1 = BBB
 230
 231   dcoeff = b2d[b5] | b2d2[b4] | b2d3[b3] | b2d4[b2] | b2d5[b1];
 232
 233   if (b0 >= 8)  // is b0 8 or 9?
 234     res =
 235       sign |
 236       ((0x1800 | ((exp >> 8) << 9) | ((b0 & 1) << 8) | (exp & 0xff)) <<
 237        50) | dcoeff;
 238   else  // else b0 is 0..7
 239     res =
 240       sign | ((((exp >> 8) << 11) | (b0 << 8) | (exp & 0xff)) << 50) |
 241       dcoeff;
 242
 243   res |= nanb;
 244
 245   BID_RETURN (res);
 246 }
 247
 248 #if DECIMAL_CALL_BY_REFERENCE
 249 void
 250 dpd_to_bid32 (UINT32 * pres, UINT32 * pda) {
 251   UINT32 da = *pda;
 252 #else
 253 UINT32
 254 dpd_to_bid32 (UINT32 da) {
 255 #endif
 256   UINT32 in = *(UINT32 *) & da;
 257   UINT32 res;
 258
 259   UINT32 sign, comb, exp;
 260   UINT32 trailing;
 261   UINT32 d0 = 0, d1, d2;
 262   UINT64 bcoeff;
 263   UINT32 nanb = 0;
 264
 265   sign = (in & 0x80000000);
 266   comb = (in & 0x7ff00000) >> 20;
 267   trailing = (in & 0x000fffff);
 268
 269   if ((comb & 0x7e0) == 0x780) {        // G0..G4 = 1111 -> Inf
 270     res = in & 0xf8000000;
 271     BID_RETURN (res);
 272   } else if ((comb & 0x7c0) == 0x7c0) { // G0..G5 = 11111 -> NaN
 273     nanb = in & 0xfe000000;
 274     exp = 0;
 275   } else {      // Normal number
 276     if ((comb & 0x600) == 0x600) {      // G0..G1 = 11 -> d0 = 8 + G4
 277       d0 = ((comb >> 6) & 1) | 8;
 278       exp = ((comb & 0x180) >> 1) | (comb & 0x3f);
 279     } else {
 280       d0 = (comb >> 6) & 0x7;
 281       exp = ((comb & 0x600) >> 3) | (comb & 0x3f);
 282     }
 283   }
 284   d1 = d2b2[(trailing >> 10) & 0x3ff];
 285   d2 = d2b[(trailing) & 0x3ff];
 286
 287   bcoeff = d2 + d1 + (1000000 * d0);
 288   if (bcoeff < 0x800000) {
 289     res = (exp << 23) | bcoeff | sign;
 290   } else {
 291     res = (exp << 21) | sign | 0x60000000 | (bcoeff & 0x1fffff);
 292   }
 293
 294   res |= nanb;
 295
 296   BID_RETURN (res);
 297 }
 298
 299 #if DECIMAL_CALL_BY_REFERENCE
 300 void
 301 dpd_to_bid64 (UINT64 * pres, UINT64 * pda) {
 302   UINT64 da = *pda;
 303 #else
 304 UINT64
 305 dpd_to_bid64 (UINT64 da) {
 306 #endif
 307   UINT64 in = *(UINT64 *) & da;
 308   UINT64 res;
 309
 310   UINT64 sign, comb, exp;
 311   UINT64 trailing;
 312   // UINT64 d0, d1, d2, d3, d4, d5;
 313
 314   UINT64 d1, d2;
 315   UINT32 d0, d3, d4, d5;
 316   UINT64 bcoeff;
 317   UINT64 nanb = 0;
 318
 319 //printf("arg dpd "FMT_LLX16" \n", in);
 320   sign = (in & 0x8000000000000000ull);
 321   comb = (in & 0x7ffc000000000000ull) >> 50;
 322   trailing = (in & 0x0003ffffffffffffull);
 323
 324   if ((comb & 0x1f00) == 0x1e00) {      // G0..G4 = 1111 -> Inf
 325     res = in & 0xf800000000000000ull;
 326     BID_RETURN (res);
 327   } else if ((comb & 0x1f00) == 0x1f00) {       // G0..G5 = 11111 -> NaN
 328     nanb = in & 0xfe00000000000000ull;
 329     exp = 0;
 330     d0 = 0;
 331   } else {      // Normal number
 332     if ((comb & 0x1800) == 0x1800) {    // G0..G1 = 11 -> d0 = 8 + G4
 333       d0 = ((comb >> 8) & 1) | 8;
 334       // d0 = (comb & 0x0100 ? 9 : 8);
 335       exp = (comb & 0x600) >> 1;
 336       // exp = (comb & 0x0400 ? 1 : 0) * 0x200 + (comb & 0x0200 ? 1 : 0) * 0x100; // exp leading bits are G2..G3
 337     } else {
 338       d0 = (comb >> 8) & 0x7;
 339       exp = (comb & 0x1800) >> 3;
 340       // exp = (comb & 0x1000 ? 1 : 0) * 0x200 + (comb & 0x0800 ? 1 : 0) * 0x100; // exp loading bits are G0..G1
 341     }
 342   }
 343   d1 = d2b5[(trailing >> 40) & 0x3ff];
 344   d2 = d2b4[(trailing >> 30) & 0x3ff];
 345   d3 = d2b3[(trailing >> 20) & 0x3ff];
 346   d4 = d2b2[(trailing >> 10) & 0x3ff];
 347   d5 = d2b[(trailing) & 0x3ff];
 348
 349   bcoeff = (d5 + d4 + d3) + d2 + d1 + (1000000000000000ull * d0);
 350   exp += (comb & 0xff);
 351   res = very_fast_get_BID64 (sign, exp, bcoeff);
 352
 353   res |= nanb;
 354
 355   BID_RETURN (res);
 356 }
 357
 358 #if DECIMAL_CALL_BY_REFERENCE
 359 void
 360 bid_to_dpd128 (UINT128 * pres, UINT128 * pba) {
 361   UINT128 ba = *pba;
 362 #else
 363 UINT128
 364 bid_to_dpd128 (UINT128 ba) {
 365 #endif
 366   UINT128 res;
 367
 368   UINT128 sign;
 369   UINT32 comb, exp;
 370   UINT128 trailing;
 371   UINT128 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
 372   UINT128 bcoeff;
 373   UINT128 dcoeff;
 374   UINT64 nanb = 0;
 375
 376   sign.w[1] = (ba.w[HIGH_128W] & 0x8000000000000000ull);
 377   sign.w[0] = 0;
 378   comb = (ba.w[HIGH_128W] & 0x7fffc00000000000ull) >> 46;
 379   trailing.w[1] = (ba.w[HIGH_128W] & 0x00003fffffffffffull);
 380   trailing.w[0] = ba.w[LOW_128W];
 381   exp = 0;
 382
 383   if ((comb & 0x1f000) == 0x1e000) {    // G0..G4 = 1111 -> Inf
 384     res.w[HIGH_128W] = ba.w[HIGH_128W] & 0xf800000000000000ull;
 385     res.w[LOW_128W] = 0;
 386     BID_RETURN (res);
 387     // Detect NaN, and canonicalize trailing
 388   } else if ((comb & 0x1f000) == 0x1f000) {
 389     if ((trailing.w[1] > 0x0000314dc6448d93ULL) ||      // significand is non-canonical
 390         ((trailing.w[1] == 0x0000314dc6448d93ULL)
 391          && (trailing.w[0] >= 0x38c15b0a00000000ULL))
 392         // significand is non-canonical
 393       ) {
 394       trailing.w[1] = trailing.w[0] = 0ull;
 395     }
 396     bcoeff.w[1] = trailing.w[1];
 397     bcoeff.w[0] = trailing.w[0];
 398     nanb = ba.w[HIGH_128W] & 0xfe00000000000000ull;
 399     exp = 0;
 400   } else {      // Normal number
 401     if ((comb & 0x18000) == 0x18000) {  // G0..G1 = 11 -> exp is G2..G11
 402       exp = (comb >> 1) & 0x3fff;
 403       bcoeff.w[1] =
 404         ((UINT64) (8 + (comb & 1)) << (UINT64) 46) | trailing.w[1];
 405       bcoeff.w[0] = trailing.w[0];
 406     } else {
 407       exp = (comb >> 3) & 0x3fff;
 408       bcoeff.w[1] =
 409         ((UINT64) (comb & 7) << (UINT64) 46) | trailing.w[1];
 410       bcoeff.w[0] = trailing.w[0];
 411     }
 412     // Zero the coefficient if non-canonical (>= 10^34)
 413     if (bcoeff.w[1] > 0x1ed09bead87c0ull ||
 414         (bcoeff.w[1] == 0x1ed09bead87c0ull
 415          && bcoeff.w[0] >= 0x378D8E6400000000ull)) {
 416       bcoeff.w[0] = bcoeff.w[1] = 0;
 417     }
 418   }
 419   // Constant 2^128 / 1000 + 1
 420   {
 421     UINT128 t;
 422     UINT64 t2;
 423     UINT128 d1000;
 424     UINT128 b11, b10, b9, b8, b7, b6, b5, b4, b3, b2, b1;
 425     d1000.w[1] = 0x4189374BC6A7EFull;
 426     d1000.w[0] = 0x9DB22D0E56041894ull;
 427     __mul_128x128_high (b11, bcoeff, d1000);
 428     __mul_128x128_high (b10, b11, d1000);
 429     __mul_128x128_high (b9, b10, d1000);
 430     __mul_128x128_high (b8, b9, d1000);
 431     __mul_128x128_high (b7, b8, d1000);
 432     __mul_128x128_high (b6, b7, d1000);
 433     __mul_128x128_high (b5, b6, d1000);
 434     __mul_128x128_high (b4, b5, d1000);
 435     __mul_128x128_high (b3, b4, d1000);
 436     __mul_128x128_high (b2, b3, d1000);
 437     __mul_128x128_high (b1, b2, d1000);
 438
 439
 440     __mul_64x128_full (t2, t, 1000ull, b11);
 441     __sub_128_128 (d11, bcoeff, t);
 442     __mul_64x128_full (t2, t, 1000ull, b10);
 443     __sub_128_128 (d10, b11, t);
 444     __mul_64x128_full (t2, t, 1000ull, b9);
 445     __sub_128_128 (d9, b10, t);
 446     __mul_64x128_full (t2, t, 1000ull, b8);
 447     __sub_128_128 (d8, b9, t);
 448     __mul_64x128_full (t2, t, 1000ull, b7);
 449     __sub_128_128 (d7, b8, t);
 450     __mul_64x128_full (t2, t, 1000ull, b6);
 451     __sub_128_128 (d6, b7, t);
 452     __mul_64x128_full (t2, t, 1000ull, b5);
 453     __sub_128_128 (d5, b6, t);
 454     __mul_64x128_full (t2, t, 1000ull, b4);
 455     __sub_128_128 (d4, b5, t);
 456     __mul_64x128_full (t2, t, 1000ull, b3);
 457     __sub_128_128 (d3, b4, t);
 458     __mul_64x128_full (t2, t, 1000ull, b2);
 459     __sub_128_128 (d2, b3, t);
 460     __mul_64x128_full (t2, t, 1000ull, b1);
 461     __sub_128_128 (d1, b2, t);
 462     d0 = b1;
 463
 464   }
 465
 466   dcoeff.w[0] = b2d[d11.w[0]] | (b2d[d10.w[0]] << 10) |
 467     (b2d[d9.w[0]] << 20) | (b2d[d8.w[0]] << 30) | (b2d[d7.w[0]] << 40) |
 468     (b2d[d6.w[0]] << 50) | (b2d[d5.w[0]] << 60);
 469   dcoeff.w[1] =
 470     (b2d[d5.w[0]] >> 4) | (b2d[d4.w[0]] << 6) | (b2d[d3.w[0]] << 16) |
 471     (b2d[d2.w[0]] << 26) | (b2d[d1.w[0]] << 36);
 472
 473   res.w[0] = dcoeff.w[0];
 474   if (d0.w[0] >= 8) {
 475     res.w[1] =
 476       sign.
 477       w[1] |
 478       ((0x18000 | ((exp >> 12) << 13) | ((d0.w[0] & 1) << 12) |
 479         (exp & 0xfff)) << 46) | dcoeff.w[1];
 480   } else {
 481     res.w[1] =
 482       sign.
 483       w[1] | ((((exp >> 12) << 15) | (d0.w[0] << 12) | (exp & 0xfff))
 484               << 46) | dcoeff.w[1];
 485   }
 486
 487   res.w[1] |= nanb;
 488
 489   BID_SWAP128 (res);
 490   BID_RETURN (res);
 491 }
 492
 493 #if DECIMAL_CALL_BY_REFERENCE
 494 void
 495 dpd_to_bid128 (UINT128 * pres, UINT128 * pda) {
 496   UINT128 da = *pda;
 497 #else
 498 UINT128
 499 dpd_to_bid128 (UINT128 da) {
 500 #endif
 501   UINT128 in = *(UINT128 *) & da;
 502   UINT128 res;
 503
 504   UINT128 sign;
 505   UINT64 exp, comb;
 506   UINT128 trailing;
 507   UINT64 d0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10, d11;
 508   UINT128 bcoeff;
 509   UINT64 tl, th;
 510   UINT64 nanb = 0;
 511
 512   sign.w[1] = (in.w[HIGH_128W] & 0x8000000000000000ull);
 513   sign.w[0] = 0;
 514   comb = (in.w[HIGH_128W] & 0x7fffc00000000000ull) >> 46;
 515   trailing.w[1] = (in.w[HIGH_128W] & 0x00003fffffffffffull);
 516   trailing.w[0] = in.w[LOW_128W];
 517   exp = 0;
 518
 519   if ((comb & 0x1f000) == 0x1e000) {    // G0..G4 = 1111 -> Inf
 520     res.w[HIGH_128W] = in.w[HIGH_128W] & 0xf800000000000000ull;
 521     res.w[LOW_128W] = 0ull;
 522     BID_RETURN (res);
 523   } else if ((comb & 0x1f000) == 0x1f000) {     // G0..G4 = 11111 -> NaN
 524     nanb = in.w[HIGH_128W] & 0xfe00000000000000ull;
 525     exp = 0;
 526     d0 = 0;
 527   } else {      // Normal number
 528     if ((comb & 0x18000) == 0x18000) {  // G0..G1 = 11 -> d0 = 8 + G4
 529       d0 = 8 + (comb & 0x01000 ? 1 : 0);
 530       exp =
 531         (comb & 0x04000 ? 1 : 0) * 0x2000 +
 532         (comb & 0x02000 ? 1 : 0) * 0x1000;
 533       // exp leading bits are G2..G3
 534     } else {
 535       d0 =
 536         4 * (comb & 0x04000 ? 1 : 0) + 2 * (comb & 0x2000 ? 1 : 0) +
 537         (comb & 0x1000 ? 1 : 0);
 538       exp =
 539         (comb & 0x10000 ? 1 : 0) * 0x2000 +
 540         (comb & 0x08000 ? 1 : 0) * 0x1000;
 541       // exp loading bits are G0..G1
 542     }
 543   }
 544
 545   d11 = d2b[(trailing.w[0]) & 0x3ff];
 546   d10 = d2b[(trailing.w[0] >> 10) & 0x3ff];
 547   d9 = d2b[(trailing.w[0] >> 20) & 0x3ff];
 548   d8 = d2b[(trailing.w[0] >> 30) & 0x3ff];
 549   d7 = d2b[(trailing.w[0] >> 40) & 0x3ff];
 550   d6 = d2b[(trailing.w[0] >> 50) & 0x3ff];
 551   d5 = d2b[(trailing.w[0] >> 60) | ((trailing.w[1] & 0x3f) << 4)];
 552   d4 = d2b[(trailing.w[1] >> 6) & 0x3ff];
 553   d3 = d2b[(trailing.w[1] >> 16) & 0x3ff];
 554   d2 = d2b[(trailing.w[1] >> 26) & 0x3ff];
 555   d1 = d2b[(trailing.w[1] >> 36) & 0x3ff];
 556
 557   tl =
 558     d11 + (d10 * 1000ull) + (d9 * 1000000ull) + (d8 * 1000000000ull) +
 559     (d7 * 1000000000000ull) + (d6 * 1000000000000000ull);
 560   th =
 561     d5 + (d4 * 1000ull) + (d3 * 1000000ull) + (d2 * 1000000000ull) +
 562     (d1 * 1000000000000ull) + (d0 * 1000000000000000ull);
 563   __mul_64x64_to_128 (bcoeff, th, 1000000000000000000ull);
 564   __add_128_64 (bcoeff, bcoeff, tl);
 565
 566   if (!nanb)
 567     exp += (comb & 0xfff);
 568
 569   res.w[0] = bcoeff.w[0];
 570   res.w[1] = (exp << 49) | sign.w[1] | bcoeff.w[1];
 571
 572   res.w[1] |= nanb;
 573
 574   BID_SWAP128 (res);
 575   BID_RETURN (res);
 576 }
 577
 578 UINT128
 579 bid_to_bid128 (UINT128 bq) {
 580   UINT128 res;
 581   UINT64 sign, comb, exp;
 582   UINT64 trailing;
 583   UINT64 bcoeff;
 584
 585   UINT128 rq;
 586   UINT128 qcoeff;
 587   UINT64 ba, bb;
 588
 589   ba = *((UINT64 *) & bq + HIGH_128W);
 590   bb = *((UINT64 *) & bq + LOW_128W);
 591
 592   sign = (ba & 0x8000000000000000ull);
 593   comb = (ba & 0x7fffc00000000000ull) >> 46;
 594   trailing = (ba & 0x00003fffffffffffull);
 595
 596   if ((comb & 0x18000) == 0x18000) {    // G0..G1 = 11 -> exp is G2..G11
 597     exp = (comb >> 1) & 0x3fff;
 598     bcoeff = ((8 + (comb & 1)) << 46) | trailing;
 599   } else {
 600     exp = (comb >> 3) & 0x3fff;
 601     bcoeff = ((comb & 7) << 46) | trailing;
 602   }
 603
 604   if ((comb & 0x1f000) == 0x1f000) {    //NaN
 605     ba &= 0xfe003fffffffffffULL;        // make exponent 0
 606     bcoeff &= 0x00003fffffffffffull;    // NaN payloat is only T.
 607     if ((bcoeff > 0x0000314dc6448d93ULL) ||     // significand is non-canonical
 608         ((bcoeff == 0x0000314dc6448d93ULL)
 609          && (bb >= 0x38c15b0a00000000ULL))
 610         // significand is non-canonical
 611       ) {
 612       bcoeff = 0ull;
 613       ba &= ~0x00003fffffffffffull;
 614       bb = 0ull;
 615     }
 616     *((UINT64 *) & rq + HIGH_128W) = ba;
 617     *((UINT64 *) & rq + LOW_128W) = bb;
 618     return rq;
 619   } else if ((comb & 0x1e000) == 0x1e000) {     //Inf
 620     ba &= 0xf800000000000000ULL;        // make exponent and significand 0
 621     bb = 0;
 622     *((UINT64 *) & rq + HIGH_128W) = ba;
 623     *((UINT64 *) & rq + LOW_128W) = bb;
 624     return rq;
 625   }
 626
 627   if ((bcoeff > 0x0001ed09bead87c0ull)
 628       || ((bcoeff == 0x0001ed09bead87c0ull)
 629           && (bb > 0x378d8e63ffffffffull))) {
 630     // significand is non-canonical
 631     bcoeff = 0ull;
 632     bb = 0ull;
 633   }
 634
 635   *((UINT64 *) & qcoeff + 1) = bcoeff;
 636   *((UINT64 *) & qcoeff + 0) = bb;
 637
 638   get_BID128_fast (&res, sign, exp, qcoeff);
 639
 640   BID_SWAP128 (res);
 641   return res;
 642 }
 643
 644 UINT32
 645 bid32_canonize (UINT32 ba) {
 646   FPSC bidrnd;
 647   unsigned int rnd = 0;
 648
 649   UINT32 res;
 650   UINT32 sign, comb, exp;
 651   UINT32 trailing;
 652   UINT32 bcoeff;
 653
 654   sign = (ba & 0x80000000);
 655   comb = (ba & 0x7ff00000) >> 20;
 656   trailing = (ba & 0x000fffff);
 657
 658   if ((comb & 0x600) == 0x600) {        // G0..G1 = 11 -> exp is G2..G11
 659     exp = (comb >> 1) & 0xff;
 660     bcoeff = ((8 + (comb & 1)) << 20) | trailing;
 661   } else {
 662     exp = (comb >> 3) & 0xff;
 663     bcoeff = ((comb & 7) << 20) | trailing;
 664   }
 665
 666   if ((comb & 0x7c0) == 0x7c0) {        //NaN
 667     ba &= 0xfe0fffff;   // make exponent 0
 668     bcoeff &= 0x000fffff;       // NaN payloat is only T.
 669     if (bcoeff >= 1000000)
 670       ba &= 0xfff00000; //treat non-canonical significand
 671     return ba;
 672   } else if ((comb & 0x780) == 0x780) { //Inf
 673     ba &= 0xf8000000;   // make exponent and significand 0
 674     return ba;
 675   }
 676
 677   if (bcoeff >= 10000000)
 678     bcoeff = 0;
 679   rnd = bidrnd = ROUNDING_TO_NEAREST;
 680   res = get_BID32 (sign, exp, bcoeff, rnd, &bidrnd);
 681   return res;
 682 }
 683
 684 UINT64
 685 bid64_canonize (UINT64 ba) {
 686   UINT64 res;
 687   UINT64 sign, comb, exp;
 688   UINT64 trailing;
 689   UINT64 bcoeff;
 690
 691   sign = (ba & 0x8000000000000000ull);
 692   comb = (ba & 0x7ffc000000000000ull) >> 50;
 693   trailing = (ba & 0x0003ffffffffffffull);
 694
 695
 696   if ((comb & 0x1800) == 0x1800) {      // G0..G1 = 11 -> exp is G2..G11
 697     exp = (comb >> 1) & 0x3ff;
 698     bcoeff = ((8 + (comb & 1)) << 50) | trailing;
 699   } else {
 700     exp = (comb >> 3) & 0x3ff;
 701     bcoeff = ((comb & 7) << 50) | trailing;
 702   }
 703
 704   if ((comb & 0x1f00) == 0x1f00) {      //NaN
 705     ba &= 0xfe03ffffffffffffULL;        // make exponent 0
 706     bcoeff &= 0x0003ffffffffffffull;    // NaN payloat is only T.
 707     if (bcoeff >= 1000000000000000ull)
 708       ba &= 0xfe00000000000000ull;      // treat non canonical significand and zero G6-G12
 709     return ba;
 710   } else if ((comb & 0x1e00) == 0x1e00) {       //Inf
 711     ba &= 0xf800000000000000ULL;        // make exponent and significand 0
 712     return ba;
 713   }
 714
 715   if (bcoeff >= 10000000000000000ull) {
 716     bcoeff = 0ull;
 717   }
 718   res = very_fast_get_BID64 (sign, exp, bcoeff);
 719   return res;
 720 }
 721
 722 UINT128
 723 bid128_canonize (UINT128 bq) {
 724   UINT128 res;
 725   UINT64 sign, comb, exp;
 726   UINT64 trailing;
 727   UINT64 bcoeff;
 728
 729   UINT128 rq;
 730   UINT128 qcoeff;
 731   UINT64 ba, bb;
 732
 733   ba = *((UINT64 *) & bq + HIGH_128W);
 734   bb = *((UINT64 *) & bq + LOW_128W);
 735
 736   sign = (ba & 0x8000000000000000ull);
 737   comb = (ba & 0x7fffc00000000000ull) >> 46;
 738   trailing = (ba & 0x00003fffffffffffull);
 739
 740   if ((comb & 0x18000) == 0x18000) {    // G0..G1 = 11 -> exp is G2..G11
 741     exp = (comb >> 1) & 0x3fff;
 742     bcoeff = ((8 + (comb & 1)) << 46) | trailing;
 743   } else {
 744     exp = (comb >> 3) & 0x3fff;
 745     bcoeff = ((comb & 7) << 46) | trailing;
 746   }
 747
 748   if ((comb & 0x1f000) == 0x1f000) {    //NaN
 749     ba &= 0xfe003fffffffffffULL;        // make exponent 0
 750     bcoeff &= 0x00003fffffffffffull;    // NaN payload is only T.
 751
 752     if ((bcoeff > 0x0000314dc6448d93ULL) ||     // significand is non-canonical
 753         ((bcoeff == 0x0000314dc6448d93ULL)
 754          && (bb >= 0x38c15b0a00000000ULL))
 755         // significand is non-canonical
 756       ) {
 757       bcoeff = 0ull;
 758       ba &= ~0x00003fffffffffffull;
 759       bb = 0ull;
 760     }
 761     *((UINT64 *) & rq + HIGH_128W) = ba;
 762     *((UINT64 *) & rq + LOW_128W) = bb;
 763     return rq;
 764   } else if ((comb & 0x1e000) == 0x1e000) {     //Inf
 765     ba &= 0xf800000000000000ULL;        // make exponent and significand 0
 766     bb = 0;
 767     *((UINT64 *) & rq + HIGH_128W) = ba;
 768     *((UINT64 *) & rq + LOW_128W) = bb;
 769     return rq;
 770   }
 771
 772   if ((bcoeff > 0x0001ed09bead87c0ull) ||       // significand is non-canonical
 773       ((bcoeff == 0x0001ed09bead87c0ull)
 774        && (bb > 0x378d8e63ffffffffull))
 775       // significand is non-canonical
 776     ) {
 777     bcoeff = 0ull;
 778     bb = 0ull;
 779   }
 780
 781   *((UINT64 *) & qcoeff + 1) = bcoeff;
 782   *((UINT64 *) & qcoeff + 0) = bb;
 783
 784   get_BID128_fast (&res, sign, exp, qcoeff);
 785   BID_SWAP128 (res);
 786   return res;
 787 }