libgcc/config/libbid/bid_dpd.c

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