PR 78534, 83704 Handle large formatted I/O
[gcc.git] / libgfortran / io / write_float.def
index 776e59119931f07943af8787c1dd751bd60a7359..177a568e041d66bf1e302170a16b3d65da6611e6 100644 (file)
@@ -1,4 +1,4 @@
-/* Copyright (C) 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
+/* Copyright (C) 2007-2018 Free Software Foundation, Inc.
    Contributed by Andy Vaught
    Write float code factoring to this file by Jerry DeLisle   
    F2003 I/O support contributed by Jerry DeLisle
@@ -59,38 +59,84 @@ calculate_sign (st_parameter_dt *dtp, int negative_flag)
 }
 
 
-/* Output a real number according to its format which is FMT_G free.  */
+/* Determine the precision except for EN format. For G format,
+   determines an upper bound to be used for sizing the buffer. */
+
+static int
+determine_precision (st_parameter_dt * dtp, const fnode * f, int len)
+{
+  int precision = f->u.real.d;
+
+  switch (f->format)
+    {
+    case FMT_F:
+    case FMT_G:
+      precision += dtp->u.p.scale_factor;
+      break;
+    case FMT_ES:
+      /* Scale factor has no effect on output.  */
+      break;
+    case FMT_E:
+    case FMT_D:
+      /* See F2008 10.7.2.3.3.6 */
+      if (dtp->u.p.scale_factor <= 0)
+       precision += dtp->u.p.scale_factor - 1;
+      break;
+    default:
+      return -1;
+    }
+
+  /* If the scale factor has a large negative value, we must do our
+     own rounding? Use ROUND='NEAREST', which should be what snprintf
+     is using as well.  */
+  if (precision < 0 && 
+      (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED 
+       || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
+    dtp->u.p.current_unit->round_status = ROUND_NEAREST;
+
+  /* Add extra guard digits up to at least full precision when we do
+     our own rounding.  */
+  if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+      && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+    {
+      precision += 2 * len + 4;
+      if (precision < 0)
+       precision = 0;
+    }
+
+  return precision;
+}
+
+
+/* Build a real number according to its format which is FMT_G free.  */
 
 static void
-output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
-             int sign_bit, bool zero_flag, int ndigits, int edigits)
+build_float_string (st_parameter_dt *dtp, const fnode *f, char *buffer,
+                   size_t size, int nprinted, int precision, int sign_bit,
+                   bool zero_flag, int npad, char *result, size_t *len)
 {
-  char *out;
+  char *put;
   char *digits;
-  int e;
+  int e, w, d, p, i;
   char expchar, rchar;
   format_token ft;
-  int w;
-  int d;
   /* Number of digits before the decimal point.  */
   int nbefore;
   /* Number of zeros after the decimal point.  */
   int nzero;
   /* Number of digits after the decimal point.  */
   int nafter;
-  /* Number of zeros after the decimal point, whatever the precision.  */
-  int nzero_real;
   int leadzero;
   int nblanks;
-  int i;
+  int ndigits, edigits;
   sign_t sign;
 
   ft = f->format;
   w = f->u.real.w;
   d = f->u.real.d;
+  p = dtp->u.p.scale_factor;
 
   rchar = '5';
-  nzero_real = -1;
 
   /* We should always know the field width and precision.  */
   if (d < 0)
@@ -98,92 +144,124 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
 
   sign = calculate_sign (dtp, sign_bit);
   
-  /* The following code checks the given string has punctuation in the correct
-     places.  Uncomment if needed for debugging.
-     if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
-                   || buffer[ndigits + 2] != 'e'))
-       internal_error (&dtp->common, "printf is broken");  */
+  /* Calculate total number of digits.  */
+  if (ft == FMT_F)
+    ndigits = nprinted - 2;
+  else
+    ndigits = precision + 1;
 
   /* Read the exponent back in.  */
-  e = atoi (&buffer[ndigits + 3]) + 1;
+  if (ft != FMT_F)
+    e = atoi (&buffer[ndigits + 3]) + 1;
+  else
+    e = 0;
 
   /* Make sure zero comes out as 0.0e0.   */
   if (zero_flag)
-    {
-      e = 0;
-      if (compile_options.sign_zero == 1)
-       sign = calculate_sign (dtp, sign_bit);
-      else
-       sign = calculate_sign (dtp, 0);
-
-      /* Handle special cases.  */
-      if (w == 0)
-       w = d + 2;
-
-      /* For this one we choose to not output a decimal point.
-        F95 10.5.1.2.1  */
-      if (w == 1 && ft == FMT_F)
-       {
-         out = write_block (dtp, w);
-         if (out == NULL)
-           return;
-
-         if (unlikely (is_char4_unit (dtp)))
-           {
-             gfc_char4_t *out4 = (gfc_char4_t *) out;
-             *out4 = '0';
-             return;
-           }
-
-         *out = '0';
-         return;
-       }
-             
-    }
+    e = 0;
 
   /* Normalize the fractional component.  */
-  buffer[2] = buffer[1];
-  digits = &buffer[2];
+  if (ft != FMT_F)
+    {
+      buffer[2] = buffer[1];
+      digits = &buffer[2];
+    }
+  else
+    digits = &buffer[1];
 
   /* Figure out where to place the decimal point.  */
   switch (ft)
     {
     case FMT_F:
-      if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0)
-       {
-         memmove (digits + 1, digits, ndigits - 1);
-         digits[0] = '0';
-         e++;
+      nbefore = ndigits - precision;
+      if ((w > 0) && (nbefore > (int) size))
+        {
+         *len = w;
+         star_fill (result, w);
+         result[w] = '\0';
+         return;
        }
-
-      nbefore = e + dtp->u.p.scale_factor;
-      if (nbefore < 0)
+      /* Make sure the decimal point is a '.'; depending on the
+        locale, this might not be the case otherwise.  */
+      digits[nbefore] = '.';
+      if (p != 0)
        {
-         nzero = -nbefore;
-          nzero_real = nzero;
-         if (nzero > d)
-           nzero = d;
-         nafter = d - nzero;
-         nbefore = 0;
+         if (p > 0)
+           {
+             memmove (digits + nbefore, digits + nbefore + 1, p);
+             digits[nbefore + p] = '.';
+             nbefore += p;
+             nafter = d;
+             nzero = 0;
+           }
+         else /* p < 0  */
+           {
+             if (nbefore + p >= 0)
+               {
+                 nzero = 0;
+                 memmove (digits + nbefore + p + 1, digits + nbefore + p, -p);
+                 nbefore += p;
+                 digits[nbefore] = '.';
+                 nafter = d;
+               }
+             else
+               {
+                 nzero = -(nbefore + p);
+                 memmove (digits + 1, digits, nbefore);
+                 nafter = d - nzero;
+                 if (nafter == 0 && d > 0)
+                   {
+                     /* This is needed to get the correct rounding. */
+                     memmove (digits + 1, digits, ndigits - 1);
+                     digits[1] = '0';
+                     nafter = 1;
+                     nzero = d - 1;
+                   }
+                 else if (nafter < 0)
+                   {
+                     /* Reset digits to 0 in order to get correct rounding
+                        towards infinity. */
+                     for (i = 0; i < ndigits; i++)
+                       digits[i] = '0';
+                     digits[ndigits - 1] = '1';
+                     nafter = d;
+                     nzero = 0;
+                   }
+                 nbefore = 0;
+               }
+           }
        }
       else
        {
          nzero = 0;
          nafter = d;
        }
+
+      while (digits[0] == '0' && nbefore > 0)
+       {
+         digits++;
+         nbefore--;
+         ndigits--;
+       }
+
       expchar = 0;
+      /* If we need to do rounding ourselves, get rid of the dot by
+        moving the fractional part.  */
+      if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+         && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+       memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore);
       break;
 
     case FMT_E:
     case FMT_D:
       i = dtp->u.p.scale_factor;
-      if (d <= 0 && i == 0)
+      if (d <= 0 && p == 0)
        {
          generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
                          "greater than zero in format specifier 'E' or 'D'");
          return;
        }
-      if (i <= -d || i >= d + 2)
+      if (p <= -d || p >= d + 2)
        {
          generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
                          "out of range in format specifier 'E' or 'D'");
@@ -191,20 +269,20 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
        }
 
       if (!zero_flag)
