Merge latest changes from http://www.netlib.org/fp/dtoa.c into dtoa.cc

base/third_party/dmg_fp/dtoa.cc

 /****************************************************************
  *
  * The author of this software is David M. Gay.
  *
  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
  *
  * Permission to use, copy, modify, and distribute this software for any
  * purpose without fee is hereby granted, provided that this entire notice
  * is included in all copies of any software which is or includes a copy
  * or modification of this software and in all copies of the supporting
@@ skipping 14 matching lines @@
  * before invoking strtod or dtoa.  If the machine uses (the equivalent
  * of) Intel 80x87 arithmetic, the call
  *      _control87(PC_53, MCW_PC);
  * does this with many compilers.  Whether this or another call is
  * appropriate depends on the compiler; for this to work, it may be
  * necessary to #include "float.h" or another system-dependent header
  * file.
  */
 
 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
+ * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.)
  *
  * This strtod returns a nearest machine number to the input decimal
  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  * biased rounding (add half and chop).
  *
  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  *
  * Modifications:
@@ skipping 18 matching lines @@
 
 /*
  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  *      significant byte has the lowest address.
  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  *      significant byte has the lowest address.
  * #define Long int on machines with 32-bit ints and 64-bit longs.
  * #define IBM for IBM mainframe-style floating-point arithmetic.
  * #define VAX for VAX-style floating-point arithmetic (D_floating).
  * #define No_leftright to omit left-right logic in fast floating-point
- *      computation of dtoa.
+ *      computation of dtoa.  This will cause dtoa modes 4 and 5 to be
+ *      treated the same as modes 2 and 3 for some inputs.
  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  *      and strtod and dtoa should round accordingly.  Unless Trust_FLT_ROUNDS
  *      is also #defined, fegetround() will be queried for the rounding mode.
  *      Note that both FLT_ROUNDS and fegetround() are specified by the C99
  *      standard (and are specified to be consistent, with fesetround()
  *      affecting the value of FLT_ROUNDS), but that some (Linux) systems
  *      do not work correctly in this regard, so using fegetround() is more
- *      portable than using FLT_FOUNDS directly.
+ *      portable than using FLT_ROUNDS directly.
  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  *      and Honor_FLT_ROUNDS is not #defined.
  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  *      that use extended-precision instructions to compute rounded
  *      products and quotients) with IBM.
- * #define ROUND_BIASED for IEEE-format with biased rounding.
+ * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
+ *      that rounds toward +Infinity.
+ * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
+ *      rounding when the underlying floating-point arithmetic uses
+ *      unbiased rounding.  This prevent using ordinary floating-point
+ *      arithmetic when the result could be computed with one rounding error.
  * #define Inaccurate_Divide for IEEE-format with correctly rounded
  *      products but inaccurate quotients, e.g., for Intel i860.
  * #define NO_LONG_LONG on machines that do not have a "long long"
  *      integer type (of >= 64 bits).  On such machines, you can
  *      #define Just_16 to store 16 bits per 32-bit Long when doing
  *      high-precision integer arithmetic.  Whether this speeds things
  *      up or slows things down depends on the machine and the number
  *      being converted.  If long long is available and the name is
  *      something other than "long long", #define Llong to be the name,
  *      and if "unsigned Llong" does not work as an unsigned version of
@@ skipping 353 matching lines @@
 #define Log2P 1
 #define Tiny0 0x80
 #define Tiny1 0
 #define Quick_max 15
 #define Int_max 15
 #endif /* IBM, VAX */
 #endif /* IEEE_Arith */
 
 #ifndef IEEE_Arith
 #define ROUND_BIASED
+#else
+#ifdef ROUND_BIASED_without_Round_Up
+#undef  ROUND_BIASED
+#define ROUND_BIASED
+#endif
 #endif
 
 #ifdef RND_PRODQUOT
 #define rounded_product(a,b) a = rnd_prod(a, b)
 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 #ifdef KR_headers
 extern double rnd_prod(), rnd_quot();
 #else
 extern double rnd_prod(double, double), rnd_quot(double, double);
 #endif
