--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libgcc/config/libbid/bid128_fma.c	Fri Jul 17 14:47:48 2009 +0900
@@ -0,0 +1,4460 @@
+/* Copyright (C) 2007, 2009  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 3, or (at your option) any later
+version.
+
+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.
+
+Under Section 7 of GPL version 3, you are granted additional
+permissions described in the GCC Runtime Library Exception, version
+3.1, as published by the Free Software Foundation.
+
+You should have received a copy of the GNU General Public License and
+a copy of the GCC Runtime Library Exception along with this program;
+see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+<http://www.gnu.org/licenses/>.  */
+
+/*****************************************************************************
+ * 
+ *  BID128 fma   x * y + z
+ * 
+ ****************************************************************************/
+
+#include "bid_internal.h"
+
+static void
+rounding_correction (unsigned int rnd_mode,
+        	     unsigned int is_inexact_lt_midpoint,
+        	     unsigned int is_inexact_gt_midpoint,
+        	     unsigned int is_midpoint_lt_even,
+        	     unsigned int is_midpoint_gt_even,
+        	     int unbexp,
+        	     UINT128 * ptrres, _IDEC_flags * ptrfpsf) {
+  // unbiased true exponent unbexp may be larger than emax
+
+  UINT128 res = *ptrres; // expected to have the correct sign and coefficient
+  // (the exponent field is ignored, as unbexp is used instead)
+  UINT64 sign, exp;
+  UINT64 C_hi, C_lo;
+
+  // general correction from RN to RA, RM, RP, RZ
+  // Note: if the result is negative, then is_inexact_lt_midpoint, 
+  // is_inexact_gt_midpoint, is_midpoint_lt_even, and is_midpoint_gt_even 
+  // have to be considered as if determined for the absolute value of the 
+  // result (so they seem to be reversed)
+
+  if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+      is_midpoint_lt_even || is_midpoint_gt_even) {
+    *ptrfpsf |= INEXACT_EXCEPTION;
+  }
+  // apply correction to result calculated with unbounded exponent
+  sign = res.w[1] & MASK_SIGN;
+  exp = (UINT64) (unbexp + 6176) << 49; // valid only if expmin<=unbexp<=expmax
+  C_hi = res.w[1] & MASK_COEFF;
+  C_lo = res.w[0];
+  if ((!sign && ((rnd_mode == ROUNDING_UP && is_inexact_lt_midpoint) ||
+      ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_UP) && 
+      is_midpoint_gt_even))) || 
+      (sign && ((rnd_mode == ROUNDING_DOWN && is_inexact_lt_midpoint) ||
+      ((rnd_mode == ROUNDING_TIES_AWAY || rnd_mode == ROUNDING_DOWN) && 
+      is_midpoint_gt_even)))) {
+    // C = C + 1
+    C_lo = C_lo + 1;
+    if (C_lo == 0)
+      C_hi = C_hi + 1;
+    if (C_hi == 0x0001ed09bead87c0ull && C_lo == 0x378d8e6400000000ull) {
+      // C = 10^34 => rounding overflow
+      C_hi = 0x0000314dc6448d93ull;
+      C_lo = 0x38c15b0a00000000ull; // 10^33
+      // exp = exp + EXP_P1;
+      unbexp = unbexp + 1;
+      exp = (UINT64) (unbexp + 6176) << 49;
+    }
+  } else if ((is_midpoint_lt_even || is_inexact_gt_midpoint) &&
+      ((sign && (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TO_ZERO)) ||
+      (!sign && (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TO_ZERO)))) {
+    // C = C - 1
+    C_lo = C_lo - 1;
+    if (C_lo == 0xffffffffffffffffull)
+      C_hi--;
+    // check if we crossed into the lower decade
+    if (C_hi == 0x0000314dc6448d93ull && C_lo == 0x38c15b09ffffffffull) { 
+      // C = 10^33 - 1
+      if (exp > 0) {
+        C_hi = 0x0001ed09bead87c0ull; // 10^34 - 1
+        C_lo = 0x378d8e63ffffffffull;
+        // exp = exp - EXP_P1;
+        unbexp = unbexp - 1;
+        exp = (UINT64) (unbexp + 6176) << 49;
+      } else { // if exp = 0
+        if (is_midpoint_lt_even || is_midpoint_lt_even ||
+            is_inexact_gt_midpoint || is_inexact_gt_midpoint) // tiny & inexact
+          *ptrfpsf |= UNDERFLOW_EXCEPTION;
+      }
+    }
+  } else {
+    ; // the result is already correct
+  }
+  if (unbexp > expmax) { // 6111
+    *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+    exp = 0;
+    if (!sign) { // result is positive
+      if (rnd_mode == ROUNDING_UP || rnd_mode == ROUNDING_TIES_AWAY) { // +inf
+        C_hi = 0x7800000000000000ull;
+        C_lo = 0x0000000000000000ull;
+      } else { // res = +MAXFP = (10^34-1) * 10^emax
+        C_hi = 0x5fffed09bead87c0ull;
+        C_lo = 0x378d8e63ffffffffull;
+      }
+    } else { // result is negative
+      if (rnd_mode == ROUNDING_DOWN || rnd_mode == ROUNDING_TIES_AWAY) { // -inf
+        C_hi = 0xf800000000000000ull;
+        C_lo = 0x0000000000000000ull;
+      } else { // res = -MAXFP = -(10^34-1) * 10^emax
+        C_hi = 0xdfffed09bead87c0ull;
+        C_lo = 0x378d8e63ffffffffull;
+      }
+    }
+  }
+  // assemble the result
+  res.w[1] = sign | exp | C_hi;
+  res.w[0] = C_lo;
+  *ptrres = res;
+}
+
+static void
+add256 (UINT256 x, UINT256 y, UINT256 * pz) {
+  // *z = x + yl assume the sum fits in 256 bits
+  UINT256 z;
+  z.w[0] = x.w[0] + y.w[0];
+  if (z.w[0] < x.w[0]) {
+    x.w[1]++;
+    if (x.w[1] == 0x0000000000000000ull) {
+      x.w[2]++;
+      if (x.w[2] == 0x0000000000000000ull) {
+        x.w[3]++;
+      }
+    }
+  }
+  z.w[1] = x.w[1] + y.w[1];
+  if (z.w[1] < x.w[1]) {
+    x.w[2]++;
+    if (x.w[2] == 0x0000000000000000ull) {
+      x.w[3]++;
+    }
+  }
+  z.w[2] = x.w[2] + y.w[2];
+  if (z.w[2] < x.w[2]) {
+    x.w[3]++;
+  }
+  z.w[3] = x.w[3] + y.w[3]; // it was assumed that no carry is possible
+  *pz = z;
+}
+
+static void
+sub256 (UINT256 x, UINT256 y, UINT256 * pz) {
+  // *z = x - y; assume x >= y
+  UINT256 z;
+  z.w[0] = x.w[0] - y.w[0];
+  if (z.w[0] > x.w[0]) {
+    x.w[1]--;
+    if (x.w[1] == 0xffffffffffffffffull) {
+      x.w[2]--;
+      if (x.w[2] == 0xffffffffffffffffull) {
+        x.w[3]--;
+      }
+    }
+  }
+  z.w[1] = x.w[1] - y.w[1];
+  if (z.w[1] > x.w[1]) {
+    x.w[2]--;
+    if (x.w[2] == 0xffffffffffffffffull) {
+      x.w[3]--;
+    }
+  }
+  z.w[2] = x.w[2] - y.w[2];
+  if (z.w[2] > x.w[2]) {
+    x.w[3]--;
+  }
+  z.w[3] = x.w[3] - y.w[3]; // no borrow possible, because x >= y
+  *pz = z;
+}
+
+
+static int
+nr_digits256 (UINT256 R256) {
+  int ind;
+  // determine the number of decimal digits in R256
+  if (R256.w[3] == 0x0 && R256.w[2] == 0x0 && R256.w[1] == 0x0) {
+    // between 1 and 19 digits
+    for (ind = 1; ind <= 19; ind++) {
+      if (R256.w[0] < ten2k64[ind]) {
+        break;
+      }
+    }
+    // ind digits
+  } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0 &&
+             (R256.w[1] < ten2k128[0].w[1] ||
+              (R256.w[1] == ten2k128[0].w[1]
+               && R256.w[0] < ten2k128[0].w[0]))) {
+    // 20 digits
+    ind = 20;
+  } else if (R256.w[3] == 0x0 && R256.w[2] == 0x0) {
+    // between 21 and 38 digits
+    for (ind = 1; ind <= 18; ind++) {
+      if (R256.w[1] < ten2k128[ind].w[1] ||
+          (R256.w[1] == ten2k128[ind].w[1] &&
+           R256.w[0] < ten2k128[ind].w[0])) {
+        break;
+      }
+    }
+    // ind + 20 digits
+    ind = ind + 20;
+  } else if (R256.w[3] == 0x0 &&
+             (R256.w[2] < ten2k256[0].w[2] ||
+              (R256.w[2] == ten2k256[0].w[2] &&
+               R256.w[1] < ten2k256[0].w[1]) ||
+              (R256.w[2] == ten2k256[0].w[2] &&
+               R256.w[1] == ten2k256[0].w[1] &&
+               R256.w[0] < ten2k256[0].w[0]))) {
+    // 39 digits
+    ind = 39;
+  } else {
+    // between 40 and 68 digits
+    for (ind = 1; ind <= 29; ind++) {
+      if (R256.w[3] < ten2k256[ind].w[3] ||
+          (R256.w[3] == ten2k256[ind].w[3] &&
+           R256.w[2] < ten2k256[ind].w[2]) ||
+          (R256.w[3] == ten2k256[ind].w[3] &&
+           R256.w[2] == ten2k256[ind].w[2] &&
+           R256.w[1] < ten2k256[ind].w[1]) ||
+          (R256.w[3] == ten2k256[ind].w[3] &&
+           R256.w[2] == ten2k256[ind].w[2] &&
+           R256.w[1] == ten2k256[ind].w[1] &&
+           R256.w[0] < ten2k256[ind].w[0])) {
+        break;
+      }
+    }
+    // ind + 39 digits
+    ind = ind + 39;
+  }
+  return (ind);
+}
+
+// add/subtract C4 and C3 * 10^scale; this may follow a previous rounding, so
+// use the rounding information from ptr_is_* to avoid a double rounding error
+static void
+add_and_round (int q3,
+               int q4,
+               int e4,
+               int delta,
+               int p34,
+               UINT64 z_sign,
+               UINT64 p_sign,
+               UINT128 C3,
+               UINT256 C4,
+               int rnd_mode,
+               int *ptr_is_midpoint_lt_even,
+               int *ptr_is_midpoint_gt_even,
+               int *ptr_is_inexact_lt_midpoint,
+               int *ptr_is_inexact_gt_midpoint,
+               _IDEC_flags * ptrfpsf, UINT128 * ptrres) {
+
+  int scale;
+  int x0;
+  int ind;
+  UINT64 R64;
+  UINT128 P128, R128;
+  UINT192 P192, R192;
+  UINT256 R256;
+  int is_midpoint_lt_even = 0;
+  int is_midpoint_gt_even = 0;
+  int is_inexact_lt_midpoint = 0;
+  int is_inexact_gt_midpoint = 0;
+  int is_midpoint_lt_even0 = 0;
+  int is_midpoint_gt_even0 = 0;
+  int is_inexact_lt_midpoint0 = 0;
+  int is_inexact_gt_midpoint0 = 0;
+  int incr_exp = 0;
+  int is_tiny = 0;
+  int lt_half_ulp = 0;
+  int eq_half_ulp = 0;
+  // int gt_half_ulp = 0;
+  UINT128 res = *ptrres;
+
+  // scale C3 up by 10^(q4-delta-q3), 0 <= q4-delta-q3 <= 2*P34-2 = 66
+  scale = q4 - delta - q3; // 0 <= scale <= 66 (or 0 <= scale <= 68 if this
+  // comes from Cases (2), (3), (4), (5), (6), with 0 <= |delta| <= 1
+
+  // calculate C3 * 10^scale in R256 (it has at most 67 decimal digits for
+  // Cases (15),(16),(17) and at most 69 for Cases (2),(3),(4),(5),(6))
+  if (scale == 0) {
+    R256.w[3] = 0x0ull;
+    R256.w[2] = 0x0ull;
+    R256.w[1] = C3.w[1];
+    R256.w[0] = C3.w[0];
+  } else if (scale <= 19) { // 10^scale fits in 64 bits
+    P128.w[1] = 0;
+    P128.w[0] = ten2k64[scale];
+    __mul_128x128_to_256 (R256, P128, C3);
+  } else if (scale <= 38) { // 10^scale fits in 128 bits
+    __mul_128x128_to_256 (R256, ten2k128[scale - 20], C3);
+  } else if (scale <= 57) { // 39 <= scale <= 57 
+    // 10^scale fits in 192 bits but C3 * 10^scale fits in 223 or 230 bits
+    // (10^67 has 223 bits; 10^69 has 230 bits);  
+    // must split the computation:  
+    // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127
+    // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty
+    // Note that 1 <= scale - 38 <= 19 => 10^(scale-38) fits in 64 bits
+    __mul_64x128_to_128 (R128, ten2k64[scale - 38], C3);
+    // now multiply R128 by 10^38
+    __mul_128x128_to_256 (R256, R128, ten2k128[18]);
+  } else { // 58 <= scale <= 66
+    // 10^scale takes between 193 and 220 bits,
+    // and C3 * 10^scale fits in 223 bits (10^67/10^69 has 223/230 bits)
+    // must split the computation: 
+    // 10^scale * C3 = 10*38 * 10^(scale-38) * C3 where 10^38 takes 127
+    // bits and so 10^(scale-38) * C3 fits in 128 bits with certainty 
+    // Note that 20 <= scale - 38 <= 30 => 10^(scale-38) fits in 128 bits
+    // Calculate first 10^(scale-38) * C3, which fits in 128 bits; because
+    // 10^(scale-38) takes more than 64 bits, C3 will take less than 64
+    __mul_64x128_to_128 (R128, C3.w[0], ten2k128[scale - 58]);
+    // now calculate 10*38 * 10^(scale-38) * C3 
+    __mul_128x128_to_256 (R256, R128, ten2k128[18]);
+  }
+  // C3 * 10^scale is now in R256 
+
+  // for Cases (15), (16), (17) C4 > C3 * 10^scale because C4 has at least 
+  // one extra digit; for Cases (2), (3), (4), (5), or (6) any order is 
+  // possible 
+  // add/subtract C4 and C3 * 10^scale; the exponent is e4
+  if (p_sign == z_sign) { // R256 = C4 + R256
+    // calculate R256 = C4 + C3 * 10^scale = C4 + R256 which is exact,
+    // but may require rounding
+    add256 (C4, R256, &R256);
+  } else { // if (p_sign != z_sign) { // R256 = C4 - R256
+    // calculate R256 = C4 - C3 * 10^scale = C4 - R256 or
+    // R256 = C3 * 10^scale - C4 = R256 - C4 which is exact,
+    // but may require rounding
+
+    // compare first R256 = C3 * 10^scale and C4 
+    if (R256.w[3] > C4.w[3] || (R256.w[3] == C4.w[3] && R256.w[2] > C4.w[2]) ||
+        (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] > C4.w[1]) ||
+        (R256.w[3] == C4.w[3] && R256.w[2] == C4.w[2] && R256.w[1] == C4.w[1] &&
+        R256.w[0] >= C4.w[0])) { // C3 * 10^scale >= C4
+      // calculate R256 = C3 * 10^scale - C4 = R256 - C4, which is exact,
+      // but may require rounding 
+      sub256 (R256, C4, &R256);
+      // flip p_sign too, because the result has the sign of z 
+      p_sign = z_sign;
+    } else { // if C4 > C3 * 10^scale
+      // calculate R256 = C4 - C3 * 10^scale = C4 - R256, which is exact,
+      // but may require rounding  
+      sub256 (C4, R256, &R256);
+    }
+    // if the result is pure zero, the sign depends on the rounding mode
+    // (x*y and z had opposite signs) 
+    if (R256.w[3] == 0x0ull && R256.w[2] == 0x0ull &&
+        R256.w[1] == 0x0ull && R256.w[0] == 0x0ull) {
+      if (rnd_mode != ROUNDING_DOWN)
+        p_sign = 0x0000000000000000ull;
+      else
+        p_sign = 0x8000000000000000ull;
+      // the exponent is max (e4, expmin)
+      if (e4 < -6176)
+        e4 = expmin;
+      // assemble result 
+      res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49);
+      res.w[0] = 0x0;
+      *ptrres = res;
+      return;
+    }
+  }
+
+  // determine the number of decimal digits in R256
+  ind = nr_digits256 (R256);
+
+  // the exact result is (-1)^p_sign * R256 * 10^e4 where q (R256) = ind;
+  // round to the destination precision, with unbounded exponent
+
+  if (ind <= p34) {
+    // result rounded to the destination precision with unbounded exponent
+    // is exact
+    if (ind + e4 < p34 + expmin) {
+      is_tiny = 1; // applies to all rounding modes
+    }
+    res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R256.w[1];
+    res.w[0] = R256.w[0];
+    // Note: res is correct only if expmin <= e4 <= expmax
+  } else { // if (ind > p34)
+    // if more than P digits, round to nearest to P digits
+    // round R256 to p34 digits
+    x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68
+    if (ind <= 38) {
+      P128.w[1] = R256.w[1];
+      P128.w[0] = R256.w[0];
+      round128_19_38 (ind, x0, P128, &R128, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+    } else if (ind <= 57) {
+      P192.w[2] = R256.w[2];
+      P192.w[1] = R256.w[1];
+      P192.w[0] = R256.w[0];
+      round192_39_57 (ind, x0, P192, &R192, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+      R128.w[1] = R192.w[1];
+      R128.w[0] = R192.w[0];
+    } else { // if (ind <= 68)
+      round256_58_76 (ind, x0, R256, &R256, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+      R128.w[1] = R256.w[1];
+      R128.w[0] = R256.w[0];
+    }
+    // the rounded result has p34 = 34 digits
+    e4 = e4 + x0 + incr_exp;
+    if (rnd_mode == ROUNDING_TO_NEAREST) {
+      if (e4 < expmin) {
+        is_tiny = 1; // for other rounding modes apply correction
+      }
+    } else {
+      // for RM, RP, RZ, RA apply correction in order to determine tininess
+      // but do not save the result; apply the correction to 
+      // (-1)^p_sign * significand * 10^0
+      P128.w[1] = p_sign | 0x3040000000000000ull | R128.w[1];
+      P128.w[0] = R128.w[0];
+      rounding_correction (rnd_mode,
+        		   is_inexact_lt_midpoint,
+        		   is_inexact_gt_midpoint, is_midpoint_lt_even,
+        		   is_midpoint_gt_even, 0, &P128, ptrfpsf);
+      scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
+      // the number of digits in the significand is p34 = 34
+      if (e4 + scale < expmin) {
+        is_tiny = 1;
+      }
+    }
+    ind = p34; // the number of decimal digits in the signifcand of res
+    res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | R128.w[1]; // RN
+    res.w[0] = R128.w[0];
+    // Note: res is correct only if expmin <= e4 <= expmax
+    // set the inexact flag after rounding with bounded exponent, if any
+  }
+  // at this point we have the result rounded with unbounded exponent in
+  // res and we know its tininess:
+  // res = (-1)^p_sign * significand * 10^e4, 
+  // where q (significand) = ind <= p34
+  // Note: res is correct only if expmin <= e4 <= expmax
+
+  // check for overflow if RN
+  if (rnd_mode == ROUNDING_TO_NEAREST && (ind + e4) > (p34 + expmax)) {
+    res.w[1] = p_sign | 0x7800000000000000ull;
+    res.w[0] = 0x0000000000000000ull;
+    *ptrres = res;
+    *ptrfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+    return; // BID_RETURN (res)
+  } // else not overflow or not RN, so continue
+
+  // if (e4 >= expmin) we have the result rounded with bounded exponent
+  if (e4 < expmin) {
+    x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res
+    // where the result rounded [at most] once is
+    //   (-1)^p_sign * significand_res * 10^e4
+
+    // avoid double rounding error
+    is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+    is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+    is_midpoint_lt_even0 = is_midpoint_lt_even;
+    is_midpoint_gt_even0 = is_midpoint_gt_even;
+    is_inexact_lt_midpoint = 0;
+    is_inexact_gt_midpoint = 0;
+    is_midpoint_lt_even = 0;
+    is_midpoint_gt_even = 0;
+
+    if (x0 > ind) {
+      // nothing is left of res when moving the decimal point left x0 digits
+      is_inexact_lt_midpoint = 1;
+      res.w[1] = p_sign | 0x0000000000000000ull;
+      res.w[0] = 0x0000000000000000ull;
+      e4 = expmin;
+    } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34
+      // this is <, =, or > 1/2 ulp
+      // compare the ind-digit value in the significand of res with
+      // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is 
+      // less than, equal to, or greater than 1/2 ulp (significand of res)
+      R128.w[1] = res.w[1] & MASK_COEFF;
+      R128.w[0] = res.w[0];
+      if (ind <= 19) {
+        if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
+          lt_half_ulp = 1;
+          is_inexact_lt_midpoint = 1;
+        } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+          is_midpoint_gt_even = 1;
+        } else { // > 1/2 ulp
+          // gt_half_ulp = 1;
+          is_inexact_gt_midpoint = 1;
+        }
+      } else { // if (ind <= 38) {
+        if (R128.w[1] < midpoint128[ind - 20].w[1] || 
+            (R128.w[1] == midpoint128[ind - 20].w[1] && 
+            R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
+          lt_half_ulp = 1;
+          is_inexact_lt_midpoint = 1;
+        } else if (R128.w[1] == midpoint128[ind - 20].w[1] && 
+            R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+          is_midpoint_gt_even = 1;
+        } else { // > 1/2 ulp
+          // gt_half_ulp = 1;
+          is_inexact_gt_midpoint = 1;
+        }
+      }
+      if (lt_half_ulp || eq_half_ulp) {
+        // res = +0.0 * 10^expmin
+        res.w[1] = 0x0000000000000000ull;
+        res.w[0] = 0x0000000000000000ull;
+      } else { // if (gt_half_ulp)
+        // res = +1 * 10^expmin
+        res.w[1] = 0x0000000000000000ull;
+        res.w[0] = 0x0000000000000001ull;
+      }
+      res.w[1] = p_sign | res.w[1];
+      e4 = expmin;
+    } else { // if (1 <= x0 <= ind - 1 <= 33)
+      // round the ind-digit result to ind - x0 digits
+
+      if (ind <= 18) { // 2 <= ind <= 18
+        round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+        res.w[1] = 0x0;
+        res.w[0] = R64;
+      } else if (ind <= 38) {
+        P128.w[1] = res.w[1] & MASK_COEFF;
+        P128.w[0] = res.w[0];
+        round128_19_38 (ind, x0, P128, &res, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+      }
+      e4 = e4 + x0; // expmin
+      // we want the exponent to be expmin, so if incr_exp = 1 then
+      // multiply the rounded result by 10 - it will still fit in 113 bits
+      if (incr_exp) {
+        // 64 x 128 -> 128
+        P128.w[1] = res.w[1] & MASK_COEFF;
+        P128.w[0] = res.w[0];
+        __mul_64x128_to_128 (res, ten2k64[1], P128);
+      }
+      res.w[1] =
+        p_sign | ((UINT64) (e4 + 6176) << 49) | (res.w[1] & MASK_COEFF);
+      // avoid a double rounding error
+      if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+          is_midpoint_lt_even) { // double rounding error upward
+        // res = res - 1
+        res.w[0]--;
+        if (res.w[0] == 0xffffffffffffffffull)
+          res.w[1]--;
+        // Note: a double rounding error upward is not possible; for this
+        // the result after the first rounding would have to be 99...95
+        // (35 digits in all), possibly followed by a number of zeros; this
+        // is not possible in Cases (2)-(6) or (15)-(17) which may get here
+        is_midpoint_lt_even = 0;
+        is_inexact_lt_midpoint = 1;
+      } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+          is_midpoint_gt_even) { // double rounding error downward
+        // res = res + 1
+        res.w[0]++;
+        if (res.w[0] == 0)
+          res.w[1]++;
+        is_midpoint_gt_even = 0;
+        is_inexact_gt_midpoint = 1;
+      } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+        // if this second rounding was exact the result may still be 
+        // inexact because of the first rounding
+        if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+          is_inexact_gt_midpoint = 1;
+        }
+        if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+          is_inexact_lt_midpoint = 1;
+        }
+      } else if (is_midpoint_gt_even &&
+        	 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+        // pulled up to a midpoint
+        is_inexact_lt_midpoint = 1;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+      } else if (is_midpoint_lt_even &&
+        	 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+        // pulled down to a midpoint
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 1;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+      } else {
+        ;
+      }
+    }
+  }
+  // res contains the correct result
+  // apply correction if not rounding to nearest
+  if (rnd_mode != ROUNDING_TO_NEAREST) {
+    rounding_correction (rnd_mode,
+        		 is_inexact_lt_midpoint, is_inexact_gt_midpoint,
+        		 is_midpoint_lt_even, is_midpoint_gt_even,
+        		 e4, &res, ptrfpsf);
+  }
+  if (is_midpoint_lt_even || is_midpoint_gt_even ||
+      is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
+    // set the inexact flag
+    *ptrfpsf |= INEXACT_EXCEPTION;
+    if (is_tiny)
+      *ptrfpsf |= UNDERFLOW_EXCEPTION;
+  }
+
+  *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+  *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+  *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+  *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+  *ptrres = res;
+  return;
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+static void
+bid128_ext_fma (int *ptr_is_midpoint_lt_even,
+        	int *ptr_is_midpoint_gt_even,
+        	int *ptr_is_inexact_lt_midpoint,
+        	int *ptr_is_inexact_gt_midpoint, UINT128 * pres,
+        	UINT128 * px, UINT128 * py,
+        	UINT128 *
+        	pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+        	_EXC_INFO_PARAM) {
+  UINT128 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+static UINT128
+bid128_ext_fma (int *ptr_is_midpoint_lt_even,
+        	int *ptr_is_midpoint_gt_even,
+        	int *ptr_is_inexact_lt_midpoint,
+        	int *ptr_is_inexact_gt_midpoint, UINT128 x, UINT128 y,
+        	UINT128 z _RND_MODE_PARAM _EXC_FLAGS_PARAM
+        	_EXC_MASKS_PARAM _EXC_INFO_PARAM) {
+#endif
+
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT64 x_sign, y_sign, z_sign, p_sign, tmp_sign;
+  UINT64 x_exp = 0, y_exp = 0, z_exp = 0, p_exp;
+  int true_p_exp;
+  UINT128 C1, C2, C3;
+  UINT256 C4;
+  int q1 = 0, q2 = 0, q3 = 0, q4;
+  int e1, e2, e3, e4;
+  int scale, ind, delta, x0;
+  int p34 = P34; // used to modify the limit on the number of digits
+  BID_UI64DOUBLE tmp;
+  int x_nr_bits, y_nr_bits, z_nr_bits;
+  unsigned int save_fpsf;
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0;
+  int is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0;
+  int is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0;
+  int incr_exp = 0;
+  int lsb;
+  int lt_half_ulp = 0;
+  int eq_half_ulp = 0;
+  int gt_half_ulp = 0;
+  int is_tiny = 0;
+  UINT64 R64, tmp64;
+  UINT128 P128, R128;
+  UINT192 P192, R192;
+  UINT256 R256;
+
+  // the following are based on the table of special cases for fma; the NaN
+  // behavior is similar to that of the IA-64 Architecture fma 
+
+  // identify cases where at least one operand is NaN
+
+  BID_SWAP128 (x);
+  BID_SWAP128 (y);
+  BID_SWAP128 (z);
+  if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
+    // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)
+    // check first for non-canonical NaN payload
+    if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
+        (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
+         (y.w[0] > 0x38c15b09ffffffffull))) {
+      y.w[1] = y.w[1] & 0xffffc00000000000ull;
+      y.w[0] = 0x0ull;
+    }
+    if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
+      // set invalid flag
+      *pfpsf |= INVALID_EXCEPTION;
+      // return quiet (y)
+      res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
+      res.w[0] = y.w[0];
+    } else { // y is QNaN
+      // return y
+      res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
+      res.w[0] = y.w[0];
+      // if z = SNaN or x = SNaN signal invalid exception
+      if ((z.w[1] & MASK_SNAN) == MASK_SNAN ||
+          (x.w[1] & MASK_SNAN) == MASK_SNAN) {
+        // set invalid flag
+        *pfpsf |= INVALID_EXCEPTION;
+      }
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } else if ((z.w[1] & MASK_NAN) == MASK_NAN) { // z is NAN
+    // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)
+    // check first for non-canonical NaN payload
+    if (((z.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
+        (((z.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
+         (z.w[0] > 0x38c15b09ffffffffull))) {
+      z.w[1] = z.w[1] & 0xffffc00000000000ull;
+      z.w[0] = 0x0ull;
+    }
+    if ((z.w[1] & MASK_SNAN) == MASK_SNAN) { // z is SNAN 
+      // set invalid flag 
+      *pfpsf |= INVALID_EXCEPTION;
+      // return quiet (z) 
+      res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
+      res.w[0] = z.w[0];
+    } else { // z is QNaN 
+      // return z  
+      res.w[1] = z.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
+      res.w[0] = z.w[0];
+      // if x = SNaN signal invalid exception
+      if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {
+        // set invalid flag
+        *pfpsf |= INVALID_EXCEPTION;
+      }
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } else if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
+    // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)
+    // check first for non-canonical NaN payload
+    if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
+        (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
+         (x.w[0] > 0x38c15b09ffffffffull))) {
+      x.w[1] = x.w[1] & 0xffffc00000000000ull;
+      x.w[0] = 0x0ull;
+    }
+    if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN 
+      // set invalid flag 
+      *pfpsf |= INVALID_EXCEPTION;
+      // return quiet (x) 
+      res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
+      res.w[0] = x.w[0];
+    } else { // x is QNaN 
+      // return x  
+      res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
+      res.w[0] = x.w[0];
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  }
+  // x, y, z are 0, f, or inf but not NaN => unpack the arguments and check
+  // for non-canonical values
+
+  x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
+  C1.w[1] = x.w[1] & MASK_COEFF;
+  C1.w[0] = x.w[0];
+  if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
+    // if x is not infinity check for non-canonical values - treated as zero
+    if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
+      // non-canonical
+      x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
+      C1.w[1] = 0; // significand high
+      C1.w[0] = 0; // significand low
+    } else { // G0_G1 != 11
+      x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
+      if (C1.w[1] > 0x0001ed09bead87c0ull ||
+          (C1.w[1] == 0x0001ed09bead87c0ull &&
+           C1.w[0] > 0x378d8e63ffffffffull)) {
+        // x is non-canonical if coefficient is larger than 10^34 -1
+        C1.w[1] = 0;
+        C1.w[0] = 0;
+      } else { // canonical          
+        ;
+      }
+    }
+  }
+  y_sign = y.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
+  C2.w[1] = y.w[1] & MASK_COEFF;
+  C2.w[0] = y.w[0];
+  if ((y.w[1] & MASK_ANY_INF) != MASK_INF) { // y != inf
+    // if y is not infinity check for non-canonical values - treated as zero
+    if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
+      // non-canonical
+      y_exp = (y.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
+      C2.w[1] = 0; // significand high
+      C2.w[0] = 0; // significand low 
+    } else { // G0_G1 != 11
+      y_exp = y.w[1] & MASK_EXP; // biased and shifted left 49 bits
+      if (C2.w[1] > 0x0001ed09bead87c0ull ||
+          (C2.w[1] == 0x0001ed09bead87c0ull &&
+           C2.w[0] > 0x378d8e63ffffffffull)) {
+        // y is non-canonical if coefficient is larger than 10^34 -1
+        C2.w[1] = 0;
+        C2.w[0] = 0;
+      } else { // canonical
+        ;
+      }
+    }
+  }
+  z_sign = z.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
+  C3.w[1] = z.w[1] & MASK_COEFF;
+  C3.w[0] = z.w[0];
+  if ((z.w[1] & MASK_ANY_INF) != MASK_INF) { // z != inf
+    // if z is not infinity check for non-canonical values - treated as zero
+    if ((z.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
+      // non-canonical
+      z_exp = (z.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
+      C3.w[1] = 0; // significand high
+      C3.w[0] = 0; // significand low 
+    } else { // G0_G1 != 11
+      z_exp = z.w[1] & MASK_EXP; // biased and shifted left 49 bits
+      if (C3.w[1] > 0x0001ed09bead87c0ull ||
+          (C3.w[1] == 0x0001ed09bead87c0ull &&
+           C3.w[0] > 0x378d8e63ffffffffull)) {
+        // z is non-canonical if coefficient is larger than 10^34 -1
+        C3.w[1] = 0;
+        C3.w[0] = 0;
+      } else { // canonical
+        ;
+      }
+    }
+  }
+
+  p_sign = x_sign ^ y_sign; // sign of the product
+
+  // identify cases where at least one operand is infinity
+
+  if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
+    if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
+      if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
+        if (p_sign == z_sign) {
+          res.w[1] = z_sign | MASK_INF;
+          res.w[0] = 0x0;
+        } else {
+          // return QNaN Indefinite
+          res.w[1] = 0x7c00000000000000ull;
+          res.w[0] = 0x0000000000000000ull;
+          // set invalid flag
+          *pfpsf |= INVALID_EXCEPTION;
+        }
+      } else { // z = 0 or z = f
+        res.w[1] = p_sign | MASK_INF;
+        res.w[0] = 0x0;
+      }
+    } else if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f
+      if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
+        if (p_sign == z_sign) {
+          res.w[1] = z_sign | MASK_INF;
+          res.w[0] = 0x0;
+        } else {
+          // return QNaN Indefinite 
+          res.w[1] = 0x7c00000000000000ull;
+          res.w[0] = 0x0000000000000000ull;
+          // set invalid flag
+          *pfpsf |= INVALID_EXCEPTION;
+        }
+      } else { // z = 0 or z = f
+        res.w[1] = p_sign | MASK_INF;
+        res.w[0] = 0x0;
+      }
+    } else { // y = 0
+      // return QNaN Indefinite
+      res.w[1] = 0x7c00000000000000ull;
+      res.w[0] = 0x0000000000000000ull;
+      // set invalid flag
+      *pfpsf |= INVALID_EXCEPTION;
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
+    if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
+      // x = f, necessarily
+      if ((p_sign != z_sign)
+          || (C1.w[1] == 0x0ull && C1.w[0] == 0x0ull)) {
+        // return QNaN Indefinite
+        res.w[1] = 0x7c00000000000000ull;
+        res.w[0] = 0x0000000000000000ull;
+        // set invalid flag
+        *pfpsf |= INVALID_EXCEPTION;
+      } else {
+        res.w[1] = z_sign | MASK_INF;
+        res.w[0] = 0x0;
+      }
+    } else if (C1.w[1] == 0x0 && C1.w[0] == 0x0) { // x = 0
+      // z = 0, f, inf
+      // return QNaN Indefinite
+      res.w[1] = 0x7c00000000000000ull;
+      res.w[0] = 0x0000000000000000ull;
+      // set invalid flag
+      *pfpsf |= INVALID_EXCEPTION;
+    } else {
+      // x = f and z = 0, f, necessarily
+      res.w[1] = p_sign | MASK_INF;
+      res.w[0] = 0x0;
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } else if ((z.w[1] & MASK_ANY_INF) == MASK_INF) { // z = inf
+    // x = 0, f and y = 0, f, necessarily
+    res.w[1] = z_sign | MASK_INF;
+    res.w[0] = 0x0;
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  }
+
+  true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
+  if (true_p_exp < -6176)
+    p_exp = 0; // cannot be less than EXP_MIN
+  else
+    p_exp = (UINT64) (true_p_exp + 6176) << 49;
+
+  if (((C1.w[1] == 0x0 && C1.w[0] == 0x0) || (C2.w[1] == 0x0 && C2.w[0] == 0x0)) && C3.w[1] == 0x0 && C3.w[0] == 0x0) { // (x = 0 or y = 0) and z = 0
+    // the result is 0
+    if (p_exp < z_exp)
+      res.w[1] = p_exp; // preferred exponent
+    else
+      res.w[1] = z_exp; // preferred exponent
+    if (p_sign == z_sign) {
+      res.w[1] |= z_sign;
+      res.w[0] = 0x0;
+    } else { // x * y and z have opposite signs
+      if (rnd_mode == ROUNDING_DOWN) {
+        // res = -0.0
+        res.w[1] |= MASK_SIGN;
+        res.w[0] = 0x0;
+      } else {
+        // res = +0.0
+        // res.w[1] |= 0x0;
+        res.w[0] = 0x0;
+      }
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  }
+  // from this point on, we may need to know the number of decimal digits
+  // in the significands of x, y, z when x, y, z != 0
+
+  if (C1.w[1] != 0 || C1.w[0] != 0) { // x = f (non-zero finite)
+    // q1 = nr. of decimal digits in x
+    // determine first the nr. of bits in x
+    if (C1.w[1] == 0) {
+      if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
+        // split the 64-bit value in two 32-bit halves to avoid rounding errors
+        if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
+          tmp.d = (double) (C1.w[0] >> 32); // exact conversion
+          x_nr_bits =
+            33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        } else { // x < 2^32
+          tmp.d = (double) (C1.w[0]); // exact conversion
+          x_nr_bits =
+            1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        }
+      } else { // if x < 2^53
+        tmp.d = (double) C1.w[0]; // exact conversion
+        x_nr_bits =
+          1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+      }
+    } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
+      tmp.d = (double) C1.w[1]; // exact conversion
+      x_nr_bits =
+        65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+    }
+    q1 = nr_digits[x_nr_bits - 1].digits;
+    if (q1 == 0) {
+      q1 = nr_digits[x_nr_bits - 1].digits1;
+      if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
+          (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
+           C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
+        q1++;
+    }
+  }
+
+  if (C2.w[1] != 0 || C2.w[0] != 0) { // y = f (non-zero finite)
+    if (C2.w[1] == 0) {
+      if (C2.w[0] >= 0x0020000000000000ull) { // y >= 2^53
+        // split the 64-bit value in two 32-bit halves to avoid rounding errors
+        if (C2.w[0] >= 0x0000000100000000ull) { // y >= 2^32
+          tmp.d = (double) (C2.w[0] >> 32); // exact conversion
+          y_nr_bits =
+            32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        } else { // y < 2^32
+          tmp.d = (double) C2.w[0]; // exact conversion
+          y_nr_bits =
+            ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        }
+      } else { // if y < 2^53
+        tmp.d = (double) C2.w[0]; // exact conversion
+        y_nr_bits =
+          ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+      }
+    } else { // C2.w[1] != 0 => nr. bits = 64 + nr_bits (C2.w[1])
+      tmp.d = (double) C2.w[1]; // exact conversion
+      y_nr_bits =
+        64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+    }
+
+    q2 = nr_digits[y_nr_bits].digits;
+    if (q2 == 0) {
+      q2 = nr_digits[y_nr_bits].digits1;
+      if (C2.w[1] > nr_digits[y_nr_bits].threshold_hi ||
+          (C2.w[1] == nr_digits[y_nr_bits].threshold_hi &&
+           C2.w[0] >= nr_digits[y_nr_bits].threshold_lo))
+        q2++;
+    }
+  }
+
+  if (C3.w[1] != 0 || C3.w[0] != 0) { // z = f (non-zero finite)
+    if (C3.w[1] == 0) {
+      if (C3.w[0] >= 0x0020000000000000ull) { // z >= 2^53
+        // split the 64-bit value in two 32-bit halves to avoid rounding errors
+        if (C3.w[0] >= 0x0000000100000000ull) { // z >= 2^32
+          tmp.d = (double) (C3.w[0] >> 32); // exact conversion
+          z_nr_bits =
+            32 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        } else { // z < 2^32
+          tmp.d = (double) C3.w[0]; // exact conversion
+          z_nr_bits =
+            ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+        }
+      } else { // if z < 2^53
+        tmp.d = (double) C3.w[0]; // exact conversion
+        z_nr_bits =
+          ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+      }
+    } else { // C3.w[1] != 0 => nr. bits = 64 + nr_bits (C3.w[1])
+      tmp.d = (double) C3.w[1]; // exact conversion
+      z_nr_bits =
+        64 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+    }
+
+    q3 = nr_digits[z_nr_bits].digits;
+    if (q3 == 0) {
+      q3 = nr_digits[z_nr_bits].digits1;
+      if (C3.w[1] > nr_digits[z_nr_bits].threshold_hi ||
+          (C3.w[1] == nr_digits[z_nr_bits].threshold_hi &&
+           C3.w[0] >= nr_digits[z_nr_bits].threshold_lo))
+        q3++;
+    }
+  }
+
+  if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
+      (C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
+    // x = 0 or y = 0
+    // z = f, necessarily; for 0 + z return z, with the preferred exponent
+    // the result is z, but need to get the preferred exponent
+    if (z_exp <= p_exp) { // the preferred exponent is z_exp
+      res.w[1] = z_sign | (z_exp & MASK_EXP) | C3.w[1];
+      res.w[0] = C3.w[0];
+    } else { // if (p_exp < z_exp) the preferred exponent is p_exp
+      // return (C3 * 10^scale) * 10^(z_exp - scale)
+      // where scale = min (p34-q3, (z_exp-p_exp) >> 49)
+      scale = p34 - q3;
+      ind = (z_exp - p_exp) >> 49;
+      if (ind < scale)
+        scale = ind;
+      if (scale == 0) {
+        res.w[1] = z.w[1]; // & MASK_COEFF, which is redundant
+        res.w[0] = z.w[0];
+      } else if (q3 <= 19) { // z fits in 64 bits 
+        if (scale <= 19) { // 10^scale fits in 64 bits
+          // 64 x 64 C3.w[0] * ten2k64[scale]
+          __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+        } else { // 10^scale fits in 128 bits
+          // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+          __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
+        }
+      } else { // z fits in 128 bits, but 10^scale must fit in 64 bits 
+        // 64 x 128 ten2k64[scale] * C3
+        __mul_128x64_to_128 (res, ten2k64[scale], C3);
+      }
+      // subtract scale from the exponent
+      z_exp = z_exp - ((UINT64) scale << 49);
+      res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+    }
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } else {
+    ; // continue with x = f, y = f, z = 0 or x = f, y = f, z = f
+  }
+
+  e1 = (x_exp >> 49) - 6176; // unbiased exponent of x 
+  e2 = (y_exp >> 49) - 6176; // unbiased exponent of y 
+  e3 = (z_exp >> 49) - 6176; // unbiased exponent of z
+  e4 = e1 + e2; // unbiased exponent of the exact x * y
+
+  // calculate C1 * C2 and its number of decimal digits, q4
+
+  // the exact product has either q1 + q2 - 1 or q1 + q2 decimal digits
+  // where 2 <= q1 + q2 <= 68
+  // calculate C4 = C1 * C2 and determine q
+  C4.w[3] = C4.w[2] = C4.w[1] = C4.w[0] = 0;
+  if (q1 + q2 <= 19) { // if 2 <= q1 + q2 <= 19, C4 = C1 * C2 fits in 64 bits
+    C4.w[0] = C1.w[0] * C2.w[0];
+    // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
+    if (C4.w[0] < ten2k64[q1 + q2 - 1])
+      q4 = q1 + q2 - 1; // q4 in [1, 18]
+    else
+      q4 = q1 + q2; // q4 in [2, 19]
+    // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
+  } else if (q1 + q2 == 20) { // C4 = C1 * C2 fits in 64 or 128 bits
+    // q1 <= 19 and q2 <= 19 so both C1 and C2 fit in 64 bits
+    __mul_64x64_to_128MACH (C4, C1.w[0], C2.w[0]);
+    // if C4 < 10^(q1+q2-1) = 10^19 then q4 = q1+q2-1 = 19 else q4 = q1+q2 = 20
+    if (C4.w[1] == 0 && C4.w[0] < ten2k64[19]) { // 19 = q1+q2-1
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
+      q4 = 19; // 19 = q1 + q2 - 1
+    } else {
+      // if (C4.w[1] == 0)
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      q4 = 20; // 20 = q1 + q2
+    }
+  } else if (q1 + q2 <= 38) { // 21 <= q1 + q2 <= 38
+    // C4 = C1 * C2 fits in 64 or 128 bits
+    // (64 bits possibly, but only when q1 + q2 = 21 and C4 has 20 digits)
+    // at least one of C1, C2 has at most 19 decimal digits & fits in 64 bits
+    if (q1 <= 19) {
+      __mul_128x64_to_128 (C4, C1.w[0], C2);
+    } else { // q2 <= 19
+      __mul_128x64_to_128 (C4, C2.w[0], C1);
+    }
+    // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
+    if (C4.w[1] < ten2k128[q1 + q2 - 21].w[1] ||
+        (C4.w[1] == ten2k128[q1 + q2 - 21].w[1] &&
+         C4.w[0] < ten2k128[q1 + q2 - 21].w[0])) {
+      // if (C4.w[1] == 0) // q4 = 20, necessarily
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 64;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      q4 = q1 + q2 - 1; // q4 in [20, 37]
+    } else {
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      q4 = q1 + q2; // q4 in [21, 38]
+    }
+  } else if (q1 + q2 == 39) { // C4 = C1 * C2 fits in 128 or 192 bits
+    // both C1 and C2 fit in 128 bits (actually in 113 bits)
+    // may replace this by 128x128_to192
+    __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] is 0
+    // if C4 < 10^(q1+q2-1) = 10^38 then q4 = q1+q2-1 = 38 else q4 = q1+q2 = 39
+    if (C4.w[2] == 0 && (C4.w[1] < ten2k128[18].w[1] ||
+        		 (C4.w[1] == ten2k128[18].w[1]
+        		  && C4.w[0] < ten2k128[18].w[0]))) {
+      // 18 = 38 - 20 = q1+q2-1 - 20
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      q4 = 38; // 38 = q1 + q2 - 1
+    } else {
+      // if (C4.w[2] == 0)
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      q4 = 39; // 39 = q1 + q2
+    }
+  } else if (q1 + q2 <= 57) { // 40 <= q1 + q2 <= 57
+    // C4 = C1 * C2 fits in 128 or 192 bits
+    // (128 bits possibly, but only when q1 + q2 = 40 and C4 has 39 digits)
+    // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
+    // may fit in 64 bits
+    if (C1.w[1] == 0) { // C1 fits in 64 bits
+      // __mul_64x128_full (REShi64, RESlo128, A64, B128)
+      __mul_64x128_full (C4.w[2], C4, C1.w[0], C2);
+    } else if (C2.w[1] == 0) { // C2 fits in 64 bits
+      // __mul_64x128_full (REShi64, RESlo128, A64, B128)
+      __mul_64x128_full (C4.w[2], C4, C2.w[0], C1);
+    } else { // both C1 and C2 require 128 bits
+      // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1);
+      __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
+    }
+    // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
+    if (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
+        (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
+         (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
+          (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
+           C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))) {
+      // if (C4.w[2] == 0) // q4 = 39, necessarily
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 128;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      q4 = q1 + q2 - 1; // q4 in [39, 56]
+    } else {
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      q4 = q1 + q2; // q4 in [40, 57]
+    }
+  } else if (q1 + q2 == 58) { // C4 = C1 * C2 fits in 192 or 256 bits
+    // both C1 and C2 fit in 128 bits (actually in 113 bits); at most one
+    // may fit in 64 bits
+    if (C1.w[1] == 0) { // C1 * C2 will fit in 192 bits
+      __mul_64x128_full (C4.w[2], C4, C1.w[0], C2); // may use 64x128_to_192
+    } else if (C2.w[1] == 0) { // C1 * C2 will fit in 192 bits
+      __mul_64x128_full (C4.w[2], C4, C2.w[0], C1); // may use 64x128_to_192
+    } else { // C1 * C2 will fit in 192 bits or in 256 bits
+      __mul_128x128_to_256 (C4, C1, C2);
+    }
+    // if C4 < 10^(q1+q2-1) = 10^57 then q4 = q1+q2-1 = 57 else q4 = q1+q2 = 58
+    if (C4.w[3] == 0 && (C4.w[2] < ten2k256[18].w[2] ||
+        		 (C4.w[2] == ten2k256[18].w[2]
+        		  && (C4.w[1] < ten2k256[18].w[1]
+        		      || (C4.w[1] == ten2k256[18].w[1]
+        			  && C4.w[0] < ten2k256[18].w[0]))))) {
+      // 18 = 57 - 39 = q1+q2-1 - 39
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      q4 = 57; // 57 = q1 + q2 - 1
+    } else {
+      // if (C4.w[3] == 0)
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
+      q4 = 58; // 58 = q1 + q2
+    }
+  } else { // if 59 <= q1 + q2 <= 68
+    // C4 = C1 * C2 fits in 192 or 256 bits
+    // (192 bits possibly, but only when q1 + q2 = 59 and C4 has 58 digits)
+    // both C1 and C2 fit in 128 bits (actually in 113 bits); none fits in
+    // 64 bits
+    // may use __mul_128x128_to_192 (C4.w[2], C4.w[0], C2.w[0], C1);
+    __mul_128x128_to_256 (C4, C1, C2); // C4.w[3] = 0
+    // if C4 < 10^(q1+q2-1) then q4 = q1 + q2 - 1 else q4 = q1 + q2
+    if (C4.w[3] < ten2k256[q1 + q2 - 40].w[3] ||
+        (C4.w[3] == ten2k256[q1 + q2 - 40].w[3] &&
+         (C4.w[2] < ten2k256[q1 + q2 - 40].w[2] ||
+          (C4.w[2] == ten2k256[q1 + q2 - 40].w[2] &&
+           (C4.w[1] < ten2k256[q1 + q2 - 40].w[1] ||
+            (C4.w[1] == ten2k256[q1 + q2 - 40].w[1] &&
+             C4.w[0] < ten2k256[q1 + q2 - 40].w[0])))))) {
+      // if (C4.w[3] == 0) // q4 = 58, necessarily
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 192;
+      // else
+      //   length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
+      q4 = q1 + q2 - 1; // q4 in [58, 67]
+    } else {
+      // length of C1 * C2 rounded up to a multiple of 64 bits is len = 256;
+      q4 = q1 + q2; // q4 in [59, 68]
+    }
+  }
+
+  if (C3.w[1] == 0x0 && C3.w[0] == 0x0) { // x = f, y = f, z = 0
+    save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
+    *pfpsf = 0;
+
+    if (q4 > p34) {
+
+      // truncate C4 to p34 digits into res
+      // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68
+      x0 = q4 - p34;
+      if (q4 <= 38) {
+        P128.w[1] = C4.w[1];
+        P128.w[0] = C4.w[0];
+        round128_19_38 (q4, x0, P128, &res, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+      } else if (q4 <= 57) { // 35 <= q4 <= 57
+        P192.w[2] = C4.w[2];
+        P192.w[1] = C4.w[1];
+        P192.w[0] = C4.w[0];
+        round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+        res.w[0] = R192.w[0];
+        res.w[1] = R192.w[1];
+      } else { // if (q4 <= 68)
+        round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+        res.w[0] = R256.w[0];
+        res.w[1] = R256.w[1];
+      }
+      e4 = e4 + x0;
+      if (incr_exp) {
+        e4 = e4 + 1;
+      }
+      q4 = p34;
+      // res is now the coefficient of the result rounded to the destination 
+      // precision, with unbounded exponent; the exponent is e4; q4=digits(res)
+    } else { // if (q4 <= p34)
+      // C4 * 10^e4 is the result rounded to the destination precision, with  
+      // unbounded exponent (which is exact)
+
+      if ((q4 + e4 <= p34 + expmax) && (e4 > expmax)) {
+        // e4 is too large, but can be brought within range by scaling up C4
+        scale = e4 - expmax; // 1 <= scale < P-q4 <= P-1 => 1 <= scale <= P-2
+        // res = (C4 * 10^scale) * 10^expmax
+        if (q4 <= 19) { // C4 fits in 64 bits
+          if (scale <= 19) { // 10^scale fits in 64 bits
+            // 64 x 64 C4.w[0] * ten2k64[scale]
+            __mul_64x64_to_128MACH (res, C4.w[0], ten2k64[scale]);
+          } else { // 10^scale fits in 128 bits
+            // 64 x 128 C4.w[0] * ten2k128[scale - 20]
+            __mul_128x64_to_128 (res, C4.w[0], ten2k128[scale - 20]);
+          }
+        } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
+          // 64 x 128 ten2k64[scale] * CC43
+          __mul_128x64_to_128 (res, ten2k64[scale], C4);
+        }
+        e4 = e4 - scale; // expmax
+        q4 = q4 + scale;
+      } else {
+        res.w[1] = C4.w[1];
+        res.w[0] = C4.w[0];
+      }
+      // res is the coefficient of the result rounded to the destination 
+      // precision, with unbounded exponent (it has q4 digits); the exponent 
+      // is e4 (exact result)
+    }
+
+    // check for overflow
+    if (q4 + e4 > p34 + expmax) {
+      if (rnd_mode == ROUNDING_TO_NEAREST) {
+        res.w[1] = p_sign | 0x7800000000000000ull; // +/-inf
+        res.w[0] = 0x0000000000000000ull;
+        *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+      } else {
+        res.w[1] = p_sign | res.w[1];
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e4, &res, pfpsf);
+      }
+      *pfpsf |= save_fpsf;
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+    }
+    // check for underflow
+    if (q4 + e4 < expmin + P34) {
+      is_tiny = 1; // the result is tiny
+      if (e4 < expmin) {
+        // if e4 < expmin, we must truncate more of res
+        x0 = expmin - e4; // x0 >= 1
+        is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+        is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+        is_midpoint_lt_even0 = is_midpoint_lt_even;
+        is_midpoint_gt_even0 = is_midpoint_gt_even;
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+        // the number of decimal digits in res is q4
+        if (x0 < q4) { // 1 <= x0 <= q4-1 => round res to q4 - x0 digits
+          if (q4 <= 18) { // 2 <= q4 <= 18, 1 <= x0 <= 17
+            round64_2_18 (q4, x0, res.w[0], &R64, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+            if (incr_exp) {
+              // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
+              R64 = ten2k64[q4 - x0];
+            }
+            // res.w[1] = 0; (from above)
+            res.w[0] = R64;
+          } else { // if (q4 <= 34)
+            // 19 <= q4 <= 38
+            P128.w[1] = res.w[1];
+            P128.w[0] = res.w[0];
+            round128_19_38 (q4, x0, P128, &res, &incr_exp,
+        		    &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		    &is_inexact_lt_midpoint,
+        		    &is_inexact_gt_midpoint);
+            if (incr_exp) {
+              // increase coefficient by a factor of 10; this will be <= 10^33
+              // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
+              if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
+        	// res.w[1] = 0;
+        	res.w[0] = ten2k64[q4 - x0];
+              } else { // 20 <= q4 - x0 <= 37
+        	res.w[0] = ten2k128[q4 - x0 - 20].w[0];
+        	res.w[1] = ten2k128[q4 - x0 - 20].w[1];
+              }
+            }
+          }
+          e4 = e4 + x0; // expmin 
+        } else if (x0 == q4) {
+          // the second rounding is for 0.d(0)d(1)...d(q4-1) * 10^emin
+          // determine relationship with 1/2 ulp
+          if (q4 <= 19) {
+            if (res.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
+              lt_half_ulp = 1;
+              is_inexact_lt_midpoint = 1;
+            } else if (res.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
+              eq_half_ulp = 1;
+              is_midpoint_gt_even = 1;
+            } else { // > 1/2 ulp
+              // gt_half_ulp = 1;
+              is_inexact_gt_midpoint = 1;
+            }
+          } else { // if (q4 <= 34)
+            if (res.w[1] < midpoint128[q4 - 20].w[1] || 
+                (res.w[1] == midpoint128[q4 - 20].w[1] && 
+                res.w[0] < midpoint128[q4 - 20].w[0])) { // < 1/2 ulp
+              lt_half_ulp = 1;
+              is_inexact_lt_midpoint = 1;
+            } else if (res.w[1] == midpoint128[q4 - 20].w[1] && 
+                res.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
+              eq_half_ulp = 1;
+              is_midpoint_gt_even = 1;
+            } else { // > 1/2 ulp
+              // gt_half_ulp = 1;
+              is_inexact_gt_midpoint = 1;
+            }
+          }
+          if (lt_half_ulp || eq_half_ulp) {
+            // res = +0.0 * 10^expmin
+            res.w[1] = 0x0000000000000000ull;
+            res.w[0] = 0x0000000000000000ull;
+          } else { // if (gt_half_ulp)
+            // res = +1 * 10^expmin
+            res.w[1] = 0x0000000000000000ull;
+            res.w[0] = 0x0000000000000001ull;
+          }
+          e4 = expmin;
+        } else { // if (x0 > q4)
+          // the second rounding is for 0.0...d(0)d(1)...d(q4-1) * 10^emin
+          res.w[1] = 0;
+          res.w[0] = 0;
+          e4 = expmin;
+          is_inexact_lt_midpoint = 1;
+        }
+        // avoid a double rounding error
+        if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+            is_midpoint_lt_even) { // double rounding error upward
+          // res = res - 1
+          res.w[0]--;
+          if (res.w[0] == 0xffffffffffffffffull)
+            res.w[1]--;
+          // Note: a double rounding error upward is not possible; for this
+          // the result after the first rounding would have to be 99...95
+          // (35 digits in all), possibly followed by a number of zeros; this
+          // not possible for f * f + 0
+          is_midpoint_lt_even = 0;
+          is_inexact_lt_midpoint = 1;
+        } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+            is_midpoint_gt_even) { // double rounding error downward
+          // res = res + 1
+          res.w[0]++;
+          if (res.w[0] == 0)
+            res.w[1]++;
+          is_midpoint_gt_even = 0;
+          is_inexact_gt_midpoint = 1;
+        } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+          // if this second rounding was exact the result may still be 
+          // inexact because of the first rounding
+          if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+            is_inexact_gt_midpoint = 1;
+          }
+          if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+            is_inexact_lt_midpoint = 1;
+          }
+        } else if (is_midpoint_gt_even &&
+        	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+          // pulled up to a midpoint
+          is_inexact_lt_midpoint = 1;
+          is_inexact_gt_midpoint = 0;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else if (is_midpoint_lt_even &&
+        	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+          // pulled down to a midpoint
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 1;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else {
+          ;
+        }
+      } else { // if e4 >= emin then q4 < P and the result is tiny and exact
+        if (e3 < e4) {
+          // if (e3 < e4) the preferred exponent is e3
+          // return (C4 * 10^scale) * 10^(e4 - scale)
+          // where scale = min (p34-q4, (e4 - e3))
+          scale = p34 - q4;
+          ind = e4 - e3;
+          if (ind < scale)
+            scale = ind;
+          if (scale == 0) {
+            ; // res and e4 are unchanged
+          } else if (q4 <= 19) { // C4 fits in 64 bits
+            if (scale <= 19) { // 10^scale fits in 64 bits
+              // 64 x 64 res.w[0] * ten2k64[scale]
+              __mul_64x64_to_128MACH (res, res.w[0], ten2k64[scale]);
+            } else { // 10^scale fits in 128 bits
+              // 64 x 128 res.w[0] * ten2k128[scale - 20]
+              __mul_128x64_to_128 (res, res.w[0], ten2k128[scale - 20]);
+            }
+          } else { // res fits in 128 bits, but 10^scale must fit in 64 bits
+            // 64 x 128 ten2k64[scale] * C3
+            __mul_128x64_to_128 (res, ten2k64[scale], res);
+          }
+          // subtract scale from the exponent
+          e4 = e4 - scale;
+        }
+      }
+
+      // check for inexact result
+      if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+          is_midpoint_lt_even || is_midpoint_gt_even) {
+        // set the inexact flag and the underflow flag
+        *pfpsf |= INEXACT_EXCEPTION;
+        *pfpsf |= UNDERFLOW_EXCEPTION;
+      }
+      res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
+      if (rnd_mode != ROUNDING_TO_NEAREST) {
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e4, &res, pfpsf);
+      }
+      *pfpsf |= save_fpsf;
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+    }
+    // no overflow, and no underflow for rounding to nearest
+    res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1];
+
+    if (rnd_mode != ROUNDING_TO_NEAREST) {
+      rounding_correction (rnd_mode,
+        		   is_inexact_lt_midpoint,
+        		   is_inexact_gt_midpoint,
+        		   is_midpoint_lt_even, is_midpoint_gt_even,
+        		   e4, &res, pfpsf);
+      // if e4 = expmin && significand < 10^33 => result is tiny (for RD, RZ)
+      if (e4 == expmin) {
+        if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
+            ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
+             res.w[0] < 0x38c15b0a00000000ull)) {
+          is_tiny = 1;
+        }
+      }
+    }
+
+    if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+        is_midpoint_lt_even || is_midpoint_gt_even) {
+      // set the inexact flag
+      *pfpsf |= INEXACT_EXCEPTION;
+      if (is_tiny)
+        *pfpsf |= UNDERFLOW_EXCEPTION;
+    }
+
+    if ((*pfpsf & INEXACT_EXCEPTION) == 0) { // x * y is exact
+      // need to ensure that the result has the preferred exponent
+      p_exp = res.w[1] & MASK_EXP;
+      if (z_exp < p_exp) { // the preferred exponent is z_exp
+        // signficand of res in C3
+        C3.w[1] = res.w[1] & MASK_COEFF;
+        C3.w[0] = res.w[0];
+        // the number of decimal digits of x * y is q4 <= 34
+        // Note: the coefficient fits in 128 bits
+
+        // return (C3 * 10^scale) * 10^(p_exp - scale)
+        // where scale = min (p34-q4, (p_exp-z_exp) >> 49)
+        scale = p34 - q4;
+        ind = (p_exp - z_exp) >> 49;
+        if (ind < scale)
+          scale = ind;
+        // subtract scale from the exponent
+        p_exp = p_exp - ((UINT64) scale << 49);
+        if (scale == 0) {
+          ; // leave res unchanged
+        } else if (q4 <= 19) { // x * y fits in 64 bits
+          if (scale <= 19) { // 10^scale fits in 64 bits
+            // 64 x 64 C3.w[0] * ten2k64[scale] 
+            __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+          } else { // 10^scale fits in 128 bits 
+            // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+            __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
+          }
+          res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
+        } else { // x * y fits in 128 bits, but 10^scale must fit in 64 bits
+          // 64 x 128 ten2k64[scale] * C3 
+          __mul_128x64_to_128 (res, ten2k64[scale], C3);
+          res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
+        }
+      } // else leave the result as it is, because p_exp <= z_exp
+    }
+    *pfpsf |= save_fpsf;
+    *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+    *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+    *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+    *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+    BID_SWAP128 (res);
+    BID_RETURN (res)
+  } // else we have f * f + f
+
+  // continue with x = f, y = f, z = f
+
+  delta = q3 + e3 - q4 - e4;
+delta_ge_zero:
+  if (delta >= 0) {
+
+    if (p34 <= delta - 1 ||	// Case (1')
+        (p34 == delta && e3 + 6176 < p34 - q3)) { // Case (1''A)
+      // check for overflow, which can occur only in Case (1')
+      if ((q3 + e3) > (p34 + expmax) && p34 <= delta - 1) {
+        // e3 > expmax implies p34 <= delta-1 and e3 > expmax is a necessary
+        // condition for (q3 + e3) > (p34 + expmax)
+        if (rnd_mode == ROUNDING_TO_NEAREST) {
+          res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+          res.w[0] = 0x0000000000000000ull;
+          *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+        } else {
+          if (p_sign == z_sign) {
+            is_inexact_lt_midpoint = 1;
+          } else {
+            is_inexact_gt_midpoint = 1;
+          }
+          // q3 <= p34; if (q3 < p34) scale C3 up by 10^(p34-q3)
+          scale = p34 - q3;
+          if (scale == 0) {
+            res.w[1] = z_sign | C3.w[1];
+            res.w[0] = C3.w[0];
+          } else {
+            if (q3 <= 19) { // C3 fits in 64 bits
+              if (scale <= 19) { // 10^scale fits in 64 bits
+        	// 64 x 64 C3.w[0] * ten2k64[scale]
+        	__mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+              } else { // 10^scale fits in 128 bits
+        	// 64 x 128 C3.w[0] * ten2k128[scale - 20]
+        	__mul_128x64_to_128 (res, C3.w[0],
+        			     ten2k128[scale - 20]);
+              }
+            } else { // C3 fits in 128 bits, but 10^scale must fit in 64 bits
+              // 64 x 128 ten2k64[scale] * C3
+              __mul_128x64_to_128 (res, ten2k64[scale], C3);
+            }
+            // the coefficient in res has q3 + scale = p34 digits
+          }
+          e3 = e3 - scale;
+          res.w[1] = z_sign | res.w[1];
+          rounding_correction (rnd_mode,
+        		       is_inexact_lt_midpoint,
+        		       is_inexact_gt_midpoint,
+        		       is_midpoint_lt_even, is_midpoint_gt_even,
+        		       e3, &res, pfpsf);
+        }
+        *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+        *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+        *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+        *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+        BID_SWAP128 (res);
+        BID_RETURN (res)
+      }
+      // res = z
+      if (q3 < p34) { // the preferred exponent is z_exp - (p34 - q3)
+        // return (C3 * 10^scale) * 10^(z_exp - scale)
+        // where scale = min (p34-q3, z_exp-EMIN)
+        scale = p34 - q3;
+        ind = e3 + 6176;
+        if (ind < scale)
+          scale = ind;
+        if (scale == 0) {
+          res.w[1] = C3.w[1];
+          res.w[0] = C3.w[0];
+        } else if (q3 <= 19) { // z fits in 64 bits
+          if (scale <= 19) { // 10^scale fits in 64 bits
+            // 64 x 64 C3.w[0] * ten2k64[scale]
+            __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+          } else { // 10^scale fits in 128 bits
+            // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+            __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
+          }
+        } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
+          // 64 x 128 ten2k64[scale] * C3
+          __mul_128x64_to_128 (res, ten2k64[scale], C3);
+        }
+        // the coefficient in res has q3 + scale digits
+        // subtract scale from the exponent
+        z_exp = z_exp - ((UINT64) scale << 49);
+        e3 = e3 - scale;
+        res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+        if (scale + q3 < p34)
+          *pfpsf |= UNDERFLOW_EXCEPTION;
+      } else {
+        scale = 0;
+        res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | C3.w[1];
+        res.w[0] = C3.w[0];
+      }
+
+      // use the following to avoid double rounding errors when operating on
+      // mixed formats in rounding to nearest, and for correcting the result
+      // if not rounding to nearest
+      if ((p_sign != z_sign) && (delta == (q3 + scale + 1))) {
+        // there is a gap of exactly one digit between the scaled C3 and C4
+        // C3 * 10^ scale = 10^(q3+scale-1) <=> C3 = 10^(q3-1) is special case
+        if ((q3 <= 19 && C3.w[0] != ten2k64[q3 - 1]) ||
+            (q3 == 20 && (C3.w[1] != 0 || C3.w[0] != ten2k64[19])) ||
+            (q3 >= 21 && (C3.w[1] != ten2k128[q3 - 21].w[1] ||
+        		  C3.w[0] != ten2k128[q3 - 21].w[0]))) {
+          // C3 * 10^ scale != 10^(q3-1)
+          // if ((res.w[1] & MASK_COEFF) != 0x0000314dc6448d93ull ||
+          // res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33
+          is_inexact_gt_midpoint = 1; // if (z_sign), set as if for abs. value
+        } else { // if C3 * 10^scale = 10^(q3+scale-1)
+          // ok from above e3 = (z_exp >> 49) - 6176;
+          // the result is always inexact
+          if (q4 == 1) {
+            R64 = C4.w[0];
+          } else {
+            // if q4 > 1 then truncate C4 from q4 digits to 1 digit; 
+            // x = q4-1, 1 <= x <= 67 and check if this operation is exact
+            if (q4 <= 18) { // 2 <= q4 <= 18
+              round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
+        		    &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		    &is_inexact_lt_midpoint,
+        		    &is_inexact_gt_midpoint);
+            } else if (q4 <= 38) {
+              P128.w[1] = C4.w[1];
+              P128.w[0] = C4.w[0];
+              round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
+        		      &is_midpoint_lt_even,
+        		      &is_midpoint_gt_even,
+        		      &is_inexact_lt_midpoint,
+        		      &is_inexact_gt_midpoint);
+              R64 = R128.w[0]; // one decimal digit
+            } else if (q4 <= 57) {
+              P192.w[2] = C4.w[2];
+              P192.w[1] = C4.w[1];
+              P192.w[0] = C4.w[0];
+              round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
+        		      &is_midpoint_lt_even,
+        		      &is_midpoint_gt_even,
+        		      &is_inexact_lt_midpoint,
+        		      &is_inexact_gt_midpoint);
+              R64 = R192.w[0]; // one decimal digit
+            } else { // if (q4 <= 68)
+              round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
+        		      &is_midpoint_lt_even,
+        		      &is_midpoint_gt_even,
+        		      &is_inexact_lt_midpoint,
+        		      &is_inexact_gt_midpoint);
+              R64 = R256.w[0]; // one decimal digit
+            }
+            if (incr_exp) {
+              R64 = 10;
+            }
+          }
+          if (q4 == 1 && C4.w[0] == 5) {
+            is_inexact_lt_midpoint = 0;
+            is_inexact_gt_midpoint = 0;
+            is_midpoint_lt_even = 1;
+            is_midpoint_gt_even = 0;
+          } else if ((e3 == expmin) ||
+        	     R64 < 5 || (R64 == 5 && is_inexact_gt_midpoint)) {
+            // result does not change
+            is_inexact_lt_midpoint = 0;
+            is_inexact_gt_midpoint = 1;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+          } else {
+            is_inexact_lt_midpoint = 1;
+            is_inexact_gt_midpoint = 0;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+            // result decremented is 10^(q3+scale) - 1
+            if ((q3 + scale) <= 19) {
+              res.w[1] = 0;
+              res.w[0] = ten2k64[q3 + scale];
+            } else { // if ((q3 + scale + 1) <= 35)
+              res.w[1] = ten2k128[q3 + scale - 20].w[1];
+              res.w[0] = ten2k128[q3 + scale - 20].w[0];
+            }
+            res.w[0] = res.w[0] - 1; // borrow never occurs
+            z_exp = z_exp - EXP_P1;
+            e3 = e3 - 1;
+            res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
+          }
+          if (e3 == expmin) {
+            if (R64 < 5 || (R64 == 5 && !is_inexact_lt_midpoint)) {
+              ; // result not tiny (in round-to-nearest mode)
+            } else {
+              *pfpsf |= UNDERFLOW_EXCEPTION;
+            }
+          }
+        } // end 10^(q3+scale-1)
+        // set the inexact flag
+        *pfpsf |= INEXACT_EXCEPTION;
+      } else {
+        if (p_sign == z_sign) {
+          // if (z_sign), set as if for absolute value
+          is_inexact_lt_midpoint = 1;
+        } else { // if (p_sign != z_sign)
+          // if (z_sign), set as if for absolute value
+          is_inexact_gt_midpoint = 1;
+        }
+        *pfpsf |= INEXACT_EXCEPTION;
+      }
+      // the result is always inexact => set the inexact flag
+      // Determine tininess:
+      //    if (exp > expmin)
+      //      the result is not tiny
+      //    else // if exp = emin
+      //      if (q3 + scale < p34)
+      //        the result is tiny
+      //      else // if (q3 + scale = p34)
+      //        if (C3 * 10^scale > 10^33)
+      //          the result is not tiny
+      //        else // if C3 * 10^scale = 10^33
+      //          if (xy * z > 0)
+      //            the result is not tiny
+      //          else // if (xy * z < 0)
+      //            if (z > 0)
+      //              if rnd_mode != RP
+      //                the result is tiny
+      //              else // if RP
+      //                the result is not tiny
+      //            else // if (z < 0)
+      //              if rnd_mode != RM
+      //                the result is tiny
+      //              else // if RM
+      //                the result is not tiny
+      //              endif
+      //            endif
+      //          endif
+      //        endif
+      //      endif
+      //    endif 
+      if ((e3 == expmin && (q3 + scale) < p34) || 
+          (e3 == expmin && (q3 + scale) == p34 && 
+          (res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&	// 10^33_high
+          res.w[0] == 0x38c15b0a00000000ull &&	// 10^33_low
+          z_sign != p_sign && ((!z_sign && rnd_mode != ROUNDING_UP) || 
+          (z_sign && rnd_mode != ROUNDING_DOWN)))) {
+        *pfpsf |= UNDERFLOW_EXCEPTION;
+      }
+      if (rnd_mode != ROUNDING_TO_NEAREST) {
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e3, &res, pfpsf);
+      }
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else if (p34 == delta) { // Case (1''B)
+
+      // because Case (1''A) was treated above, e3 + 6176 >= p34 - q3
+      // and C3 can be scaled up to p34 digits if needed
+
+      // scale C3 to p34 digits if needed
+      scale = p34 - q3; // 0 <= scale <= p34 - 1
+      if (scale == 0) {
+        res.w[1] = C3.w[1];
+        res.w[0] = C3.w[0];
+      } else if (q3 <= 19) { // z fits in 64 bits
+        if (scale <= 19) { // 10^scale fits in 64 bits
+          // 64 x 64 C3.w[0] * ten2k64[scale]
+          __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+        } else { // 10^scale fits in 128 bits
+          // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+          __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
+        }
+      } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
+        // 64 x 128 ten2k64[scale] * C3
+        __mul_128x64_to_128 (res, ten2k64[scale], C3);
+      }
+      // subtract scale from the exponent
+      z_exp = z_exp - ((UINT64) scale << 49);
+      e3 = e3 - scale;
+      // now z_sign, z_exp, and res correspond to a z scaled to p34 = 34 digits
+
+      // determine whether x * y is less than, equal to, or greater than 
+      // 1/2 ulp (z)
+      if (q4 <= 19) {
+        if (C4.w[0] < midpoint64[q4 - 1]) { // < 1/2 ulp
+          lt_half_ulp = 1;
+        } else if (C4.w[0] == midpoint64[q4 - 1]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+        } else { // > 1/2 ulp
+          gt_half_ulp = 1;
+        }
+      } else if (q4 <= 38) {
+        if (C4.w[2] == 0 && (C4.w[1] < midpoint128[q4 - 20].w[1] || 
+            (C4.w[1] == midpoint128[q4 - 20].w[1] && 
+            C4.w[0] < midpoint128[q4 - 20].w[0]))) { // < 1/2 ulp
+          lt_half_ulp = 1;
+        } else if (C4.w[2] == 0 && C4.w[1] == midpoint128[q4 - 20].w[1] && 
+            C4.w[0] == midpoint128[q4 - 20].w[0]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+        } else { // > 1/2 ulp
+          gt_half_ulp = 1;
+        }
+      } else if (q4 <= 58) {
+        if (C4.w[3] == 0 && (C4.w[2] < midpoint192[q4 - 39].w[2] || 
+            (C4.w[2] == midpoint192[q4 - 39].w[2] && 
+            C4.w[1] < midpoint192[q4 - 39].w[1]) || 
+            (C4.w[2] == midpoint192[q4 - 39].w[2] && 
+            C4.w[1] == midpoint192[q4 - 39].w[1] && 
+            C4.w[0] < midpoint192[q4 - 39].w[0]))) { // < 1/2 ulp
+          lt_half_ulp = 1;
+        } else if (C4.w[3] == 0 && C4.w[2] == midpoint192[q4 - 39].w[2] && 
+            C4.w[1] == midpoint192[q4 - 39].w[1] && 
+            C4.w[0] == midpoint192[q4 - 39].w[0]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+        } else { // > 1/2 ulp
+          gt_half_ulp = 1;
+        }
+      } else {
+        if (C4.w[3] < midpoint256[q4 - 59].w[3] || 
+            (C4.w[3] == midpoint256[q4 - 59].w[3] && 
+            C4.w[2] < midpoint256[q4 - 59].w[2]) || 
+            (C4.w[3] == midpoint256[q4 - 59].w[3] && 
+            C4.w[2] == midpoint256[q4 - 59].w[2] && 
+            C4.w[1] < midpoint256[q4 - 59].w[1]) || 
+            (C4.w[3] == midpoint256[q4 - 59].w[3] && 
+            C4.w[2] == midpoint256[q4 - 59].w[2] && 
+            C4.w[1] == midpoint256[q4 - 59].w[1] && 
+            C4.w[0] < midpoint256[q4 - 59].w[0])) { // < 1/2 ulp
+          lt_half_ulp = 1;
+        } else if (C4.w[3] == midpoint256[q4 - 59].w[3] && 
+            C4.w[2] == midpoint256[q4 - 59].w[2] && 
+            C4.w[1] == midpoint256[q4 - 59].w[1] && 
+            C4.w[0] == midpoint256[q4 - 59].w[0]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+        } else { // > 1/2 ulp
+          gt_half_ulp = 1;
+        }
+      }
+
+      if (p_sign == z_sign) {
+        if (lt_half_ulp) {
+          res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+          // use the following to avoid double rounding errors when operating on
+          // mixed formats in rounding to nearest
+          is_inexact_lt_midpoint = 1; // if (z_sign), as if for absolute value
+        } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) {
+          // add 1 ulp to the significand
+          res.w[0]++;
+          if (res.w[0] == 0x0ull)
+            res.w[1]++;
+          // check for rounding overflow, when coeff == 10^34
+          if ((res.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull && 
+              res.w[0] == 0x378d8e6400000000ull) { // coefficient = 10^34
+            e3 = e3 + 1;
+            // coeff = 10^33
+            z_exp = ((UINT64) (e3 + 6176) << 49) & MASK_EXP;
+            res.w[1] = 0x0000314dc6448d93ull;
+            res.w[0] = 0x38c15b0a00000000ull;
+          }
+          // end add 1 ulp
+          res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+          if (eq_half_ulp) {
+            is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
+          } else {
+            is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value
+          }
+        } else { // if (eq_half_ulp && !(res.w[0] & 0x01))
+          // leave unchanged 
+          res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+          is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value
+        }
+        // the result is always inexact, and never tiny
+        // set the inexact flag
+        *pfpsf |= INEXACT_EXCEPTION;
+        // check for overflow
+        if (e3 > expmax && rnd_mode == ROUNDING_TO_NEAREST) {
+          res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+          res.w[0] = 0x0000000000000000ull;
+          *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+          *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+          *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+          *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+          *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+          BID_SWAP128 (res);
+          BID_RETURN (res)
+        }
+        if (rnd_mode != ROUNDING_TO_NEAREST) {
+          rounding_correction (rnd_mode,
+        		       is_inexact_lt_midpoint,
+        		       is_inexact_gt_midpoint,
+        		       is_midpoint_lt_even, is_midpoint_gt_even,
+        		       e3, &res, pfpsf);
+          z_exp = res.w[1] & MASK_EXP;
+        }
+      } else { // if (p_sign != z_sign)
+        // consider two cases, because C3 * 10^scale = 10^33 is a special case
+        if (res.w[1] != 0x0000314dc6448d93ull || 
+            res.w[0] != 0x38c15b0a00000000ull) { // C3 * 10^scale != 10^33
+          if (lt_half_ulp) {
+            res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+            // use the following to avoid double rounding errors when operating
+            // on mixed formats in rounding to nearest
+            is_inexact_gt_midpoint = 1; // if (z_sign), as if for absolute value
+          } else if ((eq_half_ulp && (res.w[0] & 0x01)) || gt_half_ulp) {
+            // subtract 1 ulp from the significand
+            res.w[0]--;
+            if (res.w[0] == 0xffffffffffffffffull)
+              res.w[1]--;
+            res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+            if (eq_half_ulp) {
+              is_midpoint_gt_even = 1; // if (z_sign), as if for absolute value
+            } else {
+              is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value
+            }
+          } else { // if (eq_half_ulp && !(res.w[0] & 0x01))
+            // leave unchanged
+            res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+            is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
+          }
+          // the result is always inexact, and never tiny
+          // check for overflow for RN
+          if (e3 > expmax) {
+            if (rnd_mode == ROUNDING_TO_NEAREST) {
+              res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+              res.w[0] = 0x0000000000000000ull;
+              *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+            } else {
+              rounding_correction (rnd_mode,
+        			   is_inexact_lt_midpoint,
+        			   is_inexact_gt_midpoint,
+        			   is_midpoint_lt_even,
+        			   is_midpoint_gt_even, e3, &res,
+        			   pfpsf);
+            }
+            *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+            *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+            *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+            *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+            BID_SWAP128 (res);
+            BID_RETURN (res)
+          }
+          // set the inexact flag
+          *pfpsf |= INEXACT_EXCEPTION;
+          if (rnd_mode != ROUNDING_TO_NEAREST) {
+            rounding_correction (rnd_mode,
+        			 is_inexact_lt_midpoint,
+        			 is_inexact_gt_midpoint,
+        			 is_midpoint_lt_even,
+        			 is_midpoint_gt_even, e3, &res, pfpsf);
+          }
+          z_exp = res.w[1] & MASK_EXP;
+        } else { // if C3 * 10^scale = 10^33
+          e3 = (z_exp >> 49) - 6176;
+          if (e3 > expmin) {
+            // the result is exact if exp > expmin and C4 = d*10^(q4-1), 
+            // where d = 1, 2, 3, ..., 9; it could be tiny too, but exact
+            if (q4 == 1) {
+              // if q4 = 1 the result is exact
+              // result coefficient = 10^34 - C4
+              res.w[1] = 0x0001ed09bead87c0ull;
+              res.w[0] = 0x378d8e6400000000ull - C4.w[0];
+              z_exp = z_exp - EXP_P1;
+              e3 = e3 - 1;
+              res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+            } else {
+              // if q4 > 1 then truncate C4 from q4 digits to 1 digit; 
+              // x = q4-1, 1 <= x <= 67 and check if this operation is exact
+              if (q4 <= 18) { // 2 <= q4 <= 18
+        	round64_2_18 (q4, q4 - 1, C4.w[0], &R64, &incr_exp,
+        		      &is_midpoint_lt_even,
+        		      &is_midpoint_gt_even,
+        		      &is_inexact_lt_midpoint,
+        		      &is_inexact_gt_midpoint);
+              } else if (q4 <= 38) {
+        	P128.w[1] = C4.w[1];
+        	P128.w[0] = C4.w[0];
+        	round128_19_38 (q4, q4 - 1, P128, &R128, &incr_exp,
+        			&is_midpoint_lt_even,
+        			&is_midpoint_gt_even,
+        			&is_inexact_lt_midpoint,
+        			&is_inexact_gt_midpoint);
+        	R64 = R128.w[0]; // one decimal digit
+              } else if (q4 <= 57) {
+        	P192.w[2] = C4.w[2];
+        	P192.w[1] = C4.w[1];
+        	P192.w[0] = C4.w[0];
+        	round192_39_57 (q4, q4 - 1, P192, &R192, &incr_exp,
+        			&is_midpoint_lt_even,
+        			&is_midpoint_gt_even,
+        			&is_inexact_lt_midpoint,
+        			&is_inexact_gt_midpoint);
+        	R64 = R192.w[0]; // one decimal digit
+              } else { // if (q4 <= 68)
+        	round256_58_76 (q4, q4 - 1, C4, &R256, &incr_exp,
+        			&is_midpoint_lt_even,
+        			&is_midpoint_gt_even,
+        			&is_inexact_lt_midpoint,
+        			&is_inexact_gt_midpoint);
+        	R64 = R256.w[0]; // one decimal digit
+              }
+              if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	  !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+        	// the result is exact: 10^34 - R64
+        	// incr_exp = 0 with certainty
+        	z_exp = z_exp - EXP_P1;
+        	e3 = e3 - 1;
+        	res.w[1] =
+        	  z_sign | (z_exp & MASK_EXP) | 0x0001ed09bead87c0ull;
+        	res.w[0] = 0x378d8e6400000000ull - R64;
+              } else {
+        	// We want R64 to be the top digit of C4, but we actually 
+        	// obtained (C4 * 10^(-q4+1))RN; a correction may be needed,
+        	// because the top digit is (C4 * 10^(-q4+1))RZ
+        	// however, if incr_exp = 1 then R64 = 10 with certainty
+        	if (incr_exp) {
+        	  R64 = 10;
+        	}
+        	// the result is inexact as C4 has more than 1 significant digit
+        	// and C3 * 10^scale = 10^33
+        	// example of case that is treated here:
+        	// 100...0 * 10^e3 - 0.41 * 10^e3 =
+        	// 0999...9.59 * 10^e3 -> rounds to 99...96*10^(e3-1)
+        	// note that (e3 > expmin}
+        	// in order to round, subtract R64 from 10^34 and then compare
+        	// C4 - R64 * 10^(q4-1) with 1/2 ulp
+        	// calculate 10^34 - R64
+        	res.w[1] = 0x0001ed09bead87c0ull;
+        	res.w[0] = 0x378d8e6400000000ull - R64;
+        	z_exp = z_exp - EXP_P1; // will be OR-ed with sign & significand
+        	// calculate C4 - R64 * 10^(q4-1); this is a rare case and
+        	// R64 is small, 1 <= R64 <= 9
+        	e3 = e3 - 1;
+        	if (is_inexact_lt_midpoint) {
+        	  is_inexact_lt_midpoint = 0;
+        	  is_inexact_gt_midpoint = 1;
+        	} else if (is_inexact_gt_midpoint) {
+        	  is_inexact_gt_midpoint = 0;
+        	  is_inexact_lt_midpoint = 1;
+        	} else if (is_midpoint_lt_even) {
+        	  is_midpoint_lt_even = 0;
+        	  is_midpoint_gt_even = 1;
+        	} else if (is_midpoint_gt_even) {
+        	  is_midpoint_gt_even = 0;
+        	  is_midpoint_lt_even = 1;
+        	} else {
+        	  ;
+        	}
+        	// the result is always inexact, and never tiny
+        	// check for overflow for RN
+        	if (e3 > expmax) {
+        	  if (rnd_mode == ROUNDING_TO_NEAREST) {
+        	    res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+        	    res.w[0] = 0x0000000000000000ull;
+        	    *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+        	  } else {
+        	    rounding_correction (rnd_mode,
+        				 is_inexact_lt_midpoint,
+        				 is_inexact_gt_midpoint,
+        				 is_midpoint_lt_even,
+        				 is_midpoint_gt_even, e3, &res,
+        				 pfpsf);
+        	  }
+        	  *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+        	  *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+        	  *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+        	  *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+        	  BID_SWAP128 (res);
+        	  BID_RETURN (res)
+        	}
+        	// set the inexact flag
+        	*pfpsf |= INEXACT_EXCEPTION;
+        	res.w[1] =
+        	  z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
+        	if (rnd_mode != ROUNDING_TO_NEAREST) {
+        	  rounding_correction (rnd_mode,
+        			       is_inexact_lt_midpoint,
+        			       is_inexact_gt_midpoint,
+        			       is_midpoint_lt_even,
+        			       is_midpoint_gt_even, e3, &res,
+        			       pfpsf);
+        	}
+        	z_exp = res.w[1] & MASK_EXP;
+              } // end result is inexact
+            } // end q4 > 1
+          } else { // if (e3 = emin)
+            // if e3 = expmin the result is also tiny (the condition for
+            // tininess is C4 > 050...0 [q4 digits] which is met because
+            // the msd of C4 is not zero)
+            // the result is tiny and inexact in all rounding modes;
+            // it is either 100...0 or 0999...9 (use lt_half_ulp, eq_half_ulp, 
+            // gt_half_ulp to calculate)
+            // if (lt_half_ulp || eq_half_ulp) res = 10^33 stays unchanged
+
+            // p_sign != z_sign so swap gt_half_ulp and lt_half_ulp
+            if (gt_half_ulp) { // res = 10^33 - 1
+              res.w[1] = 0x0000314dc6448d93ull;
+              res.w[0] = 0x38c15b09ffffffffull;
+            } else {
+              res.w[1] = 0x0000314dc6448d93ull;
+              res.w[0] = 0x38c15b0a00000000ull;
+            }
+            res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+            *pfpsf |= UNDERFLOW_EXCEPTION; // inexact is set later
+
+            if (eq_half_ulp) {
+              is_midpoint_lt_even = 1; // if (z_sign), as if for absolute value
+            } else if (lt_half_ulp) {
+              is_inexact_gt_midpoint = 1; //if(z_sign), as if for absolute value
+            } else { // if (gt_half_ulp)
+              is_inexact_lt_midpoint = 1; //if(z_sign), as if for absolute value
+            }
+
+            if (rnd_mode != ROUNDING_TO_NEAREST) {
+              rounding_correction (rnd_mode,
+                  is_inexact_lt_midpoint,
+                  is_inexact_gt_midpoint,
+                  is_midpoint_lt_even,
+                  is_midpoint_gt_even, e3, &res,
+                  pfpsf);
+              z_exp = res.w[1] & MASK_EXP;
+            }
+          } // end e3 = emin
+          // set the inexact flag (if the result was not exact)
+          if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+              is_midpoint_lt_even || is_midpoint_gt_even)
+            *pfpsf |= INEXACT_EXCEPTION;
+        } // end 10^33
+      } // end if (p_sign != z_sign)
+      res.w[1] = z_sign | (z_exp & MASK_EXP) | res.w[1];
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else if (((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2)
+        (q3 <= delta && delta + q4 <= p34) || // Case (3)
+        (delta < q3 && p34 < delta + q4) || // Case (4)
+        (delta < q3 && q3 <= delta + q4 && delta + q4 <= p34) || // Case (5)
+        (delta + q4 < q3)) && // Case (6)
+        !(delta <= 1 && p_sign != z_sign)) { // Case (2), (3), (4), (5) or (6)
+
+      // the result has the sign of z
+
+      if ((q3 <= delta && delta < p34 && p34 < delta + q4) || // Case (2)
+          (delta < q3 && p34 < delta + q4)) { // Case (4)
+        // round first the sum x * y + z with unbounded exponent
+        // scale C3 up by scale = p34 - q3, 1 <= scale <= p34-1, 
+        // 1 <= scale <= 33
+        // calculate res = C3 * 10^scale
+        scale = p34 - q3;
+        x0 = delta + q4 - p34;
+      } else if (delta + q4 < q3) { // Case (6)
+        // make Case (6) look like Case (3) or Case (5) with scale = 0
+        // by scaling up C4 by 10^(q3 - delta - q4) 
+        scale = q3 - delta - q4; // 1 <= scale <= 33
+        if (q4 <= 19) { // 1 <= scale <= 19; C4 fits in 64 bits
+          if (scale <= 19) { // 10^scale fits in 64 bits
+            // 64 x 64 C4.w[0] * ten2k64[scale]
+            __mul_64x64_to_128MACH (P128, C4.w[0], ten2k64[scale]);
+          } else { // 10^scale fits in 128 bits
+            // 64 x 128 C4.w[0] * ten2k128[scale - 20]
+            __mul_128x64_to_128 (P128, C4.w[0], ten2k128[scale - 20]);
+          }
+        } else { // C4 fits in 128 bits, but 10^scale must fit in 64 bits
+          // 64 x 128 ten2k64[scale] * C4
+          __mul_128x64_to_128 (P128, ten2k64[scale], C4);
+        }
+        C4.w[0] = P128.w[0];
+        C4.w[1] = P128.w[1];
+        // e4 does not need adjustment, as it is not used from this point on
+        scale = 0;
+        x0 = 0;
+        // now Case (6) looks like Case (3) or Case (5) with scale = 0 
+      } else { // if Case (3) or Case (5)
+        // Note: Case (3) is similar to Case (2), but scale differs and the
+        // result is exact, unless it is tiny (so x0 = 0 when calculating the
+        // result with unbounded exponent)
+
+        // calculate first the sum x * y + z with unbounded exponent (exact)
+        // scale C3 up by scale = delta + q4 - q3, 1 <= scale <= p34-1,
+        // 1 <= scale <= 33
+        // calculate res = C3 * 10^scale
+        scale = delta + q4 - q3;
+        x0 = 0;
+        // Note: the comments which follow refer [mainly] to Case (2)]
+      }
+
+    case2_repeat:
+      if (scale == 0) { // this could happen e.g. if we return to case2_repeat
+        // or in Case (4)
+        res.w[1] = C3.w[1];
+        res.w[0] = C3.w[0];
+      } else if (q3 <= 19) { // 1 <= scale <= 19; z fits in 64 bits
+        if (scale <= 19) { // 10^scale fits in 64 bits
+          // 64 x 64 C3.w[0] * ten2k64[scale]
+          __mul_64x64_to_128MACH (res, C3.w[0], ten2k64[scale]);
+        } else { // 10^scale fits in 128 bits
+          // 64 x 128 C3.w[0] * ten2k128[scale - 20]
+          __mul_128x64_to_128 (res, C3.w[0], ten2k128[scale - 20]);
+        }
+      } else { // z fits in 128 bits, but 10^scale must fit in 64 bits
+        // 64 x 128 ten2k64[scale] * C3
+        __mul_128x64_to_128 (res, ten2k64[scale], C3);
+      }
+      // e3 is already calculated
+      e3 = e3 - scale;
+      // now res = C3 * 10^scale and e3 = e3 - scale
+      // Note: C3 * 10^scale could be 10^34 if we returned to case2_repeat
+      // because the result was too small
+
+      // round C4 to nearest to q4 - x0 digits, where x0 = delta + q4 - p34,
+      // 1 <= x0 <= min (q4 - 1, 2 * p34 - 1) <=> 1 <= x0 <= min (q4 - 1, 67)
+      // Also: 1 <= q4 - x0 <= p34 -1 => 1 <= q4 - x0 <= 33 (so the result of
+      // the rounding fits in 128 bits!)
+      // x0 = delta + q4 - p34 (calculated before reaching case2_repeat)
+      // because q3 + q4 - x0 <= P => x0 >= q3 + q4 - p34
+      if (x0 == 0) { // this could happen only if we return to case2_repeat, or
+        // for Case (3) or Case (6)
+        R128.w[1] = C4.w[1];
+        R128.w[0] = C4.w[0];
+      } else if (q4 <= 18) {
+        // 2 <= q4 <= 18, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 17
+        round64_2_18 (q4, x0, C4.w[0], &R64, &incr_exp,
+            &is_midpoint_lt_even, &is_midpoint_gt_even,
+            &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+        if (incr_exp) {
+          // R64 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 17
+          R64 = ten2k64[q4 - x0];
+        }
+        R128.w[1] = 0;
+        R128.w[0] = R64;
+      } else if (q4 <= 38) {
+        // 19 <= q4 <= 38, max(1, q3+q4-p34) <= x0 <= q4 - 1, 1 <= x0 <= 37
+        P128.w[1] = C4.w[1];
+        P128.w[0] = C4.w[0];
+        round128_19_38 (q4, x0, P128, &R128, &incr_exp,
+            &is_midpoint_lt_even, &is_midpoint_gt_even,
+            &is_inexact_lt_midpoint,
+            &is_inexact_gt_midpoint);
+        if (incr_exp) {
+          // R128 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 1, 1 <= q4 - x0 <= 37
+          if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
+            R128.w[0] = ten2k64[q4 - x0];
+            // R128.w[1] stays 0
+          } else { // 20 <= q4 - x0 <= 37
+            R128.w[0] = ten2k128[q4 - x0 - 20].w[0];
+            R128.w[1] = ten2k128[q4 - x0 - 20].w[1];
+          }
+        }
+      } else if (q4 <= 57) {
+        // 38 <= q4 <= 57, max(1, q3+q4-p34) <= x0 <= q4 - 1, 5 <= x0 <= 56
+        P192.w[2] = C4.w[2];
+        P192.w[1] = C4.w[1];
+        P192.w[0] = C4.w[0];
+        round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+            &is_midpoint_lt_even, &is_midpoint_gt_even,
+            &is_inexact_lt_midpoint,
+            &is_inexact_gt_midpoint);
+        // R192.w[2] is always 0
+        if (incr_exp) {
+          // R192 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 5, 1 <= q4 - x0 <= 52
+          if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19
+            R192.w[0] = ten2k64[q4 - x0];
+            // R192.w[1] stays 0
+            // R192.w[2] stays 0
+          } else { // 20 <= q4 - x0 <= 33
+            R192.w[0] = ten2k128[q4 - x0 - 20].w[0];
+            R192.w[1] = ten2k128[q4 - x0 - 20].w[1];
+            // R192.w[2] stays 0
+          }
+        }
+        R128.w[1] = R192.w[1];
+        R128.w[0] = R192.w[0];
+      } else {
+        // 58 <= q4 <= 68, max(1, q3+q4-p34) <= x0 <= q4 - 1, 25 <= x0 <= 67
+        round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+            &is_midpoint_lt_even, &is_midpoint_gt_even,
+            &is_inexact_lt_midpoint,
+            &is_inexact_gt_midpoint);
+        // R256.w[3] and R256.w[2] are always 0
+        if (incr_exp) {
+          // R256 = 10^(q4-x0), 1 <= q4 - x0 <= q4 - 25, 1 <= q4 - x0 <= 43
+          if (q4 - x0 <= 19) { // 1 <= q4 - x0 <= 19  
+            R256.w[0] = ten2k64[q4 - x0];
+            // R256.w[1] stays 0
+            // R256.w[2] stays 0
+            // R256.w[3] stays 0
+          } else { // 20 <= q4 - x0 <= 33 
+            R256.w[0] = ten2k128[q4 - x0 - 20].w[0];
+            R256.w[1] = ten2k128[q4 - x0 - 20].w[1];
+            // R256.w[2] stays 0
+            // R256.w[3] stays 0
+          }
+        }
+        R128.w[1] = R256.w[1];
+        R128.w[0] = R256.w[0];
+      }
+      // now add C3 * 10^scale in res and the signed top (q4-x0) digits of C4,
+      // rounded to nearest, which were copied into R128
+      if (z_sign == p_sign) {
+        lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale
+        // the sum can result in [up to] p34 or p34 + 1 digits
+        res.w[0] = res.w[0] + R128.w[0];
+        res.w[1] = res.w[1] + R128.w[1];
+        if (res.w[0] < R128.w[0])
+          res.w[1]++; // carry
+        // if res > 10^34 - 1 need to increase x0 and decrease scale by 1
+        if (res.w[1] > 0x0001ed09bead87c0ull ||
+            (res.w[1] == 0x0001ed09bead87c0ull &&
+             res.w[0] > 0x378d8e63ffffffffull)) {
+          // avoid double rounding error
+          is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+          is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+          is_midpoint_lt_even0 = is_midpoint_lt_even;
+          is_midpoint_gt_even0 = is_midpoint_gt_even;
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 0;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+          P128.w[1] = res.w[1];
+          P128.w[0] = res.w[0];
+          round128_19_38 (35, 1, P128, &res, &incr_exp,
+              &is_midpoint_lt_even, &is_midpoint_gt_even,
+              &is_inexact_lt_midpoint,
+              &is_inexact_gt_midpoint);
+          // incr_exp is 0 with certainty in this case
+          // avoid a double rounding error
+          if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+              is_midpoint_lt_even) { // double rounding error upward
+            // res = res - 1
+            res.w[0]--;
+            if (res.w[0] == 0xffffffffffffffffull)
+              res.w[1]--;
+            // Note: a double rounding error upward is not possible; for this
+            // the result after the first rounding would have to be 99...95
+            // (35 digits in all), possibly followed by a number of zeros; this
+            // not possible in Cases (2)-(6) or (15)-(17) which may get here
+            is_midpoint_lt_even = 0;
+            is_inexact_lt_midpoint = 1;
+          } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+              is_midpoint_gt_even) { // double rounding error downward
+            // res = res + 1
+            res.w[0]++;
+            if (res.w[0] == 0)
+              res.w[1]++;
+            is_midpoint_gt_even = 0;
+            is_inexact_gt_midpoint = 1;
+          } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	     !is_inexact_lt_midpoint
+        	     && !is_inexact_gt_midpoint) {
+            // if this second rounding was exact the result may still be 
+            // inexact because of the first rounding
+            if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+              is_inexact_gt_midpoint = 1;
+            }
+            if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+              is_inexact_lt_midpoint = 1;
+            }
+          } else if (is_midpoint_gt_even &&
+        	     (is_inexact_gt_midpoint0
+        	      || is_midpoint_lt_even0)) {
+            // pulled up to a midpoint
+            is_inexact_lt_midpoint = 1;
+            is_inexact_gt_midpoint = 0;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+          } else if (is_midpoint_lt_even &&
+        	     (is_inexact_lt_midpoint0
+        	      || is_midpoint_gt_even0)) {
+            // pulled down to a midpoint
+            is_inexact_lt_midpoint = 0;
+            is_inexact_gt_midpoint = 1;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+          } else {
+            ;
+          }
+          // adjust exponent
+          e3 = e3 + 1;
+          if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+              !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+            if (is_midpoint_lt_even0 || is_midpoint_gt_even0 ||
+        	is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
+              is_inexact_lt_midpoint = 1;
+            }
+          }
+        } else {
+          // this is the result rounded with unbounded exponent, unless a
+          // correction is needed
+          res.w[1] = res.w[1] & MASK_COEFF;
+          if (lsb == 1) {
+            if (is_midpoint_gt_even) {
+              // res = res + 1
+              is_midpoint_gt_even = 0;
+              is_midpoint_lt_even = 1;
+              res.w[0]++;
+              if (res.w[0] == 0x0)
+        	res.w[1]++;
+              // check for rounding overflow
+              if (res.w[1] == 0x0001ed09bead87c0ull &&
+        	  res.w[0] == 0x378d8e6400000000ull) {
+        	// res = 10^34 => rounding overflow
+        	res.w[1] = 0x0000314dc6448d93ull;
+        	res.w[0] = 0x38c15b0a00000000ull; // 10^33
+        	e3++;
+              }
+            } else if (is_midpoint_lt_even) {
+              // res = res - 1
+              is_midpoint_lt_even = 0;
+              is_midpoint_gt_even = 1;
+              res.w[0]--;
+              if (res.w[0] == 0xffffffffffffffffull)
+        	res.w[1]--;
+              // if the result is pure zero, the sign depends on the rounding 
+              // mode (x*y and z had opposite signs)
+              if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
+        	if (rnd_mode != ROUNDING_DOWN)
+        	  z_sign = 0x0000000000000000ull;
+        	else
+        	  z_sign = 0x8000000000000000ull;
+        	// the exponent is max (e3, expmin)
+        	res.w[1] = 0x0;
+        	res.w[0] = 0x0;
+        	*ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+        	*ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+        	*ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+        	*ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+        	BID_SWAP128 (res);
+        	BID_RETURN (res)
+              }
+            } else {
+              ;
+            }
+          }
+        }
+      } else { // if (z_sign != p_sign)
+        lsb = res.w[0] & 0x01; // lsb of C3 * 10^scale; R128 contains rounded C4
+        // used to swap rounding indicators if p_sign != z_sign
+        // the sum can result in [up to] p34 or p34 - 1 digits
+        tmp64 = res.w[0];
+        res.w[0] = res.w[0] - R128.w[0];
+        res.w[1] = res.w[1] - R128.w[1];
+        if (res.w[0] > tmp64)
+          res.w[1]--; // borrow
+        // if res < 10^33 and exp > expmin need to decrease x0 and 
+        // increase scale by 1
+        if (e3 > expmin && ((res.w[1] < 0x0000314dc6448d93ull ||
+        		     (res.w[1] == 0x0000314dc6448d93ull &&
+        		      res.w[0] < 0x38c15b0a00000000ull)) ||
+        		    (is_inexact_lt_midpoint
+        		     && res.w[1] == 0x0000314dc6448d93ull
+        		     && res.w[0] == 0x38c15b0a00000000ull))
+            && x0 >= 1) {
+          x0 = x0 - 1;
+          // first restore e3, otherwise it will be too small
+          e3 = e3 + scale;
+          scale = scale + 1;
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 0;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+          incr_exp = 0;
+          goto case2_repeat;
+        }
+        // else this is the result rounded with unbounded exponent;
+        // because the result has opposite sign to that of C4 which was 
+        // rounded, need to change the rounding indicators
+        if (is_inexact_lt_midpoint) {
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 1;
+        } else if (is_inexact_gt_midpoint) {
+          is_inexact_gt_midpoint = 0;
+          is_inexact_lt_midpoint = 1;
+        } else if (lsb == 0) {
+          if (is_midpoint_lt_even) {
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 1;
+          } else if (is_midpoint_gt_even) {
+            is_midpoint_gt_even = 0;
+            is_midpoint_lt_even = 1;
+          } else {
+            ;
+          }
+        } else if (lsb == 1) {
+          if (is_midpoint_lt_even) {
+            // res = res + 1
+            res.w[0]++;
+            if (res.w[0] == 0x0)
+              res.w[1]++;
+            // check for rounding overflow
+            if (res.w[1] == 0x0001ed09bead87c0ull &&
+        	res.w[0] == 0x378d8e6400000000ull) {
+              // res = 10^34 => rounding overflow
+              res.w[1] = 0x0000314dc6448d93ull;
+              res.w[0] = 0x38c15b0a00000000ull; // 10^33
+              e3++;
+            }
+          } else if (is_midpoint_gt_even) {
+            // res = res - 1
+            res.w[0]--;
+            if (res.w[0] == 0xffffffffffffffffull)
+              res.w[1]--;
+            // if the result is pure zero, the sign depends on the rounding 
+            // mode (x*y and z had opposite signs)
+            if (res.w[1] == 0x0ull && res.w[0] == 0x0ull) {
+              if (rnd_mode != ROUNDING_DOWN)
+        	z_sign = 0x0000000000000000ull;
+              else
+        	z_sign = 0x8000000000000000ull;
+              // the exponent is max (e3, expmin)
+              res.w[1] = 0x0;
+              res.w[0] = 0x0;
+              *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+              *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+              *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+              *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+              BID_SWAP128 (res);
+              BID_RETURN (res)
+            }
+          } else {
+            ;
+          }
+        } else {
+          ;
+        }
+      }
+      // check for underflow
+      if (e3 == expmin) { // and if significand < 10^33 => result is tiny
+        if ((res.w[1] & MASK_COEFF) < 0x0000314dc6448d93ull ||
+            ((res.w[1] & MASK_COEFF) == 0x0000314dc6448d93ull &&
+             res.w[0] < 0x38c15b0a00000000ull)) {
+          is_tiny = 1;
+        }
+      } else if (e3 < expmin) {
+        // the result is tiny, so we must truncate more of res
+        is_tiny = 1;
+        x0 = expmin - e3;
+        is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+        is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+        is_midpoint_lt_even0 = is_midpoint_lt_even;
+        is_midpoint_gt_even0 = is_midpoint_gt_even;
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+        // determine the number of decimal digits in res
+        if (res.w[1] == 0x0) {
+          // between 1 and 19 digits
+          for (ind = 1; ind <= 19; ind++) {
+            if (res.w[0] < ten2k64[ind]) {
+              break;
+            }
+          }
+          // ind digits
+        } else if (res.w[1] < ten2k128[0].w[1] ||
+        	   (res.w[1] == ten2k128[0].w[1]
+        	    && res.w[0] < ten2k128[0].w[0])) {
+          // 20 digits
+          ind = 20;
+        } else { // between 21 and 38 digits
+          for (ind = 1; ind <= 18; ind++) {
+            if (res.w[1] < ten2k128[ind].w[1] ||
+        	(res.w[1] == ten2k128[ind].w[1] &&
+        	 res.w[0] < ten2k128[ind].w[0])) {
+              break;
+            }
+          }
+          // ind + 20 digits
+          ind = ind + 20;
+        }
+
+        // at this point ind >= x0; because delta >= 2 on this path, the case
+        // ind = x0 can occur only in Case (2) or case (3), when C3 has one
+        // digit (q3 = 1) equal to 1 (C3 = 1), e3 is expmin (e3 = expmin), 
+        // the signs of x * y and z are opposite, and through cancellation 
+        // the most significant decimal digit in res has the weight
+        // 10^(emin-1); however, it is clear that in this case the most
+        // significant digit is 9, so the result before rounding is
+        // 0.9... * 10^emin
+        // Otherwise, ind > x0 because there are non-zero decimal digits in the
+        // result with weight of at least 10^emin, and correction for underflow
+        //  can be carried out using the round*_*_2_* () routines
+        if (x0 == ind) { // the result before rounding is 0.9... * 10^emin
+          res.w[1] = 0x0;
+          res.w[0] = 0x1;
+          is_inexact_gt_midpoint = 1;
+        } else if (ind <= 18) { // check that 2 <= ind
+          // 2 <= ind <= 18, 1 <= x0 <= 17
+          round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+          if (incr_exp) {
+            // R64 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 17
+            R64 = ten2k64[ind - x0];
+          }
+          res.w[1] = 0;
+          res.w[0] = R64;
+        } else if (ind <= 38) {
+          // 19 <= ind <= 38
+          P128.w[1] = res.w[1];
+          P128.w[0] = res.w[0];
+          round128_19_38 (ind, x0, P128, &res, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+          if (incr_exp) {
+            // R128 = 10^(ind-x0), 1 <= ind - x0 <= ind - 1, 1 <= ind - x0 <= 37
+            if (ind - x0 <= 19) { // 1 <= ind - x0 <= 19
+              res.w[0] = ten2k64[ind - x0];
+              // res.w[1] stays 0
+            } else { // 20 <= ind - x0 <= 37
+              res.w[0] = ten2k128[ind - x0 - 20].w[0];
+              res.w[1] = ten2k128[ind - x0 - 20].w[1];
+            }
+          }
+        }
+        // avoid a double rounding error
+        if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+            is_midpoint_lt_even) { // double rounding error upward
+          // res = res - 1
+          res.w[0]--;
+          if (res.w[0] == 0xffffffffffffffffull)
+            res.w[1]--;
+          // Note: a double rounding error upward is not possible; for this
+          // the result after the first rounding would have to be 99...95
+          // (35 digits in all), possibly followed by a number of zeros; this
+          // not possible in Cases (2)-(6) which may get here
+          is_midpoint_lt_even = 0;
+          is_inexact_lt_midpoint = 1;
+        } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+            is_midpoint_gt_even) { // double rounding error downward
+          // res = res + 1
+          res.w[0]++;
+          if (res.w[0] == 0)
+            res.w[1]++;
+          is_midpoint_gt_even = 0;
+          is_inexact_gt_midpoint = 1;
+        } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+          // if this second rounding was exact the result may still be 
+          // inexact because of the first rounding
+          if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+            is_inexact_gt_midpoint = 1;
+          }
+          if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+            is_inexact_lt_midpoint = 1;
+          }
+        } else if (is_midpoint_gt_even &&
+        	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+          // pulled up to a midpoint
+          is_inexact_lt_midpoint = 1;
+          is_inexact_gt_midpoint = 0;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else if (is_midpoint_lt_even &&
+        	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+          // pulled down to a midpoint
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 1;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else {
+          ;
+        }
+        // adjust exponent
+        e3 = e3 + x0;
+        if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+            !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+          if (is_midpoint_lt_even0 || is_midpoint_gt_even0 ||
+              is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0) {
+            is_inexact_lt_midpoint = 1;
+          }
+        }
+      } else {
+        ; // not underflow
+      }
+      // check for inexact result
+      if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+          is_midpoint_lt_even || is_midpoint_gt_even) {
+        // set the inexact flag
+        *pfpsf |= INEXACT_EXCEPTION;
+        if (is_tiny)
+          *pfpsf |= UNDERFLOW_EXCEPTION;
+      }
+      // now check for significand = 10^34 (may have resulted from going
+      // back to case2_repeat)
+      if (res.w[1] == 0x0001ed09bead87c0ull && 
+          res.w[0] == 0x378d8e6400000000ull) { // if  res = 10^34
+        res.w[1] = 0x0000314dc6448d93ull; // res = 10^33
+        res.w[0] = 0x38c15b0a00000000ull;
+        e3 = e3 + 1;
+      }
+      res.w[1] = z_sign | ((UINT64) (e3 + 6176) << 49) | res.w[1];
+      // check for overflow
+      if (rnd_mode == ROUNDING_TO_NEAREST && e3 > expmax) {
+        res.w[1] = z_sign | 0x7800000000000000ull; // +/-inf
+        res.w[0] = 0x0000000000000000ull;
+        *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+      }
+      if (rnd_mode != ROUNDING_TO_NEAREST) {
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e3, &res, pfpsf);
+      }
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else {
+
+      // we get here only if delta <= 1 in Cases (2), (3), (4), (5), or (6) and
+      // the signs of x*y and z are opposite; in these cases massive
+      // cancellation can occur, so it is better to scale either C3 or C4 and 
+      // to perform the subtraction before rounding; rounding is performed 
+      // next, depending on the number of decimal digits in the result and on 
+      // the exponent value
+      // Note: overlow is not possible in this case
+      // this is similar to Cases (15), (16), and (17)
+
+      if (delta + q4 < q3) { // from Case (6) 
+        // Case (6) with 0<= delta <= 1 is similar to Cases (15), (16), and 
+        // (17) if we swap (C3, C4), (q3, q4), (e3, e4), (z_sign, p_sign)
+        // and call add_and_round; delta stays positive
+        // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3
+        P128.w[1] = C3.w[1];
+        P128.w[0] = C3.w[0];
+        C3.w[1] = C4.w[1];
+        C3.w[0] = C4.w[0];
+        C4.w[1] = P128.w[1];
+        C4.w[0] = P128.w[0];
+        ind = q3;
+        q3 = q4;
+        q4 = ind;
+        ind = e3;
+        e3 = e4;
+        e4 = ind;
+        tmp_sign = z_sign;
+        z_sign = p_sign;
+        p_sign = tmp_sign;
+      } else { // from Cases (2), (3), (4), (5)
+        // In Cases (2), (3), (4), (5) with 0 <= delta <= 1 C3 has to be 
+        // scaled up by q4 + delta - q3; this is the same as in Cases (15), 
+        // (16), and (17) if we just change the sign of delta
+        delta = -delta;
+      }
+      add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
+        	     rnd_mode, &is_midpoint_lt_even,
+        	     &is_midpoint_gt_even, &is_inexact_lt_midpoint,
+        	     &is_inexact_gt_midpoint, pfpsf, &res);
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    }
+
+  } else { // if delta < 0
+
+    delta = -delta;
+
+    if (p34 < q4 && q4 <= delta) { // Case (7)
+
+      // truncate C4 to p34 digits into res
+      // x = q4-p34, 1 <= x <= 34 because 35 <= q4 <= 68
+      x0 = q4 - p34;
+      if (q4 <= 38) {
+        P128.w[1] = C4.w[1];
+        P128.w[0] = C4.w[0];
+        round128_19_38 (q4, x0, P128, &res, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+      } else if (q4 <= 57) { // 35 <= q4 <= 57
+        P192.w[2] = C4.w[2];
+        P192.w[1] = C4.w[1];
+        P192.w[0] = C4.w[0];
+        round192_39_57 (q4, x0, P192, &R192, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+        res.w[0] = R192.w[0];
+        res.w[1] = R192.w[1];
+      } else { // if (q4 <= 68)
+        round256_58_76 (q4, x0, C4, &R256, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+        res.w[0] = R256.w[0];
+        res.w[1] = R256.w[1];
+      }
+      e4 = e4 + x0;
+      if (incr_exp) {
+        e4 = e4 + 1;
+      }
+      if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+          !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+        // if C4 rounded to p34 digits is exact then the result is inexact,
+        // in a way that depends on the signs of x * y and z
+        if (p_sign == z_sign) {
+          is_inexact_lt_midpoint = 1;
+        } else { // if (p_sign != z_sign)
+          if (res.w[1] != 0x0000314dc6448d93ull || 
+              res.w[0] != 0x38c15b0a00000000ull) { // res != 10^33
+            is_inexact_gt_midpoint = 1;
+          } else { // res = 10^33 and exact is a special case
+            // if C3 < 1/2 ulp then res = 10^33 and is_inexact_gt_midpoint = 1
+            // if C3 = 1/2 ulp then res = 10^33 and is_midpoint_lt_even = 1
+            // if C3 > 1/2 ulp then res = 10^34-1 and is_inexact_lt_midpoint = 1
+            // Note: ulp is really ulp/10 (after borrow which propagates to msd)
+            if (delta > p34 + 1) { // C3 < 1/2
+              // res = 10^33, unchanged
+              is_inexact_gt_midpoint = 1;
+            } else { // if (delta == p34 + 1)
+              if (q3 <= 19) {
+        	if (C3.w[0] < midpoint64[q3 - 1]) { // C3 < 1/2 ulp
+        	  // res = 10^33, unchanged
+        	  is_inexact_gt_midpoint = 1;
+        	} else if (C3.w[0] == midpoint64[q3 - 1]) { // C3 = 1/2 ulp
+        	  // res = 10^33, unchanged
+        	  is_midpoint_lt_even = 1;
+        	} else { // if (C3.w[0] > midpoint64[q3-1]), C3 > 1/2 ulp
+        	  res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+        	  res.w[0] = 0x378d8e63ffffffffull;
+        	  e4 = e4 - 1;
+        	  is_inexact_lt_midpoint = 1;
+        	}
+              } else { // if (20 <= q3 <=34)
+        	if (C3.w[1] < midpoint128[q3 - 20].w[1] || 
+                    (C3.w[1] == midpoint128[q3 - 20].w[1] && 
+                    C3.w[0] < midpoint128[q3 - 20].w[0])) { // C3 < 1/2 ulp
+        	  // res = 10^33, unchanged
+        	  is_inexact_gt_midpoint = 1;
+        	} else if (C3.w[1] == midpoint128[q3 - 20].w[1] && 
+                    C3.w[0] == midpoint128[q3 - 20].w[0]) { // C3 = 1/2 ulp
+        	  // res = 10^33, unchanged
+        	  is_midpoint_lt_even = 1;
+        	} else { // if (C3 > midpoint128[q3-20]), C3 > 1/2 ulp
+        	  res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+        	  res.w[0] = 0x378d8e63ffffffffull;
+        	  e4 = e4 - 1;
+        	  is_inexact_lt_midpoint = 1;
+        	}
+              }
+            }
+          }
+        }
+      } else if (is_midpoint_lt_even) {
+        if (z_sign != p_sign) {
+          // needs correction: res = res - 1
+          res.w[0] = res.w[0] - 1;
+          if (res.w[0] == 0xffffffffffffffffull)
+            res.w[1]--;
+          // if it is (10^33-1)*10^e4 then the corect result is 
+          // (10^34-1)*10(e4-1)
+          if (res.w[1] == 0x0000314dc6448d93ull &&
+              res.w[0] == 0x38c15b09ffffffffull) {
+            res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+            res.w[0] = 0x378d8e63ffffffffull;
+            e4 = e4 - 1;
+          }
+          is_midpoint_lt_even = 0;
+          is_inexact_lt_midpoint = 1;
+        } else { // if (z_sign == p_sign)
+          is_midpoint_lt_even = 0;
+          is_inexact_gt_midpoint = 1;
+        }
+      } else if (is_midpoint_gt_even) {
+        if (z_sign == p_sign) {
+          // needs correction: res = res + 1 (cannot cross in the next binade)
+          res.w[0] = res.w[0] + 1;
+          if (res.w[0] == 0x0000000000000000ull)
+            res.w[1]++;
+          is_midpoint_gt_even = 0;
+          is_inexact_gt_midpoint = 1;
+        } else { // if (z_sign != p_sign)
+          is_midpoint_gt_even = 0;
+          is_inexact_lt_midpoint = 1;
+        }
+      } else {
+        ; // the rounded result is already correct
+      }
+      // check for overflow
+      if (rnd_mode == ROUNDING_TO_NEAREST && e4 > expmax) {
+        res.w[1] = p_sign | 0x7800000000000000ull;
+        res.w[0] = 0x0000000000000000ull;
+        *pfpsf |= (OVERFLOW_EXCEPTION | INEXACT_EXCEPTION);
+      } else { // no overflow or not RN
+        p_exp = ((UINT64) (e4 + 6176) << 49);
+        res.w[1] = p_sign | (p_exp & MASK_EXP) | res.w[1];
+      }
+      if (rnd_mode != ROUNDING_TO_NEAREST) {
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e4, &res, pfpsf);
+      }
+      if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+          is_midpoint_lt_even || is_midpoint_gt_even) {
+        // set the inexact flag
+        *pfpsf |= INEXACT_EXCEPTION;
+      }
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else if ((q4 <= p34 && p34 <= delta) || // Case (8)
+        (q4 <= delta && delta < p34 && p34 < delta + q3) || // Case (9)
+        (q4 <= delta && delta + q3 <= p34) || // Case (10)
+        (delta < q4 && q4 <= p34 && p34 < delta + q3) || // Case (13)
+        (delta < q4 && q4 <= delta + q3 && delta + q3 <= p34) || // Case (14)
+        (delta + q3 < q4 && q4 <= p34)) { // Case (18)
+
+      // Case (8) is similar to Case (1), with C3 and C4 swapped
+      // Case (9) is similar to Case (2), with C3 and C4 swapped
+      // Case (10) is similar to Case (3), with C3 and C4 swapped
+      // Case (13) is similar to Case (4), with C3 and C4 swapped
+      // Case (14) is similar to Case (5), with C3 and C4 swapped
+      // Case (18) is similar to Case (6), with C3 and C4 swapped
+
+      // swap (C3, C4), (q3, q4), (e3, 34), (z_sign, p_sign), (z_exp, p_exp)
+      // and go back to delta_ge_zero
+      // C4.w[3] = 0 and C4.w[2] = 0, so swap just the low part of C4 with C3
+      P128.w[1] = C3.w[1];
+      P128.w[0] = C3.w[0];
+      C3.w[1] = C4.w[1];
+      C3.w[0] = C4.w[0];
+      C4.w[1] = P128.w[1];
+      C4.w[0] = P128.w[0];
+      ind = q3;
+      q3 = q4;
+      q4 = ind;
+      ind = e3;
+      e3 = e4;
+      e4 = ind;
+      tmp_sign = z_sign;
+      z_sign = p_sign;
+      p_sign = tmp_sign;
+      tmp.ui64 = z_exp;
+      z_exp = p_exp;
+      p_exp = tmp.ui64;
+      goto delta_ge_zero;
+
+    } else if ((p34 <= delta && delta < q4 && q4 < delta + q3) || // Case (11)
+               (delta < p34 && p34 < q4 && q4 < delta + q3)) { // Case (12)
+
+      // round C3 to nearest to q3 - x0 digits, where x0 = e4 - e3,
+      // 1 <= x0 <= q3 - 1 <= p34 - 1 
+      x0 = e4 - e3; // or x0 = delta + q3 - q4
+      if (q3 <= 18) { // 2 <= q3 <= 18
+        round64_2_18 (q3, x0, C3.w[0], &R64, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+        // C3.w[1] = 0;
+        C3.w[0] = R64;
+      } else if (q3 <= 38) {
+        round128_19_38 (q3, x0, C3, &R128, &incr_exp,
+        		&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint);
+        C3.w[1] = R128.w[1];
+        C3.w[0] = R128.w[0];
+      }
+      // the rounded result has q3 - x0 digits
+      // we want the exponent to be e4, so if incr_exp = 1 then
+      // multiply the rounded result by 10 - it will still fit in 113 bits
+      if (incr_exp) {
+        // 64 x 128 -> 128
+        P128.w[1] = C3.w[1];
+        P128.w[0] = C3.w[0];
+        __mul_64x128_to_128 (C3, ten2k64[1], P128);
+      }
+      e3 = e3 + x0; // this is e4
+      // now add/subtract the 256-bit C4 and the new (and shorter) 128-bit C3; 
+      // the result will have the sign of x * y; the exponent is e4
+      R256.w[3] = 0;
+      R256.w[2] = 0;
+      R256.w[1] = C3.w[1];
+      R256.w[0] = C3.w[0];
+      if (p_sign == z_sign) { // R256 = C4 + R256
+        add256 (C4, R256, &R256);
+      } else { // if (p_sign != z_sign) { // R256 = C4 - R256
+        sub256 (C4, R256, &R256); // the result cannot be pure zero
+        // because the result has opposite sign to that of R256 which was 
+        // rounded, need to change the rounding indicators
+        lsb = C4.w[0] & 0x01;
+        if (is_inexact_lt_midpoint) {
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 1;
+        } else if (is_inexact_gt_midpoint) {
+          is_inexact_gt_midpoint = 0;
+          is_inexact_lt_midpoint = 1;
+        } else if (lsb == 0) {
+          if (is_midpoint_lt_even) {
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 1;
+          } else if (is_midpoint_gt_even) {
+            is_midpoint_gt_even = 0;
+            is_midpoint_lt_even = 1;
+          } else {
+            ;
+          }
+        } else if (lsb == 1) {
+          if (is_midpoint_lt_even) {
+            // res = res + 1
+            R256.w[0]++;
+            if (R256.w[0] == 0x0) {
+              R256.w[1]++;
+              if (R256.w[1] == 0x0) {
+        	R256.w[2]++;
+        	if (R256.w[2] == 0x0) {
+        	  R256.w[3]++;
+        	}
+              }
+            }
+            // no check for rounding overflow - R256 was a difference
+          } else if (is_midpoint_gt_even) {
+            // res = res - 1
+            R256.w[0]--;
+            if (R256.w[0] == 0xffffffffffffffffull) {
+              R256.w[1]--;
+              if (R256.w[1] == 0xffffffffffffffffull) {
+        	R256.w[2]--;
+        	if (R256.w[2] == 0xffffffffffffffffull) {
+        	  R256.w[3]--;
+        	}
+              }
+            }
+          } else {
+            ;
+          }
+        } else {
+          ;
+        }
+      }
+      // determine the number of decimal digits in R256
+      ind = nr_digits256 (R256); // ind >= p34
+      // if R256 is sum, then ind > p34; if R256 is a difference, then 
+      // ind >= p34; this means that we can calculate the result rounded to
+      // the destination precision, with unbounded exponent, starting from R256
+      // and using the indicators from the rounding of C3 to avoid a double
+      // rounding error 
+
+      if (ind < p34) {
+        ;
+      } else if (ind == p34) {
+        // the result rounded to the destination precision with 
+        // unbounded exponent
+        // is (-1)^p_sign * R256 * 10^e4
+        res.w[1] = R256.w[1];
+        res.w[0] = R256.w[0];
+      } else { // if (ind > p34)
+        // if more than P digits, round to nearest to P digits
+        // round R256 to p34 digits
+        x0 = ind - p34; // 1 <= x0 <= 34 as 35 <= ind <= 68
+        // save C3 rounding indicators to help avoid double rounding error
+        is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+        is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+        is_midpoint_lt_even0 = is_midpoint_lt_even;
+        is_midpoint_gt_even0 = is_midpoint_gt_even;
+        // initialize rounding indicators
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+        // round to p34 digits; the result fits in 113 bits
+        if (ind <= 38) {
+          P128.w[1] = R256.w[1];
+          P128.w[0] = R256.w[0];
+          round128_19_38 (ind, x0, P128, &R128, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+        } else if (ind <= 57) {
+          P192.w[2] = R256.w[2];
+          P192.w[1] = R256.w[1];
+          P192.w[0] = R256.w[0];
+          round192_39_57 (ind, x0, P192, &R192, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+          R128.w[1] = R192.w[1];
+          R128.w[0] = R192.w[0];
+        } else { // if (ind <= 68)
+          round256_58_76 (ind, x0, R256, &R256, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+          R128.w[1] = R256.w[1];
+          R128.w[0] = R256.w[0];
+        }
+        // the rounded result has p34 = 34 digits
+        e4 = e4 + x0 + incr_exp;
+
+        res.w[1] = R128.w[1];
+        res.w[0] = R128.w[0];
+
+        // avoid a double rounding error
+        if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+            is_midpoint_lt_even) { // double rounding error upward
+          // res = res - 1
+          res.w[0]--;
+          if (res.w[0] == 0xffffffffffffffffull)
+            res.w[1]--;
+          is_midpoint_lt_even = 0;
+          is_inexact_lt_midpoint = 1;
+          // Note: a double rounding error upward is not possible; for this
+          // the result after the first rounding would have to be 99...95
+          // (35 digits in all), possibly followed by a number of zeros; this
+          // not possible in Cases (2)-(6) or (15)-(17) which may get here
+          // if this is 10^33 - 1 make it 10^34 - 1 and decrement exponent
+          if (res.w[1] == 0x0000314dc6448d93ull && 
+            res.w[0] == 0x38c15b09ffffffffull) { // 10^33 - 1
+            res.w[1] = 0x0001ed09bead87c0ull; // 10^34 - 1
+            res.w[0] = 0x378d8e63ffffffffull;
+            e4--;
+          }
+        } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+            is_midpoint_gt_even) { // double rounding error downward
+          // res = res + 1 
+          res.w[0]++;
+          if (res.w[0] == 0)
+            res.w[1]++;
+          is_midpoint_gt_even = 0;
+          is_inexact_gt_midpoint = 1;
+        } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	   !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+          // if this second rounding was exact the result may still be
+          // inexact because of the first rounding
+          if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+            is_inexact_gt_midpoint = 1;
+          }
+          if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+            is_inexact_lt_midpoint = 1;
+          }
+        } else if (is_midpoint_gt_even &&
+        	   (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+          // pulled up to a midpoint
+          is_inexact_lt_midpoint = 1;
+          is_inexact_gt_midpoint = 0;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else if (is_midpoint_lt_even &&
+        	   (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+          // pulled down to a midpoint
+          is_inexact_lt_midpoint = 0;
+          is_inexact_gt_midpoint = 1;
+          is_midpoint_lt_even = 0;
+          is_midpoint_gt_even = 0;
+        } else {
+          ;
+        }
+      }
+
+      // determine tininess
+      if (rnd_mode == ROUNDING_TO_NEAREST) {
+        if (e4 < expmin) {
+          is_tiny = 1; // for other rounding modes apply correction
+        }
+      } else {
+        // for RM, RP, RZ, RA apply correction in order to determine tininess
+        // but do not save the result; apply the correction to 
+        // (-1)^p_sign * res * 10^0
+        P128.w[1] = p_sign | 0x3040000000000000ull | res.w[1];
+        P128.w[0] = res.w[0];
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     0, &P128, pfpsf);
+        scale = ((P128.w[1] & MASK_EXP) >> 49) - 6176; // -1, 0, or +1
+        // the number of digits in the significand is p34 = 34
+        if (e4 + scale < expmin) {
+          is_tiny = 1;
+        }
+      }
+
+      // the result rounded to the destination precision with unbounded exponent
+      // is (-1)^p_sign * res * 10^e4
+      res.w[1] = p_sign | ((UINT64) (e4 + 6176) << 49) | res.w[1]; // RN
+      // res.w[0] unchanged;
+      // Note: res is correct only if expmin <= e4 <= expmax
+      ind = p34; // the number of decimal digits in the signifcand of res
+
+      // at this point we have the result rounded with unbounded exponent in
+      // res and we know its tininess:
+      // res = (-1)^p_sign * significand * 10^e4, 
+      // where q (significand) = ind = p34
+      // Note: res is correct only if expmin <= e4 <= expmax
+
+      // check for overflow if RN
+      if (rnd_mode == ROUNDING_TO_NEAREST
+          && (ind + e4) > (p34 + expmax)) {
+        res.w[1] = p_sign | 0x7800000000000000ull;
+        res.w[0] = 0x0000000000000000ull;
+        *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+        *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+        *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+        *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+        *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+        BID_SWAP128 (res);
+        BID_RETURN (res)
+      } // else not overflow or not RN, so continue
+
+      // from this point on this is similar to the last part of the computation
+      // for Cases (15), (16), (17)
+
+      // if (e4 >= expmin) we have the result rounded with bounded exponent
+      if (e4 < expmin) {
+        x0 = expmin - e4; // x0 >= 1; the number of digits to chop off of res
+        // where the result rounded [at most] once is
+        //   (-1)^p_sign * significand_res * 10^e4
+
+        // avoid double rounding error
+        is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+        is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+        is_midpoint_lt_even0 = is_midpoint_lt_even;
+        is_midpoint_gt_even0 = is_midpoint_gt_even;
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+
+        if (x0 > ind) {
+          // nothing is left of res when moving the decimal point left x0 digits
+          is_inexact_lt_midpoint = 1;
+          res.w[1] = p_sign | 0x0000000000000000ull;
+          res.w[0] = 0x0000000000000000ull;
+          e4 = expmin;
+        } else if (x0 == ind) { // 1 <= x0 = ind <= p34 = 34
+          // this is <, =, or > 1/2 ulp
+          // compare the ind-digit value in the significand of res with
+          // 1/2 ulp = 5*10^(ind-1), i.e. determine whether it is 
+          // less than, equal to, or greater than 1/2 ulp (significand of res)
+          R128.w[1] = res.w[1] & MASK_COEFF;
+          R128.w[0] = res.w[0];
+          if (ind <= 19) {
+            if (R128.w[0] < midpoint64[ind - 1]) { // < 1/2 ulp
+              lt_half_ulp = 1;
+              is_inexact_lt_midpoint = 1;
+            } else if (R128.w[0] == midpoint64[ind - 1]) { // = 1/2 ulp
+              eq_half_ulp = 1;
+              is_midpoint_gt_even = 1;
+            } else { // > 1/2 ulp
+              gt_half_ulp = 1;
+              is_inexact_gt_midpoint = 1;
+            }
+          } else { // if (ind <= 38)
+            if (R128.w[1] < midpoint128[ind - 20].w[1] || 
+                (R128.w[1] == midpoint128[ind - 20].w[1] && 
+                R128.w[0] < midpoint128[ind - 20].w[0])) { // < 1/2 ulp
+              lt_half_ulp = 1;
+              is_inexact_lt_midpoint = 1;
+            } else if (R128.w[1] == midpoint128[ind - 20].w[1] && 
+                R128.w[0] == midpoint128[ind - 20].w[0]) { // = 1/2 ulp
+              eq_half_ulp = 1;
+              is_midpoint_gt_even = 1;
+            } else { // > 1/2 ulp
+              gt_half_ulp = 1;
+              is_inexact_gt_midpoint = 1;
+            }
+          }
+          if (lt_half_ulp || eq_half_ulp) {
+            // res = +0.0 * 10^expmin
+            res.w[1] = 0x0000000000000000ull;
+            res.w[0] = 0x0000000000000000ull;
+          } else { // if (gt_half_ulp)
+            // res = +1 * 10^expmin
+            res.w[1] = 0x0000000000000000ull;
+            res.w[0] = 0x0000000000000001ull;
+          }
+          res.w[1] = p_sign | res.w[1];
+          e4 = expmin;
+        } else { // if (1 <= x0 <= ind - 1 <= 33)
+          // round the ind-digit result to ind - x0 digits
+
+          if (ind <= 18) { // 2 <= ind <= 18
+            round64_2_18 (ind, x0, res.w[0], &R64, &incr_exp,
+        		  &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		  &is_inexact_lt_midpoint,
+        		  &is_inexact_gt_midpoint);
+            res.w[1] = 0x0;
+            res.w[0] = R64;
+          } else if (ind <= 38) {
+            P128.w[1] = res.w[1] & MASK_COEFF;
+            P128.w[0] = res.w[0];
+            round128_19_38 (ind, x0, P128, &res, &incr_exp,
+        		    &is_midpoint_lt_even, &is_midpoint_gt_even,
+        		    &is_inexact_lt_midpoint,
+        		    &is_inexact_gt_midpoint);
+          }
+          e4 = e4 + x0; // expmin
+          // we want the exponent to be expmin, so if incr_exp = 1 then
+          // multiply the rounded result by 10 - it will still fit in 113 bits
+          if (incr_exp) {
+            // 64 x 128 -> 128
+            P128.w[1] = res.w[1] & MASK_COEFF;
+            P128.w[0] = res.w[0];
+            __mul_64x128_to_128 (res, ten2k64[1], P128);
+          }
+          res.w[1] =
+            p_sign | ((UINT64) (e4 + 6176) << 49) | (res.
+        					     w[1] & MASK_COEFF);
+          // avoid a double rounding error
+          if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+                is_midpoint_lt_even) { // double rounding error upward
+            // res = res - 1
+            res.w[0]--;
+            if (res.w[0] == 0xffffffffffffffffull)
+              res.w[1]--;
+            // Note: a double rounding error upward is not possible; for this
+            // the result after the first rounding would have to be 99...95
+            // (35 digits in all), possibly followed by a number of zeros; this
+            // not possible in this underflow case
+            is_midpoint_lt_even = 0;
+            is_inexact_lt_midpoint = 1;
+          } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+                is_midpoint_gt_even) { // double rounding error downward
+            // res = res + 1
+            res.w[0]++;
+            if (res.w[0] == 0)
+              res.w[1]++;
+            is_midpoint_gt_even = 0;
+            is_inexact_gt_midpoint = 1;
+          } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	     !is_inexact_lt_midpoint
+        	     && !is_inexact_gt_midpoint) {
+            // if this second rounding was exact the result may still be 
+            // inexact because of the first rounding
+            if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+              is_inexact_gt_midpoint = 1;
+            }
+            if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+              is_inexact_lt_midpoint = 1;
+            }
+          } else if (is_midpoint_gt_even &&
+        	     (is_inexact_gt_midpoint0
+        	      || is_midpoint_lt_even0)) {
+            // pulled up to a midpoint
+            is_inexact_lt_midpoint = 1;
+            is_inexact_gt_midpoint = 0;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+          } else if (is_midpoint_lt_even &&
+        	     (is_inexact_lt_midpoint0
+        	      || is_midpoint_gt_even0)) {
+            // pulled down to a midpoint
+            is_inexact_lt_midpoint = 0;
+            is_inexact_gt_midpoint = 1;
+            is_midpoint_lt_even = 0;
+            is_midpoint_gt_even = 0;
+          } else {
+            ;
+          }
+        }
+      }
+      // res contains the correct result
+      // apply correction if not rounding to nearest
+      if (rnd_mode != ROUNDING_TO_NEAREST) {
+        rounding_correction (rnd_mode,
+        		     is_inexact_lt_midpoint,
+        		     is_inexact_gt_midpoint,
+        		     is_midpoint_lt_even, is_midpoint_gt_even,
+        		     e4, &res, pfpsf);
+      }
+      if (is_midpoint_lt_even || is_midpoint_gt_even ||
+          is_inexact_lt_midpoint || is_inexact_gt_midpoint) {
+        // set the inexact flag
+        *pfpsf |= INEXACT_EXCEPTION;
+        if (is_tiny)
+          *pfpsf |= UNDERFLOW_EXCEPTION;
+      }
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else if ((p34 <= delta && delta + q3 <= q4) || // Case (15)
+        (delta < p34 && p34 < delta + q3 && delta + q3 <= q4) || //Case (16)
+        (delta + q3 <= p34 && p34 < q4)) { // Case (17)
+
+      // calculate first the result rounded to the destination precision, with
+      // unbounded exponent
+
+      add_and_round (q3, q4, e4, delta, p34, z_sign, p_sign, C3, C4,
+              rnd_mode, &is_midpoint_lt_even,
+              &is_midpoint_gt_even, &is_inexact_lt_midpoint,
+              &is_inexact_gt_midpoint, pfpsf, &res);
+      *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+      *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+      *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+      *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+      BID_SWAP128 (res);
+      BID_RETURN (res)
+
+    } else {
+      ;
+    }
+
+  } // end if delta < 0
+
+  *ptr_is_midpoint_lt_even = is_midpoint_lt_even;
+  *ptr_is_midpoint_gt_even = is_midpoint_gt_even;
+  *ptr_is_inexact_lt_midpoint = is_inexact_lt_midpoint;
+  *ptr_is_inexact_gt_midpoint = is_inexact_gt_midpoint;
+  BID_SWAP128 (res);
+  BID_RETURN (res)
+
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
+            _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+            _EXC_INFO_PARAM) {
+  UINT128 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128_fma (UINT128 x, UINT128 y, UINT128 z
+            _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+            _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even, is_midpoint_gt_even,
+    is_inexact_lt_midpoint, is_inexact_gt_midpoint;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, &x, &y, &z
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x, y,
+        		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128ddd_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT64 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128ddd_fma (UINT64 x, UINT64 y, UINT64 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 x1, y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, &x1, &y1, &z1
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x1, y1,
+        		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128ddq_fma (UINT128 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 x = *px, y = *py;
+  UINT128 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128ddq_fma (UINT64 x, UINT64 y, UINT128 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 x1, y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, &x1, &y1, &z
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x1, y1,
+        		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128dqd_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 x = *px, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128dqd_fma (UINT64 x, UINT128 y, UINT64 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 x1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, &x1, py, &z1
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x1, y,
+        		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128dqq_fma (UINT128 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 x = *px;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128dqq_fma (UINT64 x, UINT128 y, UINT128 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 x1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, &x1, py, pz
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x1, y,
+        		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qdd_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qdd_fma (UINT128 x, UINT64 y, UINT64 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, px, &y1, &z1
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x, y1,
+        		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qdq_fma (UINT128 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qdq_fma (UINT128 x, UINT64 y, UINT128 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, px, &y1, pz
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x, y1,
+        		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid128qqd_fma (UINT128 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+  UINT64 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT128
+bid128qqd_fma (UINT128 x, UINT128 y, UINT64 z
+               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+               _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        	  &is_inexact_lt_midpoint, &is_inexact_gt_midpoint,
+        	  &res, px, py, &z1
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res = bid128_ext_fma (&is_midpoint_lt_even, &is_midpoint_gt_even,
+        		&is_inexact_lt_midpoint,
+        		&is_inexact_gt_midpoint, x, y,
+        		z1 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+  BID_RETURN (res);
+}
+
+// Note: bid128qqq_fma is represented by bid128_fma
+
+// Note: bid64ddd_fma is represented by bid64_fma
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64ddq_fma (UINT64 * pres, UINT64 * px, UINT64 * py, UINT128 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 x = *px, y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64ddq_fma (UINT64 x, UINT64 y, UINT128 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 x1, y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, &x1, &y1, pz
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x1, y1, z
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64dqd_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT64 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 x = *px, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64dqd_fma (UINT64 x, UINT128 y, UINT64 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 x1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, &x1, py, &z1
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x1, y, z1
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64dqq_fma (UINT64 * pres, UINT64 * px, UINT128 * py, UINT128 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 x = *px;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64dqq_fma (UINT64 x, UINT128 y, UINT128 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 x1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, &x1, py, pz
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x1, y, z
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qdd_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT64 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qdd_fma (UINT128 x, UINT64 y, UINT64 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 y1, z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, px, &y1, &z1
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x, y1, z1
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qdq_fma (UINT64 * pres, UINT128 * px, UINT64 * py, UINT128 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 y = *py;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qdq_fma (UINT128 x, UINT64 y, UINT128 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 y1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, px, &y1, pz
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x, y1, z
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qqd_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT64 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT64 z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qqd_fma (UINT128 x, UINT128 y, UINT64 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  UINT128 z1;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid64_to_bid128 (&z1, &z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  bid64qqq_fma (&res1, px, py, &z1
+        	_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	_EXC_INFO_ARG);
+#else
+  z1 = bid64_to_bid128 (z _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
+  res1 = bid64qqq_fma (x, y, z1
+        	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	       _EXC_INFO_ARG);
+#endif
+  BID_RETURN (res1);
+}
+
+
+#if DECIMAL_CALL_BY_REFERENCE
+void
+bid64qqq_fma (UINT64 * pres, UINT128 * px, UINT128 * py, UINT128 * pz
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+  UINT128 x = *px, y = *py, z = *pz;
+#if !DECIMAL_GLOBAL_ROUNDING
+  unsigned int rnd_mode = *prnd_mode;
+#endif
+#else
+UINT64
+bid64qqq_fma (UINT128 x, UINT128 y, UINT128 z
+              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+              _EXC_INFO_PARAM) {
+#endif
+  int is_midpoint_lt_even0 = 0, is_midpoint_gt_even0 = 0,
+    is_inexact_lt_midpoint0 = 0, is_inexact_gt_midpoint0 = 0;
+  int is_midpoint_lt_even = 0, is_midpoint_gt_even = 0,
+    is_inexact_lt_midpoint = 0, is_inexact_gt_midpoint = 0;
+  int incr_exp;
+  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT128 res128 = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull} };
+  UINT64 res1 = 0xbaddbaddbaddbaddull;
+  unsigned int save_fpsf; // needed because of the call to bid128_ext_fma
+  UINT64 sign;
+  UINT64 exp;
+  int unbexp;
+  UINT128 C;
+  BID_UI64DOUBLE tmp;
+  int nr_bits;
+  int q, x0;
+  int scale;
+  int lt_half_ulp = 0, eq_half_ulp = 0;
+
+  // Note: for rounding modes other than RN or RA, the result can be obtained
+  // by rounding first to BID128 and then to BID64
+
+  save_fpsf = *pfpsf; // sticky bits - caller value must be preserved
+  *pfpsf = 0;
+
+#if DECIMAL_CALL_BY_REFERENCE
+  bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
+        	  &is_inexact_lt_midpoint0, &is_inexact_gt_midpoint0,
+        	  &res, &x, &y, &z
+        	  _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        	  _EXC_INFO_ARG);
+#else
+  res = bid128_ext_fma (&is_midpoint_lt_even0, &is_midpoint_gt_even0,
+        		&is_inexact_lt_midpoint0,
+        		&is_inexact_gt_midpoint0, x, y,
+        		z _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
+        		_EXC_INFO_ARG);
+#endif
+
+  if ((rnd_mode == ROUNDING_DOWN) || (rnd_mode == ROUNDING_UP) || 
+      (rnd_mode == ROUNDING_TO_ZERO) || // no double rounding error is possible
+      ((res.w[HIGH_128W] & MASK_NAN) == MASK_NAN) || //res=QNaN (cannot be SNaN)
+      ((res.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF)) { // result is infinity  
+#if DECIMAL_CALL_BY_REFERENCE
+    bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#else
+    res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#endif
+    // determine the unbiased exponent of the result
+    unbexp = ((res1 >> 53) & 0x3ff) - 398;
+
+    // if subnormal, res1  must have exp = -398
+    // if tiny and inexact set underflow and inexact status flags
+    if (!((res1 & MASK_NAN) == MASK_NAN) &&	// res1 not NaN
+        (unbexp == -398)
+        && ((res1 & MASK_BINARY_SIG1) < 1000000000000000ull)
+        && (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0
+            || is_midpoint_lt_even0 || is_midpoint_gt_even0)) {
+      // set the inexact flag and the underflow flag
+      *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
+    } else if (is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 ||
+               is_midpoint_lt_even0 || is_midpoint_gt_even0) {
+      // set the inexact flag and the underflow flag
+      *pfpsf |= INEXACT_EXCEPTION;
+    }
+
+    *pfpsf |= save_fpsf;
+    BID_RETURN (res1);
+  } // else continue, and use rounding to nearest to round to 16 digits
+
+  // at this point the result is rounded to nearest (even or away) to 34 digits
+  // (or less if exact), and it is zero or finite non-zero canonical [sub]normal
+  sign = res.w[HIGH_128W] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
+  exp = res.w[HIGH_128W] & MASK_EXP; // biased and shifted left 49 bits
+  unbexp = (exp >> 49) - 6176;
+  C.w[1] = res.w[HIGH_128W] & MASK_COEFF;
+  C.w[0] = res.w[LOW_128W];
+
+  if ((C.w[1] == 0x0 && C.w[0] == 0x0) ||	// result is zero
+      (unbexp <= (-398 - 35)) || (unbexp >= (369 + 16))) { 
+      // clear under/overflow
+#if DECIMAL_CALL_BY_REFERENCE
+    bid128_to_bid64 (&res1, &res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#else
+    res1 = bid128_to_bid64 (res _RND_MODE_ARG _EXC_FLAGS_ARG);
+#endif
+    *pfpsf |= save_fpsf;
+    BID_RETURN (res1);
+  } // else continue
+
+  // -398 - 34 <= unbexp <= 369 + 15
+  if (rnd_mode == ROUNDING_TIES_AWAY) {
+    // apply correction, if needed, to make the result rounded to nearest-even
+    if (is_midpoint_gt_even) {
+      // res = res - 1
+      res1--; // res1 is now even
+    } // else the result is already correctly rounded to nearest-even
+  }
+  // at this point the result is finite, non-zero canonical normal or subnormal,
+  // and in most cases overflow or underflow will not occur
+
+  // determine the number of digits q in the result
+  // q = nr. of decimal digits in x
+  // determine first the nr. of bits in x
+  if (C.w[1] == 0) {
+    if (C.w[0] >= 0x0020000000000000ull) { // x >= 2^53
+      // split the 64-bit value in two 32-bit halves to avoid rounding errors
+      if (C.w[0] >= 0x0000000100000000ull) { // x >= 2^32
+        tmp.d = (double) (C.w[0] >> 32); // exact conversion
+        nr_bits =
+          33 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+      } else { // x < 2^32
+        tmp.d = (double) (C.w[0]); // exact conversion
+        nr_bits =
+          1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+      }
+    } else { // if x < 2^53
+      tmp.d = (double) C.w[0]; // exact conversion
+      nr_bits =
+        1 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+    }
+  } else { // C.w[1] != 0 => nr. bits = 64 + nr_bits (C.w[1])
+    tmp.d = (double) C.w[1]; // exact conversion
+    nr_bits =
+      65 + ((((unsigned int) (tmp.ui64 >> 52)) & 0x7ff) - 0x3ff);
+  }
+  q = nr_digits[nr_bits - 1].digits;
+  if (q == 0) {
+    q = nr_digits[nr_bits - 1].digits1;
+    if (C.w[1] > nr_digits[nr_bits - 1].threshold_hi ||
+        (C.w[1] == nr_digits[nr_bits - 1].threshold_hi &&
+         C.w[0] >= nr_digits[nr_bits - 1].threshold_lo))
+      q++;
+  }
+  // if q > 16, round to nearest even to 16 digits (but for underflow it may 
+  // have to be truncated even more)
+  if (q > 16) {
+    x0 = q - 16;
+    if (q <= 18) {
+      round64_2_18 (q, x0, C.w[0], &res1, &incr_exp,
+        	    &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	    &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+    } else { // 19 <= q <= 34
+      round128_19_38 (q, x0, C, &res128, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+      res1 = res128.w[0]; // the result fits in 64 bits
+    }
+    unbexp = unbexp + x0;
+    if (incr_exp)
+      unbexp++;
+    q = 16; // need to set in case denormalization is necessary
+  } else {
+    // the result does not require a second rounding (and it must have 
+    // been exact in the first rounding, since q <= 16)
+    res1 = C.w[0];
+  }
+
+  // avoid a double rounding error
+  if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+      is_midpoint_lt_even) { // double rounding error upward
+    // res = res - 1 
+    res1--; // res1 becomes odd 
+    is_midpoint_lt_even = 0;
+    is_inexact_lt_midpoint = 1;
+    if (res1 == 0x00038d7ea4c67fffull) { // 10^15 - 1
+      res1 = 0x002386f26fc0ffffull; // 10^16 - 1 
+      unbexp--;
+    }
+  } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+      is_midpoint_gt_even) { // double rounding error downward
+    // res = res + 1
+    res1++; // res1 becomes odd (so it cannot be 10^16)
+    is_midpoint_gt_even = 0;
+    is_inexact_gt_midpoint = 1;
+  } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+             !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+    // if this second rounding was exact the result may still be 
+    // inexact because of the first rounding
+    if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+      is_inexact_gt_midpoint = 1;
+    }
+    if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+      is_inexact_lt_midpoint = 1;
+    }
+  } else if (is_midpoint_gt_even &&
+             (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+    // pulled up to a midpoint 
+    is_inexact_lt_midpoint = 1;
+    is_inexact_gt_midpoint = 0;
+    is_midpoint_lt_even = 0;
+    is_midpoint_gt_even = 0;
+  } else if (is_midpoint_lt_even &&
+             (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+    // pulled down to a midpoint 
+    is_inexact_lt_midpoint = 0;
+    is_inexact_gt_midpoint = 1;
+    is_midpoint_lt_even = 0;
+    is_midpoint_gt_even = 0;
+  } else {
+    ;
+  }
+  // this is the result rounded correctly to nearest even, with unbounded exp. 
+
+  // check for overflow
+  if (q + unbexp > P16 + expmax16) {
+    res1 = sign | 0x7800000000000000ull;
+    *pfpsf |= (INEXACT_EXCEPTION | OVERFLOW_EXCEPTION);
+    *pfpsf |= save_fpsf;
+    BID_RETURN (res1)
+  } else if (unbexp > expmax16) { // q + unbexp <= P16 + expmax16
+    // not overflow; the result must be exact, and we can multiply res1 by
+    // 10^(unbexp - expmax16) and the product will fit in 16 decimal digits
+    scale = unbexp - expmax16;
+    res1 = res1 * ten2k64[scale]; // res1 * 10^scale
+    unbexp = expmax16; // unbexp - scale 
+  } else {
+    ; // continue
+  }
+
+  // check for underflow
+  if (q + unbexp < P16 + expmin16) {
+    if (unbexp < expmin16) {
+      // we must truncate more of res
+      x0 = expmin16 - unbexp; // x0 >= 1
+      is_inexact_lt_midpoint0 = is_inexact_lt_midpoint;
+      is_inexact_gt_midpoint0 = is_inexact_gt_midpoint;
+      is_midpoint_lt_even0 = is_midpoint_lt_even;
+      is_midpoint_gt_even0 = is_midpoint_gt_even;
+      is_inexact_lt_midpoint = 0;
+      is_inexact_gt_midpoint = 0;
+      is_midpoint_lt_even = 0;
+      is_midpoint_gt_even = 0;
+      // the number of decimal digits in res1 is q
+      if (x0 < q) { // 1 <= x0 <= q-1 => round res to q - x0 digits
+        // 2 <= q <= 16, 1 <= x0 <= 15
+        round64_2_18 (q, x0, res1, &res1, &incr_exp,
+        	      &is_midpoint_lt_even, &is_midpoint_gt_even,
+        	      &is_inexact_lt_midpoint, &is_inexact_gt_midpoint);
+        if (incr_exp) {
+          // res1 = 10^(q-x0), 1 <= q - x0 <= q - 1, 1 <= q - x0 <= 15
+          res1 = ten2k64[q - x0];
+        }
+        unbexp = unbexp + x0; // expmin16
+      } else if (x0 == q) {
+        // the second rounding is for 0.d(0)d(1)...d(q-1) * 10^emin
+        // determine relationship with 1/2 ulp
+        // q <= 16
+        if (res1 < midpoint64[q - 1]) { // < 1/2 ulp
+          lt_half_ulp = 1;
+          is_inexact_lt_midpoint = 1;
+        } else if (res1 == midpoint64[q - 1]) { // = 1/2 ulp
+          eq_half_ulp = 1;
+          is_midpoint_gt_even = 1;
+        } else { // > 1/2 ulp
+          // gt_half_ulp = 1;
+          is_inexact_gt_midpoint = 1;
+        }
+        if (lt_half_ulp || eq_half_ulp) {
+          // res = +0.0 * 10^expmin16
+          res1 = 0x0000000000000000ull;
+        } else { // if (gt_half_ulp)
+          // res = +1 * 10^expmin16
+          res1 = 0x0000000000000001ull;
+        }
+        unbexp = expmin16;
+      } else { // if (x0 > q)
+        // the second rounding is for 0.0...d(0)d(1)...d(q-1) * 10^emin
+        res1 = 0x0000000000000000ull;
+        unbexp = expmin16;
+        is_inexact_lt_midpoint = 1;
+      }
+      // avoid a double rounding error
+      if ((is_inexact_gt_midpoint0 || is_midpoint_lt_even0) && 
+          is_midpoint_lt_even) { // double rounding error upward
+        // res = res - 1
+        res1--; // res1 becomes odd
+        is_midpoint_lt_even = 0;
+        is_inexact_lt_midpoint = 1;
+      } else if ((is_inexact_lt_midpoint0 || is_midpoint_gt_even0) && 
+          is_midpoint_gt_even) { // double rounding error downward
+        // res = res + 1
+        res1++; // res1 becomes odd
+        is_midpoint_gt_even = 0;
+        is_inexact_gt_midpoint = 1;
+      } else if (!is_midpoint_lt_even && !is_midpoint_gt_even &&
+        	 !is_inexact_lt_midpoint && !is_inexact_gt_midpoint) {
+        // if this rounding was exact the result may still be 
+        // inexact because of the previous roundings
+        if (is_inexact_gt_midpoint0 || is_midpoint_lt_even0) {
+          is_inexact_gt_midpoint = 1;
+        }
+        if (is_inexact_lt_midpoint0 || is_midpoint_gt_even0) {
+          is_inexact_lt_midpoint = 1;
+        }
+      } else if (is_midpoint_gt_even &&
+        	 (is_inexact_gt_midpoint0 || is_midpoint_lt_even0)) {
+        // pulled up to a midpoint
+        is_inexact_lt_midpoint = 1;
+        is_inexact_gt_midpoint = 0;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+      } else if (is_midpoint_lt_even &&
+        	 (is_inexact_lt_midpoint0 || is_midpoint_gt_even0)) {
+        // pulled down to a midpoint
+        is_inexact_lt_midpoint = 0;
+        is_inexact_gt_midpoint = 1;
+        is_midpoint_lt_even = 0;
+        is_midpoint_gt_even = 0;
+      } else {
+        ;
+      }
+    }
+    // else if unbexp >= emin then q < P (because q + unbexp < P16 + expmin16)
+    // and the result is tiny and exact
+
+    // check for inexact result
+    if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+        is_midpoint_lt_even || is_midpoint_gt_even ||
+        is_inexact_lt_midpoint0 || is_inexact_gt_midpoint0 ||
+        is_midpoint_lt_even0 || is_midpoint_gt_even0) {
+      // set the inexact flag and the underflow flag
+      *pfpsf |= (INEXACT_EXCEPTION | UNDERFLOW_EXCEPTION);
+    }
+  } else if (is_inexact_lt_midpoint || is_inexact_gt_midpoint ||
+             is_midpoint_lt_even || is_midpoint_gt_even) {
+    *pfpsf |= INEXACT_EXCEPTION;
+  }
+  // this is the result rounded correctly to nearest, with bounded exponent
+
+  if (rnd_mode == ROUNDING_TIES_AWAY && is_midpoint_gt_even) { // correction
+    // res = res + 1
+    res1++; // res1 is now odd
+  } // else the result is already correct
+
+  // assemble the result
+  if (res1 < 0x0020000000000000ull) { // res < 2^53
+    res1 = sign | ((UINT64) (unbexp + 398) << 53) | res1;
+  } else { // res1 >= 2^53
+    res1 = sign | MASK_STEERING_BITS |
+      ((UINT64) (unbexp + 398) << 51) | (res1 & MASK_BINARY_SIG2);
+  }
+  *pfpsf |= save_fpsf;
+  BID_RETURN (res1);
+}