-       e -= i;
-      if (i < 0)
+       e -= p;
+      if (p < 0)
        {
          nbefore = 0;
-         nzero = -i;
-         nafter = d + i;
+         nzero = -p;
+         nafter = d + p;
        }
-      else if (i > 0)
+      else if (p > 0)
        {
-         nbefore = i;
+         nbefore = p;
          nzero = 0;
-         nafter = (d - i) + 1;
+         nafter = (d - p) + 1;
        }
-      else /* i == 0 */
+      else /* p == 0 */
        {
          nbefore = 0;
          nzero = 0;
@@ -251,22 +329,29 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
       internal_error (&dtp->common, "Unexpected format token");
     }
 
-  /* Round the value.  The value being rounded is an unsigned magnitude.
-     The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
+  if (zero_flag)
+    goto skip;
+
+  /* Round the value.  The value being rounded is an unsigned magnitude.  */
   switch (dtp->u.p.current_unit->round_status)
     {
+      /* For processor defined and unspecified rounding we use
+        snprintf to print the exact number of digits needed, and thus
+        let snprintf handle the rounding.  On system claiming support
+        for IEEE 754, this ought to be round to nearest, ties to
+        even, corresponding to the Fortran ROUND='NEAREST'.  */
+      case ROUND_PROCDEFINED: 
+      case ROUND_UNSPECIFIED:
       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
        goto skip;
       case ROUND_UP:
        if (sign_bit)
          goto skip;
-       rchar = '0';
-       break;
+       goto updown;
       case ROUND_DOWN:
        if (!sign_bit)
          goto skip;
-       rchar = '0';
-       break;
+       goto updown;
       case ROUND_NEAREST:
        /* Round compatible unless there is a tie. A tie is a 5 with
           all trailing zero's.  */
@@ -278,7 +363,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
                if (digits[i] != '0')
                  goto do_rnd;
              }