@@ skipping 185 matching lines @@
                 b->wds = wds;
                 }
         return b;
         }
 
  static Bigint *
 s2b
 #ifdef KR_headers
         (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
 #else
-        (CONST char *s, int nd0, int nd, ULong y9, int dplen)
+        (const char *s, int nd0, int nd, ULong y9, int dplen)
 #endif
 {
         Bigint *b;
         int i, k;
         Long x, y;
 
         x = (nd + 8) / 9;
         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
 #ifdef Pack_32
         b = Balloc(k);
@@ skipping 816 matching lines @@
 #endif
 #endif
 
 #ifndef Need_Hexdig
 #ifndef NO_HEX_FP
 #define Need_Hexdig
 #endif
 #endif
 
 #ifdef Need_Hexdig /*{*/
+#if 0
 static unsigned char hexdig[256];
 
  static void
-#ifdef KR_headers
-htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
-#else
 htinit(unsigned char *h, unsigned char *s, int inc)
-#endif
 {
         int i, j;
         for(i = 0; (j = s[i]) !=0; i++)
                 h[j] = (unsigned char)(i + inc);
         }
 
  static void
-#ifdef KR_headers
-hexdig_init()
-#else
-hexdig_init(void)
-#endif
+hexdig_init(void)       /* Use of hexdig_init omitted 20121220 to avoid a */
+                        /* race condition when multiple threads are used. */
 {
 #define USC (unsigned char *)
         htinit(hexdig, USC "0123456789", 0x10);
         htinit(hexdig, USC "abcdef", 0x10 + 10);
         htinit(hexdig, USC "ABCDEF", 0x10 + 10);
         }
+#else
+static unsigned char hexdig[256] = {

[Lei Zhang, 2016/09/23 19:32:16]

Should this be const?

[kcwu, 2016/09/23 19:48:03]

Done.

[Lei Zhang, 2016/09/23 20:33:12]

Thanks, but now we need another .patch file.

[kcwu, 2016/09/24 03:18:01]

Done.

+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
+        0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+        };
+#endif
 #endif /* } Need_Hexdig */
 
 #ifdef INFNAN_CHECK
 
 #ifndef NAN_WORD0
 #define NAN_WORD0 0x7ff80000
 #endif
 
 #ifndef NAN_WORD1
 #define NAN_WORD1 0
 #endif
 
  static int
 match
 #ifdef KR_headers
         (sp, t) char **sp, *t;
 #else
-        (CONST char **sp, CONST char *t)
+        (const char **sp, const char *t)
 #endif
 {
         int c, d;
         CONST char *s = *sp;
 
         for(d = *t++; d; d = *t++) {
                 if ((c = *++s) >= 'A' && c <= 'Z')
                         c += 'a' - 'A';
                 if (c != d)
                         return 0;
                 }
         *sp = s + 1;
         return 1;
         }
 
 #ifndef No_Hex_NaN
  static void
 hexnan
 #ifdef KR_headers
         (rvp, sp) U *rvp; CONST char **sp;
 #else
-        (U *rvp, CONST char **sp)
+        (U *rvp, const char **sp)
 #endif
 {
         ULong c, x[2];
         CONST char *s;
         int c1, havedig, udx0, xshift;
 
-        if (!hexdig['0'])
-                hexdig_init();
+        /**** if (!hexdig['0']) hexdig_init(); ****/

[Lei Zhang, 2016/09/23 19:32:16]

Hmm

[kcwu, 2016/09/23 19:48:03]

This is copied from upstream. I'd like to keep it in order to make future merge
easier.

         x[0] = x[1] = 0;
         havedig = xshift = 0;
         udx0 = 1;
         s = *sp;
         /* allow optional initial 0x or 0X */
         for(c = *(CONST unsigned char*)(s+1); c && c <= ' '; c = *(CONST unsigned char*)(s+1))
                 ++s;
         if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
                 s += 2;
         for(c = *(CONST unsigned char*)++s; c; c = *(CONST unsigned char*)++s) {
@@ skipping 46 matching lines @@
 
 #ifdef Pack_32
 #define ULbits 32
 #define kshift 5
 #define kmask 31
 #else
 #define ULbits 16
 #define kshift 4
 #define kmask 15
 #endif
+
+#if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
+ static Bigint *
+#ifdef KR_headers
+increment(b) Bigint *b;
+#else
+increment(Bigint *b)
+#endif
+{
+        ULong *x, *xe;
+        Bigint *b1;
+
+        x = b->x;
+        xe = x + b->wds;
+        do {
+                if (*x < (ULong)0xffffffffL) {
+                        ++*x;
+                        return b;
+                        }
+                *x++ = 0;
+                } while(x < xe);
+        {
+                if (b->wds >= b->maxwds) {
+                        b1 = Balloc(b->k+1);
+                        Bcopy(b1,b);
+                        Bfree(b);
+                        b = b1;
+                        }
+                b->x[b->wds++] = 1;
+                }
+        return b;
+        }
+
+#endif /*}*/
+
 #ifndef NO_HEX_FP /*{*/
 
  static void
 #ifdef KR_headers
 rshift(b, k) Bigint *b; int k;
 #else
 rshift(Bigint *b, int k)
 #endif
 {
         ULong *x, *x1, *xe, y;
@@ skipping 52 matching lines @@
         return 0;
         }
 
 enum {  /* rounding values: same as FLT_ROUNDS */
         Round_zero = 0,
         Round_near = 1,
         Round_up = 2,
         Round_down = 3
         };
 
- static Bigint *
-#ifdef KR_headers
-increment(b) Bigint *b;
-#else
-increment(Bigint *b)
-#endif
-{
-        ULong *x, *xe;
-        Bigint *b1;
-
-        x = b->x;
-        xe = x + b->wds;
-        do {
-                if (*x < (ULong)0xffffffffL) {
-                        ++*x;
-                        return b;
-                        }
-                *x++ = 0;
-                } while(x < xe);
-        {
-                if (b->wds >= b->maxwds) {
-                        b1 = Balloc(b->k+1);
-                        Bcopy(b1,b);
-                        Bfree(b);
-                        b = b1;
-                        }
-                b->x[b->wds++] = 1;
-                }
-        return b;
-        }
-
  void
 #ifdef KR_headers
 gethex(sp, rvp, rounding, sign)
         CONST char **sp; U *rvp; int rounding, sign;
 #else
 gethex( CONST char **sp, U *rvp, int rounding, int sign)
 #endif
 {
         Bigint *b;
         CONST unsigned char *decpt, *s0, *s, *s1;
@@ skipping 30 matching lines @@
                 if ((decimalpoint_cache = (unsigned char*)
                                 MALLOC(strlen((CONST char*)s0) + 1))) {
                         strcpy((char*)decimalpoint_cache, (CONST char*)s0);
                         s0 = decimalpoint_cache;
                         }
                 }
         decimalpoint = s0;
 #endif
 #endif
 
-        if (!hexdig['0'])
-                hexdig_init();
+        /**** if (!hexdig['0']) hexdig_init(); ****/
         havedig = 0;
         s0 = *(CONST unsigned char **)sp + 2;
         while(s0[havedig] == '0')
                 havedig++;
         s0 += havedig;
         s = s0;
         decpt = 0;
         zret = 0;
         e = 0;
         if (hexdig[*s])
@@ skipping 79 matching lines @@
                                         break;
                                 goto ret_tiny;
                           case Round_down:
                                 if (!sign)
                                         break;
                                 goto ret_tiny;
                           }
 #endif
                         goto retz;
 #ifdef IEEE_Arith
+ ret_tinyf:
+                        Bfree(b);
  ret_tiny:
 #ifndef NO_ERRNO
                         errno = ERANGE;
 #endif
                         word0(rvp) = 0;
                         word1(rvp) = 1;
                         return;
 #endif /* IEEE_Arith */
                         }
                 switch(rounding) {
@@ skipping 80 matching lines @@
                 }
         denorm = 0;
         if (e < emin) {
                 denorm = 1;
                 n = emin - e;
                 if (n >= nbits) {
 #ifdef IEEE_Arith /*{*/
                         switch (rounding) {
                           case Round_near:
                                 if (n == nbits && (n < 2 || any_on(b,n-1)))
-                                        goto ret_tiny;
+                                        goto ret_tinyf;
                                 break;
                           case Round_up:
                                 if (!sign)
-                                        goto ret_tiny;
+                                        goto ret_tinyf;
                                 break;
                           case Round_down:
                                 if (sign)
-                                        goto ret_tiny;
+                                        goto ret_tinyf;
                           }
 #endif /* } IEEE_Arith */
                         Bfree(b);
  retz:
 #ifndef NO_ERRNO
                         errno = ERANGE;
 #endif
  retz1:
                         rvp->d = 0.;
                         return;
@@ skipping 79 matching lines @@
 #endif
 #ifdef VAX
         /* The next two lines ignore swap of low- and high-order 2 bytes. */
         /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
         /* word1(rvp) = b->x[0]; */
         word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
         word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
 #endif
         Bfree(b);
         }
-#endif /*}!NO_HEX_FP*/
+#endif /*!NO_HEX_FP}*/
 
  static int
 #ifdef KR_headers
 dshift(b, p2) Bigint *b; int p2;
 #else
 dshift(Bigint *b, int p2)
 #endif
 {
         int rv = hi0bits(b->x[b->wds-1]) - 4;
         if (p2 > 0)
@@ skipping 26 matching lines @@
         /*debug*/       Bug("oversize b in quorem");
 #endif
         if (b->wds < n)
                 return 0;
         sx = S->x;
         sxe = sx + --n;
         bx = b->x;
         bxe = bx + n;
         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
 #ifdef DEBUG
+#ifdef NO_STRTOD_BIGCOMP
         /*debug*/ if (q > 9)
+#else
+        /* An oversized q is possible when quorem is called from bigcomp and */
+        /* the input is near, e.g., twice the smallest denormalized number. */
+        /*debug*/ if (q > 15)
+#endif
         /*debug*/       Bug("oversized quotient in quorem");
 #endif
         if (q) {
                 borrow = 0;
                 carry = 0;
                 do {
 #ifdef ULLong
                         ys = *sx++ * (ULLong)q + carry;
                         carry = ys >> 32;
                         y = *bx - (ys & FFFFFFFF) - borrow;
@@ skipping 65 matching lines @@
                 bxe = bx + n;
                 if (!*bxe) {
                         while(--bxe > bx && !*bxe)
                                 --n;
                         b->wds = n;
                         }
                 }
         return q;
         }
 
+#if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
+ static double
+sulp
+#ifdef KR_headers
+        (x, bc) U *x; BCinfo *bc;
+#else
+        (U *x, BCinfo *bc)
+#endif
+{
+        U u;
+        double rv;
+        int i;
+
+        rv = ulp(x);
+        if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
+                return rv; /* Is there an example where i <= 0 ? */
+        word0(&u) = Exp_1 + (i << Exp_shift);
+        word1(&u) = 0;
+        return rv * u.d;
+        }
+#endif /*}*/
+
 #ifndef NO_STRTOD_BIGCOMP
-
  static void
 bigcomp
 #ifdef KR_headers
         (rv, s0, bc)
         U *rv; CONST char *s0; BCinfo *bc;
 #else
-        (U *rv, CONST char *s0, BCinfo *bc)
+        (U *rv, const char *s0, BCinfo *bc)
 #endif
 {
         Bigint *b, *d;
         int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
 
         dsign = bc->dsign;
         nd = bc->nd;
         nd0 = bc->nd0;
         p5 = nd + bc->e0 - 1;
         dd = speccase = 0;
@@ skipping 111 matching lines @@
                 if (dd)
                         goto ret;
                 if (!b->x[0] && b->wds == 1) {
                         if (i < nd)
                                 dd = 1;
                         goto ret;
                         }
                 b = multadd(b, 10, 0);
                 dig = quorem(b,d);
                 }
-        if (b->x[0] || b->wds > 1)
+        if (dig > 0 || b->x[0] || b->wds > 1)
                 dd = -1;
  ret:
         Bfree(b);
         Bfree(d);
 #ifdef Honor_FLT_ROUNDS
         if (bc->rounding != 1) {
                 if (dd < 0) {
                         if (bc->rounding == 0) {
                                 if (!dsign)
                                         goto retlow1;
                                 }
                         else if (dsign)
                                 goto rethi1;
                         }
                 else if (dd > 0) {
                         if (bc->rounding == 0) {
                                 if (dsign)
                                         goto rethi1;
                                 goto ret1;
                                 }
                         if (!dsign)
                                 goto rethi1;
-                        dval(rv) += 2.*ulp(rv);
+                        dval(rv) += 2.*sulp(rv,bc);
                         }
                 else {
                         bc->inexact = 0;
                         if (dsign)
                                 goto rethi1;
                         }
                 }
         else
 #endif
         if (speccase) {
                 if (dd <= 0)
                         rv->d = 0.;
                 }
         else if (dd < 0) {
                 if (!dsign)     /* does not happen for round-near */
 retlow1:
-                        dval(rv) -= ulp(rv);
+                        dval(rv) -= sulp(rv,bc);
                 }
         else if (dd > 0) {
                 if (dsign) {
  rethi1:
-                        dval(rv) += ulp(rv);
+                        dval(rv) += sulp(rv,bc);
                         }
                 }
         else {
                 /* Exact half-way case:  apply round-even rule. */
-                if (word1(rv) & 1) {
+                if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
+                        i = 1 - j;
+                        if (i <= 31) {
+                                if (word1(rv) & (0x1 << i))
+                                        goto odd;
+                                }
+                        else if (word0(rv) & (0x1 << (i-32)))
+                                goto odd;
+                        }
+                else if (word1(rv) & 1) {
+ odd:
                         if (dsign)
                                 goto rethi1;
                         goto retlow1;
                         }
                 }
 
 #ifdef Honor_FLT_ROUNDS
  ret1:
 #endif
         return;
         }
 #endif /* NO_STRTOD_BIGCOMP */
 
  double
 strtod
 #ifdef KR_headers
         (s00, se) CONST char *s00; char **se;
 #else
-        (CONST char *s00, char **se)
+        (const char *s00, char **se)
 #endif
 {
         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
-        int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
+        int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
         CONST char *s, *s0, *s1;
         double aadj, aadj1;
         Long L;
         U aadj2, adj, rv, rv0;
         ULong y, z;
         BCinfo bc;
-        Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
+        Bigint *bb = nullptr, *bb1, *bd = nullptr, *bd0, *bs = nullptr, *delta = nullptr;
+#ifdef Avoid_Underflow
+        ULong Lsb, Lsb1;
+#endif
 #ifdef SET_INEXACT
         int oldinexact;
 #endif
+#ifndef NO_STRTOD_BIGCOMP
+        int req_bigcomp = 0;
+#endif
 #ifdef Honor_FLT_ROUNDS /*{*/
 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
         bc.rounding = Flt_Rounds;
 #else /*}{*/
         bc.rounding = 1;
         switch(fegetround()) {
           case FE_TOWARDZERO:   bc.rounding = 0; break;
           case FE_UPWARD:       bc.rounding = 2; break;
           case FE_DOWNWARD:     bc.rounding = 3;
           }
 #endif /*}}*/
 #endif /*}*/
 #ifdef USE_LOCALE
         CONST char *s2;
 #endif
 
-        sign = nz0 = nz = bc.dplen = bc.uflchk = 0;
+        sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
         dval(&rv) = 0.;
         for(s = s00;;s++) switch(*s) {
                 case '-':
                         sign = 1;
                         /* no break */
                 case '+':
                         if (*++s)
                                 goto break2;
                         /* no break */
                 case 0:
@@ skipping 25 matching lines @@
                 nz0 = 1;
                 while(*++s == '0') ;
                 if (!*s)
                         goto ret;
                 }
         s0 = s;
         y = z = 0;
         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
                 if (nd < 9)
                         y = 10*y + c - '0';
-                else if (nd < 16)
+                else if (nd < DBL_DIG + 2)
                         z = 10*z + c - '0';
         nd0 = nd;
         bc.dp0 = bc.dp1 = s - s0;
+        for(s1 = s; s1 > s0 && *--s1 == '0'; )
+                ++nz1;
 #ifdef USE_LOCALE
         s1 = localeconv()->decimal_point;
         if (c == *s1) {
                 c = '.';
                 if (*++s1) {
                         s2 = s;
                         for(;;) {
                                 if (*++s2 != *s1) {
                                         c = 0;
                                         break;
                                         }
                                 if (!*++s1) {
                                         s = s2;
                                         break;
                                         }
                                 }
                         }
                 }
 #endif
         if (c == '.') {
                 c = *++s;
                 bc.dp1 = s - s0;
                 bc.dplen = bc.dp1 - bc.dp0;
                 if (!nd) {
                         for(; c == '0'; c = *++s)
                                 nz++;
                         if (c > '0' && c <= '9') {
+                                bc.dp0 = s0 - s;
+                                bc.dp1 = bc.dp0 + bc.dplen;
                                 s0 = s;
                                 nf += nz;
                                 nz = 0;
                                 goto have_dig;
                                 }
                         goto dig_done;
                         }
                 for(; c >= '0' && c <= '9'; c = *++s) {
  have_dig:
                         nz++;
                         if (c -= '0') {
                                 nf += nz;
                                 for(i = 1; i < nz; i++)
                                         if (nd++ < 9)
                                                 y *= 10;
-                                        else if (nd <= DBL_DIG + 1)
+                                        else if (nd <= DBL_DIG + 2)
                                                 z *= 10;
                                 if (nd++ < 9)
                                         y = 10*y + c;
-                                else if (nd <= DBL_DIG + 1)
+                                else if (nd <= DBL_DIG + 2)
                                         z = 10*z + c;
-                                nz = 0;
+                                nz = nz1 = 0;
                                 }
                         }
                 }
  dig_done:
         e = 0;
         if (c == 'e' || c == 'E') {
                 if (!nd && !nz && !nz0) {
                         goto ret0;
                         }
                 s00 = s;
@@ skipping 69 matching lines @@
                 }
         bc.e0 = e1 = e -= nf;
 
         /* Now we have nd0 digits, starting at s0, followed by a
          * decimal point, followed by nd-nd0 digits.  The number we're
          * after is the integer represented by those digits times
          * 10**e */
 
         if (!nd0)
                 nd0 = nd;
-        k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
+        k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
         dval(&rv) = y;
         if (k > 9) {
 #ifdef SET_INEXACT
                 if (k > DBL_DIG)
                         oldinexact = get_inexact();
 #endif
                 dval(&rv) = tens[k - 9] * dval(&rv) + z;
                 }
         bd0 = 0;
         if (nd <= DBL_DIG
 #ifndef RND_PRODQUOT
 #ifndef Honor_FLT_ROUNDS
                 && Flt_Rounds == 1
 #endif
 #endif
                         ) {
                 if (!e)
                         goto ret;
+#ifndef ROUND_BIASED_without_Round_Up
                 if (e > 0) {
                         if (e <= Ten_pmax) {
 #ifdef VAX
                                 goto vax_ovfl_check;
 #else
 #ifdef Honor_FLT_ROUNDS
                                 /* round correctly FLT_ROUNDS = 2 or 3 */
                                 if (sign) {
                                         rv.d = -rv.d;
                                         sign = 0;
@@ skipping 40 matching lines @@
                         /* round correctly FLT_ROUNDS = 2 or 3 */
                         if (sign) {
                                 rv.d = -rv.d;
                                 sign = 0;
                                 }
 #endif
                         /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
                         goto ret;
                         }
 #endif
+#endif /* ROUND_BIASED_without_Round_Up */
                 }
         e1 += nd - k;
 
 #ifdef IEEE_Arith
 #ifdef SET_INEXACT
         bc.inexact = 1;
         if (k <= DBL_DIG)
                 oldinexact = get_inexact();
 #endif
 #ifdef Avoid_Underflow
@@ skipping 12 matching lines @@
 
         /* Get starting approximation = rv * 10**e1 */
 
         if (e1 > 0) {
                 i = e1 & 15;
                 if (i)
                         dval(&rv) *= tens[i];
                 if (e1 &= ~15) {
                         if (e1 > DBL_MAX_10_EXP) {
  ovfl:
-#ifndef NO_ERRNO
-                                errno = ERANGE;
-#endif
                                 /* Can't trust HUGE_VAL */
 #ifdef IEEE_Arith
 #ifdef Honor_FLT_ROUNDS
                                 switch(bc.rounding) {
                                   case 0: /* toward 0 */
                                   case 3: /* toward -infinity */
                                         word0(&rv) = Big0;
                                         word1(&rv) = Big1;
                                         break;
                                   default:
                                         word0(&rv) = Exp_mask;
                                         word1(&rv) = 0;
                                   }
 #else /*Honor_FLT_ROUNDS*/
                                 word0(&rv) = Exp_mask;
                                 word1(&rv) = 0;
 #endif /*Honor_FLT_ROUNDS*/
 #ifdef SET_INEXACT
                                 /* set overflow bit */
                                 dval(&rv0) = 1e300;
                                 dval(&rv0) *= dval(&rv0);
 #endif
 #else /*IEEE_Arith*/
                                 word0(&rv) = Big0;
                                 word1(&rv) = Big1;
 #endif /*IEEE_Arith*/
+ range_err:
+                                if (bd0) {
+                                        Bfree(bb);
+                                        Bfree(bd);
+                                        Bfree(bs);
+                                        Bfree(bd0);
+                                        Bfree(delta);
+                                        }
+#ifndef NO_ERRNO
+                                errno = ERANGE;
+#endif
                                 goto ret;
                                 }
                         e1 >>= 4;
                         for(j = 0; e1 > 1; j++, e1 >>= 1)
                                 if (e1 & 1)
                                         dval(&rv) *= bigtens[j];
                 /* The last multiplication could overflow. */
                         word0(&rv) -= P*Exp_msk1;
                         dval(&rv) *= bigtens[j];
                         if ((z = word0(&rv) & Exp_mask)
@@ skipping 20 matching lines @@
 #ifdef Avoid_Underflow
                         if (e1 & Scale_Bit)
                                 bc.scale = 2*P;
                         for(j = 0; e1 > 0; j++, e1 >>= 1)
                                 if (e1 & 1)
                                         dval(&rv) *= tinytens[j];
                         if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
                                                 >> Exp_shift)) > 0) {
                                 /* scaled rv is denormal; clear j low bits */
                                 if (j >= 32) {
+                                        if (j > 54)
+                                                goto undfl;
                                         word1(&rv) = 0;
                                         if (j >= 53)
                                          word0(&rv) = (P+2)*Exp_msk1;
                                         else
                                          word0(&rv) &= 0xffffffff << (j-32);
                                         }
                                 else
                                         word1(&rv) &= 0xffffffff << j;
                                 }
 #else
                         for(j = 0; e1 > 1; j++, e1 >>= 1)
                                 if (e1 & 1)
                                         dval(&rv) *= tinytens[j];
                         /* The last multiplication could underflow. */
                         dval(&rv0) = dval(&rv);
                         dval(&rv) *= tinytens[j];
                         if (!dval(&rv)) {
                                 dval(&rv) = 2.*dval(&rv0);
                                 dval(&rv) *= tinytens[j];
 #endif
                                 if (!dval(&rv)) {
  undfl:
                                         dval(&rv) = 0.;
-#ifndef NO_ERRNO
-                                        errno = ERANGE;
-#endif
-                                        goto ret;
+                                        goto range_err;
                                         }
 #ifndef Avoid_Underflow
                                 word0(&rv) = Tiny0;
                                 word1(&rv) = Tiny1;
                                 /* The refinement below will clean
                                  * this approximation up.
                                  */
                                 }
 #endif
                         }
                 }
 
         /* Now the hard part -- adjusting rv to the correct value.*/
 
         /* Put digits into bd: true value = bd * 10^e */
 
-        bc.nd = nd;
+        bc.nd = nd - nz1;
 #ifndef NO_STRTOD_BIGCOMP
         bc.nd0 = nd0;   /* Only needed if nd > strtod_diglim, but done here */
                         /* to silence an erroneous warning about bc.nd0 */
                         /* possibly not being initialized. */
         if (nd > strtod_diglim) {
                 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
                 /* minimum number of decimal digits to distinguish double values */
                 /* in IEEE arithmetic. */
                 i = j = 18;
                 if (i > nd0)
                         j += bc.dplen;
                 for(;;) {
-                        if (--j <= bc.dp1 && j >= bc.dp0)
+                        if (--j < bc.dp1 && j >= bc.dp0)
                                 j = bc.dp0 - 1;
                         if (s0[j] != '0')
                                 break;
                         --i;
                         }
                 e += nd - i;
                 nd = i;
                 if (nd0 > nd)
                         nd0 = nd;
                 if (nd < 9) { /* must recompute y */
@@ skipping 24 matching lines @@
                 if (bbe >= 0)
                         bb2 += bbe;
                 else
                         bd2 -= bbe;
                 bs2 = bb2;
 #ifdef Honor_FLT_ROUNDS
                 if (bc.rounding != 1)
                         bs2++;
 #endif
 #ifdef Avoid_Underflow
+                Lsb = LSB;
+                Lsb1 = 0;
                 j = bbe - bc.scale;
                 i = j + bbbits - 1;     /* logb(rv) */
-                if (i < Emin)   /* denormal */
-                        j += P - Emin;
-                else
-                        j = P + 1 - bbbits;
+                j = P + 1 - bbbits;
+                if (i < Emin) { /* denormal */
+                        i = Emin - i;
+                        j -= i;
+                        if (i < 32)
+                                Lsb <<= i;
+                        else if (i < 52)
+                                Lsb1 = Lsb << (i-32);
+                        else
+                                Lsb1 = Exp_mask;
+                        }
 #else /*Avoid_Underflow*/
 #ifdef Sudden_Underflow
 #ifdef IBM
                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
 #else
                 j = P + 1 - bbbits;
 #endif
 #else /*Sudden_Underflow*/
                 j = bbe;
                 i = j + bbbits - 1;     /* logb(rv) */
@@ skipping 27 matching lines @@
                 if (bd5 > 0)
                         bd = pow5mult(bd, bd5);
                 if (bd2 > 0)
                         bd = lshift(bd, bd2);
                 if (bs2 > 0)
                         bs = lshift(bs, bs2);
                 delta = diff(bb, bd);
                 bc.dsign = delta->sign;
                 delta->sign = 0;
                 i = cmp(delta, bs);
-#ifndef NO_STRTOD_BIGCOMP
+#ifndef NO_STRTOD_BIGCOMP /*{*/
                 if (bc.nd > nd && i <= 0) {
-                        if (bc.dsign)
-                                break;  /* Must use bigcomp(). */
+                        if (bc.dsign) {
+                                /* Must use bigcomp(). */
+                                req_bigcomp = 1;
+                                break;
+                                }
 #ifdef Honor_FLT_ROUNDS
                         if (bc.rounding != 1) {
-                                if (i < 0)
+                                if (i < 0) {
+                                        req_bigcomp = 1;
                                         break;
+                                        }
                                 }
                         else
 #endif
-                                {
-                                bc.nd = nd;
                                 i = -1; /* Discarded digits make delta smaller. */
-                                }
                         }
-#endif
-#ifdef Honor_FLT_ROUNDS
+#endif /*}*/
+#ifdef Honor_FLT_ROUNDS /*{*/
                 if (bc.rounding != 1) {
                         if (i < 0) {
                                 /* Error is less than an ulp */
                                 if (!delta->x[0] && delta->wds <= 1) {
                                         /* exact */
 #ifdef SET_INEXACT
                                         bc.inexact = 0;
 #endif
                                         break;
                                         }
@@ skipping 13 matching lines @@
 #else
                                                 if (y)
 #endif
                                                   {
                                                   delta = lshift(delta,Log2P);
                                                   if (cmp(delta, bs) <= 0)
                                                         adj.d = -0.5;
                                                   }
                                                 }
  apply_adj:
-#ifdef Avoid_Underflow
+#ifdef Avoid_Underflow /*{*/
                                         if (bc.scale && (y = word0(&rv) & Exp_mask)
                                                 <= 2*P*Exp_msk1)
                                           word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
                                         if ((word0(&rv) & Exp_mask) <=
                                                         P*Exp_msk1) {
                                                 word0(&rv) += P*Exp_msk1;
                                                 dval(&rv) += adj.d*ulp(dval(&rv));
                                                 word0(&rv) -= P*Exp_msk1;
                                                 }
                                         else
 #endif /*Sudden_Underflow*/
-#endif /*Avoid_Underflow*/
+#endif /*Avoid_Underflow}*/
                                         dval(&rv) += adj.d*ulp(&rv);
                                         }
                                 break;
                                 }
                         adj.d = ratio(delta, bs);
                         if (adj.d < 1.)
                                 adj.d = 1.;
                         if (adj.d <= 0x7ffffffe) {
                                 /* adj = rounding ? ceil(adj) : floor(adj); */
                                 y = adj.d;
                                 if (y != adj.d) {
                                         if (!((bc.rounding>>1) ^ bc.dsign))
                                                 y++;
                                         adj.d = y;
                                         }
                                 }
-#ifdef Avoid_Underflow
+#ifdef Avoid_Underflow /*{*/
                         if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
                                 word0(&adj) += (2*P+1)*Exp_msk1 - y;
 #else
 #ifdef Sudden_Underflow
                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
                                 word0(&rv) += P*Exp_msk1;
                                 adj.d *= ulp(dval(&rv));
                                 if (bc.dsign)
                                         dval(&rv) += adj.d;
                                 else
                                         dval(&rv) -= adj.d;
                                 word0(&rv) -= P*Exp_msk1;
                                 goto cont;
                                 }
 #endif /*Sudden_Underflow*/
-#endif /*Avoid_Underflow*/
+#endif /*Avoid_Underflow}*/
                         adj.d *= ulp(&rv);
                         if (bc.dsign) {
                                 if (word0(&rv) == Big0 && word1(&rv) == Big1)
                                         goto ovfl;
                                 dval(&rv) += adj.d;
                                 }
                         else
                                 dval(&rv) -= adj.d;
                         goto cont;
                         }
-#endif /*Honor_FLT_ROUNDS*/
+#endif /*}Honor_FLT_ROUNDS*/
 
                 if (i < 0) {
                         /* Error is less than half an ulp -- check for
                          * special case of mantissa a power of two.
                          */
                         if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
-#ifdef IEEE_Arith
+#ifdef IEEE_Arith /*{*/
 #ifdef Avoid_Underflow
                          || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
 #else
                          || (word0(&rv) & Exp_mask) <= Exp_msk1
 #endif
-#endif
+#endif /*}*/
                                 ) {
 #ifdef SET_INEXACT
                                 if (!delta->x[0] && delta->wds <= 1)
                                         bc.inexact = 0;
 #endif
                                 break;
                                 }
                         if (!delta->x[0] && delta->wds <= 1) {
                                 /* exact result */
 #ifdef SET_INEXACT
@@ skipping 10 matching lines @@
                         /* exactly half-way between */
                         if (bc.dsign) {
                                 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
                                  &&  word1(&rv) == (
 #ifdef Avoid_Underflow
                         (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
 #endif
                                                    0xffffffff)) {
                                         /*boundary case -- increment exponent*/
+                                        if (word0(&rv) == Big0 && word1(&rv) == Big1)
+                                                goto ovfl;
                                         word0(&rv) = (word0(&rv) & Exp_mask)
                                                 + Exp_msk1
 #ifdef IBM
                                                 | Exp_msk1 >> 4
 #endif
                                                 ;
                                         word1(&rv) = 0;
 #ifdef Avoid_Underflow
                                         bc.dsign = 0;
 #endif
@@ skipping 40 matching lines @@
                                                 }
                                         }
 #endif /*Avoid_Underflow*/
                                 L = (word0(&rv) & Exp_mask) - Exp_msk1;
 #endif /*Sudden_Underflow}}*/
                                 word0(&rv) = L | Bndry_mask1;
                                 word1(&rv) = 0xffffffff;
 #ifdef IBM
                                 goto cont;
 #else
+#ifndef NO_STRTOD_BIGCOMP
+                                if (bc.nd > nd)
+                                        goto cont;
+#endif
                                 break;
 #endif
                                 }
 #ifndef ROUND_BIASED
+#ifdef Avoid_Underflow
+                        if (Lsb1) {
+                                if (!(word0(&rv) & Lsb1))
+                                        break;
+                                }
+                        else if (!(word1(&rv) & Lsb))
+                                break;
+#else
                         if (!(word1(&rv) & LSB))
                                 break;
 #endif
+#endif
                         if (bc.dsign)
+#ifdef Avoid_Underflow
+                                dval(&rv) += sulp(&rv, &bc);
+#else
                                 dval(&rv) += ulp(&rv);
+#endif
 #ifndef ROUND_BIASED
                         else {
+#ifdef Avoid_Underflow
+                                dval(&rv) -= sulp(&rv, &bc);
+#else
                                 dval(&rv) -= ulp(&rv);
+#endif
 #ifndef Sudden_Underflow
                                 if (!dval(&rv)) {
                                         if (bc.nd >nd) {
                                                 bc.uflchk = 1;
                                                 break;
                                                 }
                                         goto undfl;
                                         }
 #endif
                                 }
@@ skipping 72 matching lines @@
                         if (bc.scale && y <= 2*P*Exp_msk1) {
                                 if (aadj <= 0x7fffffff) {
                                         if ((z = (ULong)aadj) <= 0)
                                                 z = 1;
                                         aadj = z;
                                         aadj1 = bc.dsign ? aadj : -aadj;
                                         }
                                 dval(&aadj2) = aadj1;
                                 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
                                 aadj1 = dval(&aadj2);
+                                adj.d = aadj1 * ulp(&rv);
+                                dval(&rv) += adj.d;
+                                if (rv.d == 0.)
+#ifdef NO_STRTOD_BIGCOMP
+                                        goto undfl;
+#else
+                                        {
+                                        req_bigcomp = 1;
+                                        break;
+                                        }
+#endif
                                 }
-                        adj.d = aadj1 * ulp(&rv);
-                        dval(&rv) += adj.d;
+                        else {
+                                adj.d = aadj1 * ulp(&rv);
+                                dval(&rv) += adj.d;
+                                }
 #else
 #ifdef Sudden_Underflow
                         if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
                                 dval(&rv0) = dval(&rv);
                                 word0(&rv) += P*Exp_msk1;
                                 adj.d = aadj1 * ulp(&rv);
                                 dval(&rv) += adj.d;
 #ifdef IBM
                                 if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
 #else
@@ skipping 62 matching lines @@
                 Bfree(bd);
                 Bfree(bs);
                 Bfree(delta);
                 }
         Bfree(bb);
         Bfree(bd);
         Bfree(bs);
         Bfree(bd0);
         Bfree(delta);
 #ifndef NO_STRTOD_BIGCOMP
-        if (bc.nd > nd)
+        if (req_bigcomp) {
+                bd0 = 0;
+                bc.e0 += nz1;
                 bigcomp(&rv, s0, &bc);
+                y = word0(&rv) & Exp_mask;
+                if (y == Exp_mask)
+                        goto ovfl;
+                if (y == 0 && rv.d == 0.)
+                        goto undfl;
+                }
 #endif
 #ifdef SET_INEXACT
         if (bc.inexact) {
                 if (!oldinexact) {
                         word0(&rv0) = Exp_1 + (70 << Exp_shift);
                         word1(&rv0) = 0;
                         dval(&rv0) += 1.;
                         }
                 }
         else if (!oldinexact)
@@ skipping 52 matching lines @@
 #ifndef MULTIPLE_THREADS
         dtoa_result =
 #endif
                 (char *)(r+1);
         }
 
  static char *
 #ifdef KR_headers
 nrv_alloc(s, rve, n) char *s, **rve; int n;
 #else
-nrv_alloc(CONST char *s, char **rve, int n)
+nrv_alloc(const char *s, char **rve, int n)
 #endif
 {
         char *rv, *t;
 
         t = rv = rv_alloc(n);
         for(*t = *s++; *t; *t = *s++) t++;
         if (rve)
                 *rve = t;
         return rv;
         }
@@ skipping 91 matching lines @@
                 6-9 ==> Debugging modes similar to mode - 4:  don't try
                         fast floating-point estimate (if applicable).
 
                 Values of mode other than 0-9 are treated as mode 0.
 
                 Sufficient space is allocated to the return value
                 to hold the suppressed trailing zeros.
         */
 
         int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
-                j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
+                j, j1 = 0, k, k0, k_check, leftright, m2, m5, s2, s5,
                 spec_case, try_quick;
         Long L;
 #ifndef Sudden_Underflow
         int denorm;
         ULong x;
 #endif
         Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
         U d2, eps, u;
         double ds;
         char *s, *s0;
+#ifndef No_leftright
+#ifdef IEEE_Arith
+        U eps1;
+#endif
+#endif
 #ifdef SET_INEXACT
         int inexact, oldinexact;
 #endif
 #ifdef Honor_FLT_ROUNDS /*{*/
         int Rounding;
 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
         Rounding = Flt_Rounds;
 #else /*}{*/
         Rounding = 1;
         switch(fegetround()) {
@@ skipping 245 matching lines @@
                         if (dval(&u) < -dval(&eps))
                                 goto no_digits;
                         goto fast_failed;
                         }
 #ifndef No_leftright
                 if (leftright) {
                         /* Use Steele & White method of only
                          * generating digits needed.
                          */
                         dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
+#ifdef IEEE_Arith
+                        if (k0 < 0 && j1 >= 307) {
+                                eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
+                                word0(&eps1) -= Exp_msk1 * (Bias+P-1);
+                                dval(&eps1) *= tens[j1 & 0xf];
+                                for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
+                                        if (j & 1)
+                                                dval(&eps1) *= bigtens[i];
+                                if (eps.d < eps1.d)
+                                        eps.d = eps1.d;
+                                }
+#endif
                         for(i = 0;;) {
-                                L = (long)dval(&u);
+                                L = dval(&u);
                                 dval(&u) -= L;
-                                *s++ = '0' + (char)L;
+                                *s++ = '0' + (int)L;
+                                if (1. - dval(&u) < dval(&eps))
+                                        goto bump_up;
                                 if (dval(&u) < dval(&eps))
                                         goto ret1;
-                                if (1. - dval(&u) < dval(&eps))
-                                        goto bump_up;
                                 if (++i >= ilim)
                                         break;
                                 dval(&eps) *= 10.;
                                 dval(&u) *= 10.;
                                 }
                         }
                 else {
 #endif
                         /* Generate ilim digits, then fix them up. */
                         dval(&eps) *= tens[ilim-1];
@@ skipping 27 matching lines @@
 
         if (be >= 0 && k <= Int_max) {
                 /* Yes. */
                 ds = tens[k];
                 if (ndigits < 0 && ilim <= 0) {
                         S = mhi = 0;
                         if (ilim < 0 || dval(&u) <= 5*ds)
                                 goto no_digits;
                         goto one_digit;
                         }
-                for(i = 1; i <= k + 1; i++, dval(&u) *= 10.) {
+                for(i = 1;; i++, dval(&u) *= 10.) {
                         L = (Long)(dval(&u) / ds);
                         dval(&u) -= L*ds;
 #ifdef Check_FLT_ROUNDS
                         /* If FLT_ROUNDS == 2, L will usually be high by 1 */
                         if (dval(&u) < 0) {
                                 L--;
                                 dval(&u) += ds;
                                 }
 #endif
                         *s++ = '0' + (char)L;
                         if (!dval(&u)) {
 #ifdef SET_INEXACT
                                 inexact = 0;
 #endif
                                 break;
                                 }
                         if (i == ilim) {
 #ifdef Honor_FLT_ROUNDS
                                 if (mode > 1)
                                 switch(Rounding) {
                                   case 0: goto ret1;
                                   case 2: goto bump_up;
                                   }
 #endif
                                 dval(&u) += dval(&u);
-                                if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
+#ifdef ROUND_BIASED
+                                if (dval(&u) >= ds)
+#else
+                                if (dval(&u) > ds || (dval(&u) == ds && L & 1))
+#endif
+                                        {
  bump_up:
                                         while(*--s == '9')
                                                 if (s == s0) {
                                                         k++;
                                                         *s = '0';
                                                         break;
                                                         }
                                         ++*s++;
                                         }
                                 break;
@@ skipping 64 matching lines @@
                         }
                 }
 
         /* Arrange for convenient computation of quotients:
          * shift left if necessary so divisor has 4 leading 0 bits.
          *
          * Perhaps we should just compute leading 28 bits of S once
          * and for all and pass them and a shift to quorem, so it
          * can do shifts and ors to compute the numerator for q.
          */
-#ifdef Pack_32
-        i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f;
-        if (i)
-                i = 32 - i;
-#define iInc 28
-#else
-        if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
-                i = 16 - i;
-#define iInc 12
-#endif
         i = dshift(S, s2);
         b2 += i;
         m2 += i;
         s2 += i;
         if (b2 > 0)
                 b = lshift(b, b2);
         if (s2 > 0)
                 S = lshift(S, s2);
         if (k_check) {
                 if (cmp(b,S) < 0) {
@@ skipping 72 matching lines @@
 #ifdef Honor_FLT_ROUNDS
                                 if (mode > 1)
                                  switch(Rounding) {
                                   case 0: goto accept_dig;
                                   case 2: goto keep_dig;
                                   }
 #endif /*Honor_FLT_ROUNDS*/
                                 if (j1 > 0) {
                                         b = lshift(b, 1);
                                         j1 = cmp(b, S);
+#ifdef ROUND_BIASED
+                                        if (j1 >= 0 /*)*/
+#else
                                         if ((j1 > 0 || (j1 == 0 && dig & 1))
+#endif
                                         && dig++ == '9')
                                                 goto round_9_up;
                                         }
  accept_dig:
                                 *s++ = (char)dig;
                                 goto ret;
                                 }
                         if (j1 > 0) {
 #ifdef Honor_FLT_ROUNDS
                                 if (!Rounding)
@@ skipping 40 matching lines @@
         /* Round off last digit */
 
 #ifdef Honor_FLT_ROUNDS
         switch(Rounding) {
           case 0: goto trimzeros;
           case 2: goto roundoff;
           }
 #endif
         b = lshift(b, 1);
         j = cmp(b, S);
-        if (j > 0 || (j == 0 && dig & 1)) {
+#ifdef ROUND_BIASED
+        if (j >= 0)
+#else
+        if (j > 0 || (j == 0 && dig & 1))
+#endif
+                {
  roundoff:
                 while(*--s == '9')
                         if (s == s0) {
                                 k++;
                                 *s++ = '1';
                                 goto ret;
                                 }
                 ++*s++;
                 }
         else {
@@ skipping 24 matching lines @@
 #endif
         Bfree(b);
         *s = 0;
         *decpt = k + 1;
         if (rve)
                 *rve = s;
         return s0;
         }
 
 }  // namespace dmg_fp