-           /* It is a  tie so round to even.  */
+           /* It is a tie so round to even.  */
            switch (digits[nafter + nbefore - 1])
              {
                case '1':
@@ -293,32 +378,44 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
                  goto skip;
              }
          }
-        /* Fall through.  */ 
-      case ROUND_PROCDEFINED:
-      case ROUND_UNSPECIFIED:
+       /* Fall through.  */
+       /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
       case ROUND_COMPATIBLE:
        rchar = '5';
-       /* Just fall through and do the actual rounding.  */
+       goto do_rnd;
+    }
+
+  updown:
+
+  rchar = '0';
+  if (ft != FMT_F && w > 0 && d == 0 && p == 0)
+    nbefore = 1;
+  /* Scan for trailing zeros to see if we really need to round it.  */
+  for(i = nbefore + nafter; i < ndigits; i++)
+    {
+      if (digits[i] != '0')
+       goto do_rnd;
     }
+  goto skip;
     
   do_rnd:
  
   if (nbefore + nafter == 0)
+    /* Handle the case Fw.0 and value < 1.0 */
     {
       ndigits = 0;
-      if (nzero_real == d && digits[0] >= rchar)
+      if (digits[0] >= rchar)
        {
          /* We rounded to zero but shouldn't have */
-         nzero--;
-         nafter = 1;
+         nbefore = 1;
+         digits--;
          digits[0] = '1';
          ndigits = 1;
        }
     }
   else if (nbefore + nafter < ndigits)
     {
-      ndigits = nbefore + nafter;
-      i = ndigits;
+      i = ndigits = nbefore + nafter;
       if (digits[i] >= rchar)
        {
          /* Propagate the carry.  */
@@ -368,7 +465,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
   skip:
 
   /* Calculate the format of the exponent field.  */
-  if (expchar)
+  if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
     {
       edigits = 1;
       for (i = abs (e); i >= 10; i /= 10)
@@ -398,15 +495,24 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
   else
     edigits = 0;
 
-  /* Zero values always output as positive, even if the value was negative
-     before rounding.  */
-  for (i = 0; i < ndigits; i++)
+  /* Scan the digits string and count the number of zeros.  If we make it
+     all the way through the loop, we know the value is zero after the
+     rounding completed above.  */
+  int hasdot = 0;
+  for (i = 0; i < ndigits + hasdot; i++)
     {
-      if (digits[i] != '0')
+      if (digits[i] == '.')
+       hasdot = 1;
+      else if (digits[i] != '0')
        break;
     }
-  if (i == ndigits)
+
+  /* To format properly, we need to know if the rounded result is zero and if
+     so, we set the zero_flag which may have been already set for
+     actual zero.  */
+  if (i == ndigits + hasdot)
     {
+      zero_flag = true;
       /* The output is zero, so set the sign according to the sign bit unless
         -fno-sign-zero was specified.  */
       if (compile_options.sign_zero == 1)
@@ -415,39 +521,24 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
        sign = calculate_sign (dtp, 0);
     }
 
-  /* Pick a field size if none was specified.  */
+  /* Pick a field size if none was specified, taking into account small
+     values that may have been rounded to zero.  */
   if (w <= 0)
-    w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
-  
-  /* Work out how much padding is needed.  */
-  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
-  if (sign != S_NONE)
-    nblanks--;
-
-  if (dtp->u.p.g0_no_blanks)
     {
-      w -= nblanks;
-      nblanks = 0;
-    }
-
-  /* Create the ouput buffer.  */
-  out = write_block (dtp, w);
-  if (out == NULL)
-    return;
-
-  /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1)
-    {
-      if (unlikely (is_char4_unit (dtp)))
+      if (zero_flag)
+       w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+      else
        {
-         gfc_char4_t *out4 = (gfc_char4_t *) out;
-         memset4 (out4, '*', w);
-         return;
+         w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+         w = w == 1 ? 2 : w;
        }
-      star_fill (out, w);
-      return;
     }
 
+  /* Work out how much padding is needed.  */
+  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
+  if (sign != S_NONE)
+    nblanks--;
+
   /* See if we have space for a zero before the decimal point.  */
   if (nbefore == 0 && nblanks > 0)
     {
@@ -457,156 +548,77 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
   else
     leadzero = 0;
 
-  /* For internal character(kind=4) units, we duplicate the code used for
-     regular output slightly modified.  This needs to be maintained
-     consistent with the regular code that follows this block.  */
-  if (unlikely (is_char4_unit (dtp)))
+  if (dtp->u.p.g0_no_blanks)
     {
-      gfc_char4_t *out4 = (gfc_char4_t *) out;
-      /* Pad to full field width.  */
-
-      if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
-       {
-         memset4 (out4, ' ', nblanks);
-         out4 += nblanks;
-       }
-
-      /* Output the initial sign (if any).  */
-      if (sign == S_PLUS)
-       *(out4++) = '+';
-      else if (sign == S_MINUS)
-       *(out4++) = '-';
-
-      /* Output an optional leading zero.  */
-      if (leadzero)
-       *(out4++) = '0';
-
-      /* Output the part before the decimal point, padding with zeros.  */
-      if (nbefore > 0)
-       {
-         if (nbefore > ndigits)
-           {
-             i = ndigits;
-             memcpy4 (out4, digits, i);
-             ndigits = 0;
-             while (i < nbefore)
-               out4[i++] = '0';
-           }
-         else
-           {
-             i = nbefore;
-             memcpy4 (out4, digits, i);
-             ndigits -= i;
-           }
-
-         digits += i;
-         out4 += nbefore;
-       }
-
-      /* Output the decimal point.  */
-      *(out4++) = dtp->u.p.current_unit->decimal_status
-                   == DECIMAL_POINT ? '.' : ',';
-
-      /* Output leading zeros after the decimal point.  */
-      if (nzero > 0)
-       {
-         for (i = 0; i < nzero; i++)
-           *(out4++) = '0';
-       }
-
-      /* Output digits after the decimal point, padding with zeros.  */
-      if (nafter > 0)
-       {
-         if (nafter > ndigits)
-           i = ndigits;
-         else
-           i = nafter;
-
-         memcpy4 (out4, digits, i);
-         while (i < nafter)
-           out4[i++] = '0';
+      w -= nblanks;
+      nblanks = 0;
+    }
 
-         digits += i;
-         ndigits -= i;
-         out4 += nafter;
-       }
+  /* Create the final float string.  */
+  *len = w + npad;
+  put = result;
 
-      /* Output the exponent.  */
-      if (expchar)
-       {
-         if (expchar != ' ')
-           {
-             *(out4++) = expchar;
-             edigits--;
-           }
-#if HAVE_SNPRINTF
-         snprintf (buffer, size, "%+0*d", edigits, e);
-#else
-         sprintf (buffer, "%+0*d", edigits, e);
-#endif
-         memcpy4 (out4, buffer, edigits);
-       }
-
-      if (dtp->u.p.no_leading_blank)
-       {
-         out4 += edigits;
-         memset4 (out4, ' ' , nblanks);
-         dtp->u.p.no_leading_blank = 0;
-       }
+  /* Check the value fits in the specified field width.  */
+  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
+    {
+      star_fill (put, *len);
       return;
-    } /* End of character(kind=4) internal unit code.  */
+    }
 
   /* Pad to full field width.  */
-
   if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
     {
-      memset (out, ' ', nblanks);
-      out += nblanks;
+      memset (put, ' ', nblanks);
+      put += nblanks;
     }
 
-  /* Output the initial sign (if any).  */
+  /* Set the initial sign (if any).  */
   if (sign == S_PLUS)
-    *(out++) = '+';
+    *(put++) = '+';
   else if (sign == S_MINUS)
-    *(out++) = '-';
+    *(put++) = '-';
 
-  /* Output an optional leading zero.  */
+  /* Set an optional leading zero.  */
   if (leadzero)
-    *(out++) = '0';
+    *(put++) = '0';
 
-  /* Output the part before the decimal point, padding with zeros.  */
+  /* Set the part before the decimal point, padding with zeros.  */
   if (nbefore > 0)
     {
       if (nbefore > ndigits)
        {
          i = ndigits;
-         memcpy (out, digits, i);
+         memcpy (put, digits, i);
          ndigits = 0;
          while (i < nbefore)
-           out[i++] = '0';
+           put[i++] = '0';
        }
       else
        {
          i = nbefore;
-         memcpy (out, digits, i);
+         memcpy (put, digits, i);
          ndigits -= i;
        }
 
       digits += i;
-      out += nbefore;
+      put += nbefore;
     }
 
-  /* Output the decimal point.  */
-  *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
+  /* Set the decimal point.  */
+  *(put++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
+  if (ft == FMT_F
+         && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED 
+             || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
+    digits++;
 
-  /* Output leading zeros after the decimal point.  */
+  /* Set leading zeros after the decimal point.  */
   if (nzero > 0)
     {
       for (i = 0; i < nzero; i++)
-       *(out++) = '0';
+       *(put++) = '0';
     }
 
-  /* Output digits after the decimal point, padding with zeros.  */
+  /* Set digits after the decimal point, padding with zeros.  */
   if (nafter > 0)
     {
       if (nafter > ndigits)
@@ -614,82 +626,87 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
       else
        i = nafter;
 
-      memcpy (out, digits, i);
+      memcpy (put, digits, i);
       while (i < nafter)
-       out[i++] = '0';
+       put[i++] = '0';
 
       digits += i;
       ndigits -= i;
-      out += nafter;
+      put += nafter;
     }
 
-  /* Output the exponent.  */
-  if (expchar)
+  /* Set the exponent.  */
+  if (expchar && !(dtp->u.p.g0_no_blanks && e == 0))
     {
       if (expchar != ' ')
        {
-         *(out++) = expchar;
+         *(put++) = expchar;
          edigits--;
        }
-#if HAVE_SNPRINTF
       snprintf (buffer, size, "%+0*d", edigits, e);
-#else
-      sprintf (buffer, "%+0*d", edigits, e);
-#endif
-      memcpy (out, buffer, edigits);
+      memcpy (put, buffer, edigits);
+      put += edigits;
     }
 
   if (dtp->u.p.no_leading_blank)
     {
-      out += edigits;
-      memset( out , ' ' , nblanks );
+      memset (put , ' ' , nblanks);
       dtp->u.p.no_leading_blank = 0;
+      put += nblanks;
+    }
+
+  if (npad > 0 && !dtp->u.p.g0_no_blanks)
+    {
+      memset (put , ' ' , npad);
+      put += npad;
     }
 
-#undef STR
-#undef STR1
-#undef MIN_FIELD_WIDTH
+  /* NULL terminate the string.  */
+  *put = '\0';
+  
+  return;
 }
 
 
 /* Write "Infinite" or "Nan" as appropriate for the given format.  */
 
 static void
-write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
+build_infnan_string (st_parameter_dt *dtp, const fnode *f, int isnan_flag,
+                   int sign_bit, char *p, size_t *len)
 {
-  char * p, fin;
+  char fin;
   int nb = 0;
+  sign_t sign;
+  int mark;
 
   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
     {
+      sign = calculate_sign (dtp, sign_bit);
+      mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
+
       nb =  f->u.real.w;
-  
+      *len = nb;
+
       /* If the field width is zero, the processor must select a width 
         not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
      
-      if (nb == 0) nb = 4;
-      p = write_block (dtp, nb);
-      if (p == NULL)
-       return;
-      if (nb < 3)
+      if ((nb == 0) || dtp->u.p.g0_no_blanks)
        {
-         if (unlikely (is_char4_unit (dtp)))
-           {
-             gfc_char4_t *p4 = (gfc_char4_t *) p;
-             memset4 (p4, '*', nb);
-           }
+         if (isnan_flag)
+           nb = 3;
          else
-           memset (p, '*', nb);
-         return;
+           nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
+         *len = nb;
        }
 
-      if (unlikely (is_char4_unit (dtp)))
+      p[*len] = '\0';
+      if (nb < 3)
        {
-         gfc_char4_t *p4 = (gfc_char4_t *) p;
-         memset4 (p4, ' ', nb);
+         memset (p, '*', nb);
+         return;
        }
-      else
-       memset(p, ' ', nb);
+
+      memset(p, ' ', nb);
 
       if (!isnan_flag)
        {
@@ -699,13 +716,7 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
                 insufficient room to output '-Inf', so output asterisks */
              if (nb == 3)
                {
-                 if (unlikely (is_char4_unit (dtp)))
-                   {
-                     gfc_char4_t *p4 = (gfc_char4_t *) p;
-                     memset4 (p4, '*', nb);
-                   }
-                 else
-                   memset (p, '*', nb);
+                 memset (p, '*', nb);
                  return;
                }
              /* The negative sign is mandatory */
@@ -716,27 +727,7 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
               consistency */
            fin = '+';
            
-         if (unlikely (is_char4_unit (dtp)))
-           {
-             gfc_char4_t *p4 = (gfc_char4_t *) p;
-             if (nb > 8)
-               /* We have room, so output 'Infinity' */
-               memcpy4 (p4 + nb - 8, "Infinity", 8);
-             else
-               /* For the case of width equals 8, there is not enough room
-                  for the sign and 'Infinity' so we go with 'Inf' */
-               memcpy4 (p4 + nb - 3, "Inf", 3);
-
-             if (nb < 9 && nb > 3)
-               /* Put the sign in front of Inf */
-               p4[nb - 4] = (gfc_char4_t) fin;
-             else if (nb > 8)
-               /* Put the sign in front of Infinity */
-               p4[nb - 9] = (gfc_char4_t) fin;
-             return;
-           }
-
-         if (nb > 8)
+         if (nb > mark)
            /* We have room, so output 'Infinity' */
            memcpy(p + nb - 8, "Infinity", 8);
          else
@@ -744,22 +735,16 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
               for the sign and 'Infinity' so we go with 'Inf' */
            memcpy(p + nb - 3, "Inf", 3);
 
-         if (nb < 9 && nb > 3)
-           p[nb - 4] = fin;  /* Put the sign in front of Inf */
-         else if (nb > 8)
-           p[nb - 9] = fin;  /* Put the sign in front of Infinity */
-       }
-      else
-        {
-         if (unlikely (is_char4_unit (dtp)))
+         if (sign == S_PLUS || sign == S_MINUS)
            {
-             gfc_char4_t *p4 = (gfc_char4_t *) p;
-             memcpy4 (p4 + nb - 3, "NaN", 3);
+             if (nb < 9 && nb > 3)
+               p[nb - 4] = fin;  /* Put the sign in front of Inf */
+             else if (nb > 8)
+               p[nb - 9] = fin;  /* Put the sign in front of Infinity */
            }
-         else
-           memcpy(p + nb - 3, "NaN", 3);
        }
-      return;
+      else
+       memcpy(p + nb - 3, "NaN", 3);
     }
 }
 
@@ -767,7 +752,7 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
 /* Returns the value of 10**d.  */
 
 #define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
+static GFC_REAL_ ## x \
 calculate_exp_ ## x  (int d)\
 {\
   int i;\
@@ -791,149 +776,11 @@ CALCULATE_EXP(16)
 #endif
 #undef CALCULATE_EXP
 
-/* Generate corresponding I/O format for FMT_G and output.
-   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
-   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 
-   Data Magnitude                              Equivalent Conversion
-   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
-   m = 0                                       F(w-n).(d-1), n' '
-   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
-   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
-   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
-   ................                           ..........
-   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
-   m >= 10**d-0.5                              Ew.d[Ee]
+/* Define macros to build code for format_float.  */
 
-   notes: for Gw.d ,  n' ' means 4 blanks
-          for Gw.dEe, n' ' means e+2 blanks  */
-
-#define OUTPUT_FLOAT_FMT_G(x) \
-static void \
-output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
-                     GFC_REAL_ ## x m, char *buffer, size_t size, \
-                     int sign_bit, bool zero_flag, int ndigits, int edigits) \
-{ \
-  int e = f->u.real.e;\
-  int d = f->u.real.d;\
-  int w = f->u.real.w;\
-  fnode *newf;\
-  GFC_REAL_ ## x rexp_d;\
-  int low, high, mid;\
-  int ubound, lbound;\
-  char *p;\
-  int save_scale_factor, nb = 0;\
-\
-  save_scale_factor = dtp->u.p.scale_factor;\
-  newf = (fnode *) get_mem (sizeof (fnode));\
-\
-  rexp_d = calculate_exp_ ## x (-d);\
-  if ((m > 0.0 && m < 0.1 - 0.05 * rexp_d) || (rexp_d * (m + 0.5) >= 1.0) ||\
-      ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
-    { \
-      newf->format = FMT_E;\
-      newf->u.real.w = w;\
-      newf->u.real.d = d;\
-      newf->u.real.e = e;\
-      nb = 0;\
-      goto finish;\
-    }\
-\
-  mid = 0;\
-  low = 0;\
-  high = d + 1;\
-  lbound = 0;\
-  ubound = d + 1;\
-\
-  while (low <= high)\
-    { \
-      GFC_REAL_ ## x temp;\
-      mid = (low + high) / 2;\
-\
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
-\
-      if (m < temp)\
-        { \
-          ubound = mid;\
-          if (ubound == lbound + 1)\
-            break;\
-          high = mid - 1;\
-        }\
-      else if (m > temp)\
-        { \
-          lbound = mid;\
-          if (ubound == lbound + 1)\
-            { \
-              mid ++;\
-              break;\
-            }\
-          low = mid + 1;\
-        }\
-      else\
-       {\
-         mid++;\
-         break;\
-       }\
-    }\
-\
-  if (e < 0)\
-    nb = 4;\
-  else\
-    nb = e + 2;\
-\
-  newf->format = FMT_F;\
-  newf->u.real.w = f->u.real.w - nb;\
-\
-  if (m == 0.0)\
-    newf->u.real.d = d - 1;\
-  else\
-    newf->u.real.d = - (mid - d - 1);\
-\
-  dtp->u.p.scale_factor = 0;\
-\
- finish:\
-  output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
-               edigits);\
-  dtp->u.p.scale_factor = save_scale_factor;\
-\
-  free (newf);\
-\
-  if (nb > 0 && !dtp->u.p.g0_no_blanks)\
-    {\
-      p = write_block (dtp, nb);\
-      if (p == NULL)\
-       return;\
-      if (unlikely (is_char4_unit (dtp)))\
-       {\
-         gfc_char4_t *p4 = (gfc_char4_t *) p;\
-         memset4 (p4, ' ', nb);\
-       }\
-      else\
-       memset (p, ' ', nb);\
-    }\
-}\
-
-OUTPUT_FLOAT_FMT_G(4)
-
-OUTPUT_FLOAT_FMT_G(8)
-
-#ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
-#endif
-
-#undef OUTPUT_FLOAT_FMT_G
-
-
-/* Define a macro to build code for write_float.  */
-
-  /* Note: Before output_float is called, sprintf is used to print to buffer the
-     number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
-     (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
-     before the decimal point.
+  /* Note: Before output_float is called, snprintf is used to print to buffer the
+     number in the format +D.DDDDe+ddd. 
 
      #   The result will always contain a decimal point, even if no
         digits follow it
@@ -942,125 +789,301 @@ OUTPUT_FLOAT_FMT_G(16)
 
      +   A sign (+ or -) always be placed before a number
 
-     MIN_FIELD_WIDTH  minimum field width
-
-     *   (ndigits-1) is used as the precision
+     *   prec is used as the precision
 
      e format: [-]d.ddde±dd where there is one digit before the
        decimal-point character and the number of digits after it is
        equal to the precision. The exponent always contains at least two
        digits; if the value is zero, the exponent is 00.  */
 
-#ifdef HAVE_SNPRINTF
 
-#define DTOA \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-         "e", ndigits - 1, tmp);
+#define TOKENPASTE(x, y) TOKENPASTE2(x, y)
+#define TOKENPASTE2(x, y) x ## y
+
+#define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
+
+#define DTOA2(prec,val) \
+snprintf (buffer, size, "%+-#.*e", (prec), (val))
 
-#define DTOAL \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-         "Le", ndigits - 1, tmp);
+#define DTOA2L(prec,val) \
+snprintf (buffer, size, "%+-#.*Le", (prec), (val))
 
-#else
 
-#define DTOA \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-        "e", ndigits - 1, tmp);
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define DTOA2Q(prec,val) \
+quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val))
+#endif
+
+#define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
+
+/* For F format, we print to the buffer with f format.  */
+#define FDTOA2(prec,val) \
+snprintf (buffer, size, "%+-#.*f", (prec), (val))
 
-#define DTOAL \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-        "Le", ndigits - 1, tmp);
+#define FDTOA2L(prec,val) \
+snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
 
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define FDTOA2Q(prec,val) \
+quadmath_snprintf (buffer, size, "%+-#.*Qf", \
+                            (prec), (val))
 #endif
 
-#define WRITE_FLOAT(x,y)\
+
+/* EN format is tricky since the number of significant digits depends
+   on the magnitude.  Solve it by first printing a temporary value and
+   figure out the number of significant digits from the printed
+   exponent.  Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
+   10.0**e even when the final result will not be rounded to 10.0**e.
+   For these values the exponent returned by atoi has to be decremented
+   by one. The values y in the ranges
+       (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))  
+        (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
+         (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
+   are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
+   100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
+   represents d zeroes, by the lines 279 to 297. */
+#define EN_PREC(x,y)\
 {\
-       GFC_REAL_ ## x tmp;\
-       tmp = * (GFC_REAL_ ## x *)source;\
-       sign_bit = __builtin_signbit (tmp);\
-       if (!isfinite (tmp))\
-         { \
-           write_infnan (dtp, f, isnan (tmp), sign_bit);\
-           return;\
+    volatile GFC_REAL_ ## x tmp, one = 1.0;\
+    tmp = * (GFC_REAL_ ## x *)source;\
+    if (isfinite (tmp))\
+      {\
+       nprinted = DTOA(y,0,tmp);\
+       int e = atoi (&buffer[4]);\
+       if (buffer[1] == '1')\
+         {\
+           tmp = (calculate_exp_ ## x (-e)) * tmp;\
+           tmp = one - (tmp < 0 ? -tmp : tmp);\
+           if (tmp > 0)\
+             e = e - 1;\
          }\
-       tmp = sign_bit ? -tmp : tmp;\
-       zero_flag = (tmp == 0.0);\
-\
-       DTOA ## y\
-\
-       if (f->format != FMT_G)\
-         output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
-                       edigits);\
-       else \
-         output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
-                                   zero_flag, ndigits, edigits);\
+       nbefore = e%3;\
+       if (nbefore < 0)\
+         nbefore = 3 + nbefore;\
+      }\
+    else\
+      nprinted = -1;\
 }\
 
-/* Output a real number according to its format.  */
-
-static void
-write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
+static int
+determine_en_precision (st_parameter_dt *dtp, const fnode *f, 
+                       const char *source, int len)
 {
+  int nprinted;
+  char buffer[10];
+  const size_t size = 10;
+  int nbefore; /* digits before decimal point - 1.  */
+
+  switch (len)
+    {
+    case 4:
+      EN_PREC(4,)
+      break;
 
-#if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
-# define MIN_FIELD_WIDTH 46
-#else
-# define MIN_FIELD_WIDTH 31
+    case 8:
+      EN_PREC(8,)
+      break;
+
+#ifdef HAVE_GFC_REAL_10
+    case 10:
+      EN_PREC(10,L)
+      break;
+#endif
+#ifdef HAVE_GFC_REAL_16
+    case 16:
+# ifdef GFC_REAL_16_IS_FLOAT128
+      EN_PREC(16,Q)
+# else
+      EN_PREC(16,L)
+# endif
+      break;
 #endif
-#define STR(x) STR1(x)
-#define STR1(x) #x
+    default:
+      internal_error (NULL, "bad real kind");
+    }
 
-  /* This must be large enough to accurately hold any value.  */
-  char buffer[MIN_FIELD_WIDTH+1];
-  int sign_bit, ndigits, edigits;
-  bool zero_flag;
-  size_t size;
+  if (nprinted == -1)
+    return -1;
 
-  size = MIN_FIELD_WIDTH+1;
+  int prec = f->u.real.d + nbefore;
+  if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
+      && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
+    prec += 2 * len + 4;
+  return prec;
+}
+  
 
-  /* printf pads blanks for us on the exponent so we just need it big enough
-     to handle the largest number of exponent digits expected.  */
-  edigits=4;
+/* Generate corresponding I/O format. and output.
+   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
+   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
 
-  if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G 
-      || ((f->format == FMT_D || f->format == FMT_E)
-      && dtp->u.p.scale_factor != 0))
-    {
-      /* Always convert at full precision to avoid double rounding.  */
-      ndigits = MIN_FIELD_WIDTH - 4 - edigits;
-    }
-  else
-    {
-      /* The number of digits is known, so let printf do the rounding.  */
-      if (f->format == FMT_ES)
-       ndigits = f->u.real.d + 1;
-      else
-       ndigits = f->u.real.d;
-      if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
-       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
-    }
+   Data Magnitude                              Equivalent Conversion
+   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
+   m = 0                                       F(w-n).(d-1), n' '
+   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
+   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
+   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
+   ................                           ..........
+   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
+   m >= 10**d-0.5                              Ew.d[Ee]
 
-  switch (len)
+   notes: for Gw.d ,  n' ' means 4 blanks
+         for Gw.dEe, n' ' means e+2 blanks
+         for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
+         the asm volatile is required for 32-bit x86 platforms.  */
+#define FORMAT_FLOAT(x,y)\
+{\
+  int npad = 0;\
+  GFC_REAL_ ## x m;\
+  m = * (GFC_REAL_ ## x *)source;\
+  sign_bit = signbit (m);\
+  if (!isfinite (m))\
+    { \
+      build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\
+      return;\
+    }\
+  m = sign_bit ? -m : m;\
+  zero_flag = (m == 0.0);\
+  if (f->format == FMT_G)\
+    {\
+      int e = f->u.real.e;\
+      int d = f->u.real.d;\
+      int w = f->u.real.w;\
+      fnode newf;\
+      GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
+      int low, high, mid;\
+      int ubound, lbound;\
+      int save_scale_factor;\
+      volatile GFC_REAL_ ## x temp;\
+      save_scale_factor = dtp->u.p.scale_factor;\
+      switch (dtp->u.p.current_unit->round_status)\
+       {\
+         case ROUND_ZERO:\
+           r = sign_bit ? 1.0 : 0.0;\
+           break;\
+         case ROUND_UP:\
+           r = 1.0;\
+           break;\
+         case ROUND_DOWN:\
+           r = 0.0;\
+           break;\
+         default:\
+           break;\
+       }\
+      exp_d = calculate_exp_ ## x (d);\
+      r_sc = (1 - r / exp_d);\
+      temp = 0.1 * r_sc;\
+      if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
+         || ((m == 0.0) && !(compile_options.allow_std\
+                             & (GFC_STD_F2003 | GFC_STD_F2008)))\
+         ||  d == 0)\
+       { \
+         newf.format = FMT_E;\
+         newf.u.real.w = w;\
+         newf.u.real.d = d - comp_d;\
+         newf.u.real.e = e;\
+         npad = 0;\
+         precision = determine_precision (dtp, &newf, x);\
+         nprinted = DTOA(y,precision,m);\
+       }\
+      else \
+       {\
+         mid = 0;\
+         low = 0;\
+         high = d + 1;\
+         lbound = 0;\
+         ubound = d + 1;\
+         while (low <= high)\
+           {\
+             mid = (low + high) / 2;\
+             temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
+             if (m < temp)\
+               { \
+                 ubound = mid;\
+                 if (ubound == lbound + 1)\
+                   break;\
+                 high = mid - 1;\
+               }\
+             else if (m > temp)\
+               { \
+                 lbound = mid;\
+                 if (ubound == lbound + 1)\
+                   { \
+                     mid ++;\
+                     break;\
+                   }\
+                 low = mid + 1;\
+               }\
+             else\
+               {\
+                 mid++;\
+                 break;\
+               }\
+           }\
+         npad = e <= 0 ? 4 : e + 2;\
+         npad = npad >= w ? w - 1 : npad;\
+         npad = dtp->u.p.g0_no_blanks ? 0 : npad;\
+         newf.format = FMT_F;\
+         newf.u.real.w = w - npad;\
+         newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
+         dtp->u.p.scale_factor = 0;\
+         precision = determine_precision (dtp, &newf, x);\
+         nprinted = FDTOA(y,precision,m);\
+       }\
+      build_float_string (dtp, &newf, buffer, size, nprinted, precision,\
+                                  sign_bit, zero_flag, npad, result, res_len);\
+      dtp->u.p.scale_factor = save_scale_factor;\
+    }\
+  else\
+    {\
+      if (f->format == FMT_F)\
+       nprinted = FDTOA(y,precision,m);\
+      else\
+       nprinted = DTOA(y,precision,m);\
+      build_float_string (dtp, f, buffer, size, nprinted, precision,\
+                                  sign_bit, zero_flag, npad, result, res_len);\
+    }\
+}\
+
+/* Output a real number according to its format.  */
+
+
+static void
+get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source,
+                 int kind, int comp_d, char *buffer, int precision,
+                 size_t size, char *result, size_t *res_len)
+{
+  int sign_bit, nprinted;
+  bool zero_flag;
+
+  switch (kind)
     {
     case 4:
-      WRITE_FLOAT(4,)
+      FORMAT_FLOAT(4,)
       break;
 
     case 8:
-      WRITE_FLOAT(8,)
+      FORMAT_FLOAT(8,)
       break;
 
 #ifdef HAVE_GFC_REAL_10
     case 10:
-      WRITE_FLOAT(10,L)
+      FORMAT_FLOAT(10,L)
       break;
 #endif
 #ifdef HAVE_GFC_REAL_16
     case 16:
-      WRITE_FLOAT(16,L)
+# ifdef GFC_REAL_16_IS_FLOAT128
+      FORMAT_FLOAT(16,Q)
+# else
+      FORMAT_FLOAT(16,L)
+# endif
       break;
 #endif
     default:
       internal_error (NULL, "bad real kind");
     }
+  return;
 }