PR libfortran/64770 Segfault when trying to open existing file with status="new".
[gcc.git] / libgfortran / io / write_float.def
index 02e1b8b9b13fcd5cb3db872d5234a3942778b873..1bcd8159a22085f064db98dea32a2fcf82016c85 100644 (file)
@@ -1,4 +1,4 @@
-/* Copyright (C) 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
+/* Copyright (C) 2007-2015 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,83 @@ calculate_sign (st_parameter_dt *dtp, int negative_flag)
 }
 
 
+/* 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;
+}
+
+
 /* Output 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)
+static bool
+output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
+             int nprinted, int precision, int sign_bit, bool zero_flag)
 {
   char *out;
   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,113 +143,127 @@ 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)
+      nbefore = ndigits - precision;
+      /* Make sure the decimal point is a '.'; depending on the
+        locale, this might not be the case otherwise.  */
+      digits[nbefore] = '.';
+      if (p != 0)
        {
-         memmove (digits + 1, digits, ndigits - 1);
-         digits[0] = '0';
-         e++;
-       }
-
-      nbefore = e + dtp->u.p.scale_factor;
-      if (nbefore < 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 - p;
+             if (nafter < 0)
+               nafter = 0;
+             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);
+                 digits++;
+                 nafter = d + nbefore;
+                 nbefore = 0;
+               }
+             if (nzero > d)
+               nzero = d;
+           }
        }
       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;
+         return false;
        }
-      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'");
-         return;
+         return false;
        }
 
       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 +310,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 +344,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 +359,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.  */
@@ -398,15 +476,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,9 +502,18 @@ 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);
+    {
+      if (zero_flag)
+       w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
+      else
+       {
+         w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
+         w = w == 1 ? 2 : w;
+       }
+    }
   
   /* Work out how much padding is needed.  */
   nblanks = w - (nbefore + nzero + nafter + edigits + 1);
@@ -433,18 +529,19 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
   /* Create the ouput buffer.  */
   out = write_block (dtp, w);
   if (out == NULL)
-    return;
+    return false;
 
   /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1)
+  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
     {
       if (unlikely (is_char4_unit (dtp)))
        {
-         memset4 (out, 0, '*', w);
-         return;
+         gfc_char4_t *out4 = (gfc_char4_t *) out;
+         memset4 (out4, '*', w);
+         return false;
        }
       star_fill (out, w);
-      return;
+      return false;
     }
 
   /* See if we have space for a zero before the decimal point.  */
@@ -466,7 +563,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
 
       if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
        {
-         memset4 (out, 0, ' ', nblanks);
+         memset4 (out4, ' ', nblanks);
          out4 += nblanks;
        }
 
@@ -486,7 +583,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
          if (nbefore > ndigits)
            {
              i = ndigits;
-             memcpy4 (out4, 0, digits, i);
+             memcpy4 (out4, digits, i);
              ndigits = 0;
              while (i < nbefore)
                out4[i++] = '0';
@@ -494,7 +591,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
          else
            {
              i = nbefore;
-             memcpy4 (out4, 0, digits, i);
+             memcpy4 (out4, digits, i);
              ndigits -= i;
            }
 
@@ -505,6 +602,10 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
       /* Output the decimal point.  */
       *(out4++) = 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.  */
       if (nzero > 0)
@@ -521,7 +622,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
          else
            i = nafter;
 
-         memcpy4 (out4, 0, digits, i);
+         memcpy4 (out4, digits, i);
          while (i < nafter)
            out4[i++] = '0';
 
@@ -538,21 +639,17 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
              *(out4++) = expchar;
              edigits--;
            }
-#if HAVE_SNPRINTF
          snprintf (buffer, size, "%+0*d", edigits, e);
-#else
-         sprintf (buffer, "%+0*d", edigits, e);
-#endif
-         memcpy4 (out4, 0, buffer, edigits);
+         memcpy4 (out4, buffer, edigits);
        }
 
       if (dtp->u.p.no_leading_blank)
        {
          out4 += edigits;
-         memset4 (out4 , 0, ' ' , nblanks);
+         memset4 (out4, ' ' , nblanks);
          dtp->u.p.no_leading_blank = 0;
        }
-      return;
+      return true;
     } /* End of character(kind=4) internal unit code.  */
 
   /* Pad to full field width.  */
@@ -597,6 +694,10 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
 
   /* Output the decimal point.  */
   *(out++) = 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.  */
   if (nzero > 0)
@@ -630,11 +731,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
          *(out++) = expchar;
          edigits--;
        }
-#if HAVE_SNPRINTF
       snprintf (buffer, size, "%+0*d", edigits, e);
-#else
-      sprintf (buffer, "%+0*d", edigits, e);
-#endif
       memcpy (out, buffer, edigits);
     }
 
@@ -645,9 +742,7 @@ output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
       dtp->u.p.no_leading_blank = 0;
     }
 
-#undef STR
-#undef STR1
-#undef MIN_FIELD_WIDTH
+  return true;
 }
 
 
@@ -658,29 +753,46 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
 {
   char * p, 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;
-  
+
       /* 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;
+      if ((nb == 0) || dtp->u.p.g0_no_blanks)
+       {
+         if (isnan_flag)
+           nb = 3;
+         else
+           nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
+       }
       p = write_block (dtp, nb);
       if (p == NULL)
        return;
       if (nb < 3)
        {
          if (unlikely (is_char4_unit (dtp)))
-           memset4 (p, 0, '*', nb);
+           {
+             gfc_char4_t *p4 = (gfc_char4_t *) p;
+             memset4 (p4, '*', nb);
+           }
          else
            memset (p, '*', nb);
          return;
        }
 
       if (unlikely (is_char4_unit (dtp)))
-        memset4 (p, 0, ' ', nb);
+       {
+         gfc_char4_t *p4 = (gfc_char4_t *) p;
+         memset4 (p4, ' ', nb);
+       }
       else
        memset(p, ' ', nb);
 
@@ -693,7 +805,10 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
              if (nb == 3)
                {
                  if (unlikely (is_char4_unit (dtp)))
-                   memset4 (p, 0, '*', nb);
+                   {
+                     gfc_char4_t *p4 = (gfc_char4_t *) p;
+                     memset4 (p4, '*', nb);
+                   }
                  else
                    memset (p, '*', nb);
                  return;
@@ -709,24 +824,28 @@ write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit
          if (unlikely (is_char4_unit (dtp)))
            {
              gfc_char4_t *p4 = (gfc_char4_t *) p;
-             if (nb > 8)
+
+             if (nb > mark)
                /* We have room, so output 'Infinity' */
-               memcpy4 (p4, nb - 8, "Infinity", 8);
+               memcpy4 (p4 + nb - 8, "Infinity", 8);
              else
-               /* For the case of width equals 8, there is not enough room
+               /* For the case of width equals mark, there is not enough room
                   for the sign and 'Infinity' so we go with 'Inf' */
-               memcpy4 (p4, nb - 3, "Inf", 3);
+               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;
+             if (sign == S_PLUS || sign == S_MINUS)
+               {
+                 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
@@ -734,15 +853,21 @@ 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 */
+         if (sign == S_PLUS || sign == S_MINUS)
+           {
+             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)))
-           memcpy4 (p, nb - 3, "NaN", 3);
+           {
+             gfc_char4_t *p4 = (gfc_char4_t *) p;
+             memcpy4 (p4 + nb - 3, "NaN", 3);
+           }
          else
            memcpy(p + nb - 3, "NaN", 3);
        }
@@ -754,7 +879,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;\
@@ -778,6 +903,89 @@ CALCULATE_EXP(16)
 #endif
 #undef CALCULATE_EXP
 
+
+/* Define a macro to build code for write_float.  */
+
+  /* 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
+
+     -   The converted value is to be left adjusted on the field boundary
+
+     +   A sign (+ or -) always be placed before a number
+
+     *   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.  */
+
+
+#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 DTOA2L(prec,val)                               \
+snprintf (buffer, size, "%+-#.*Le", (prec), (val))
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define DTOA2Q(prec,val)                                                       \
+__qmath_(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 FDTOA2L(prec,val)                                              \
+snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define FDTOA2Q(prec,val)                             \
+__qmath_(quadmath_snprintf) (buffer, size, "%+-#.*Qf", \
+                            (prec), (val))
+#endif
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define ISFINITE2Q(val) finiteq(val)
+#endif
+#define ISFINITE2(val) isfinite(val)
+#define ISFINITE2L(val) isfinite(val)
+
+#define ISFINITE(suff,val) TOKENPASTE(ISFINITE2,suff)(val)
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define SIGNBIT2Q(val) signbitq(val)
+#endif
+#define SIGNBIT2(val) signbit(val)
+#define SIGNBIT2L(val) signbit(val)
+
+#define SIGNBIT(suff,val) TOKENPASTE(SIGNBIT2,suff)(val)
+
+
+#if defined(GFC_REAL_16_IS_FLOAT128)
+#define ISNAN2Q(val) isnanq(val)
+#endif
+#define ISNAN2(val) isnan(val)
+#define ISNAN2L(val) isnan(val)
+
+#define ISNAN(suff,val) TOKENPASTE(ISNAN2,suff)(val)
+
+
+
 /* 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:
@@ -793,36 +1001,61 @@ CALCULATE_EXP(16)
    m >= 10**d-0.5                              Ew.d[Ee]
 
    notes: for Gw.d ,  n' ' means 4 blanks
-          for Gw.dEe, n' ' means e+2 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 OUTPUT_FLOAT_FMT_G(x) \
+#define OUTPUT_FLOAT_FMT_G(x,y)                        \
 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 sign_bit, bool zero_flag, int comp_d) \
 { \
   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;\
+  fnode newf;\
+  GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
   int low, high, mid;\
   int ubound, lbound;\
-  char *p;\
+  char *p, pad = ' ';\
   int save_scale_factor, nb = 0;\
+  bool result;\
+  int nprinted, precision;\
+  volatile GFC_REAL_ ## x temp;\
 \
   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)))\
+  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;\
-      newf->u.real.e = e;\
+      newf.format = FMT_E;\
+      newf.u.real.w = w;\
+      newf.u.real.d = d - comp_d;\
+      newf.u.real.e = e;\
       nb = 0;\
+      precision = determine_precision (dtp, &newf, x);\
+      nprinted = DTOA(y,precision,m);                       \
       goto finish;\
     }\
 \
@@ -834,10 +1067,9 @@ output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
 \
   while (low <= high)\
     { \
-      GFC_REAL_ ## x temp;\
       mid = (low + high) / 2;\
 \
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
+      temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
 \
       if (m < temp)\
         { \
@@ -863,166 +1095,195 @@ output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
        }\
     }\
 \
-  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);\
-\
+  nb = e <= 0 ? 4 : e + 2;\
+  nb = nb >= w ? w - 1 : nb;\
+  newf.format = FMT_F;\
+  newf.u.real.w = w - nb;\
+  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);                                     \
 \
  finish:\
-  output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
-               edigits);\
+    result = output_float (dtp, &newf, buffer, size, nprinted, precision,\
+                          sign_bit, zero_flag);\
   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 (!result)\
+        pad = '*';\
       if (unlikely (is_char4_unit (dtp)))\
-       memset4 (p, 0, ' ', nb);\
-      else\
-       memset (p, ' ', nb);\
+       {\
+         gfc_char4_t *p4 = (gfc_char4_t *) p;\
+         memset4 (p4, pad, nb);\
+       }\
+      else \
+       memset (p, pad, nb);\
     }\
 }\
 
-OUTPUT_FLOAT_FMT_G(4)
+OUTPUT_FLOAT_FMT_G(4,)
 
-OUTPUT_FLOAT_FMT_G(8)
+OUTPUT_FLOAT_FMT_G(8,)
 
 #ifdef HAVE_GFC_REAL_10
-OUTPUT_FLOAT_FMT_G(10)
+OUTPUT_FLOAT_FMT_G(10,L)
 #endif
 
 #ifdef HAVE_GFC_REAL_16
-OUTPUT_FLOAT_FMT_G(16)
+# ifdef GFC_REAL_16_IS_FLOAT128
+OUTPUT_FLOAT_FMT_G(16,Q)
+#else
+OUTPUT_FLOAT_FMT_G(16,L)
+#endif
 #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.
-
-     #   The result will always contain a decimal point, even if no
-        digits follow it
-
-     -   The converted value is to be left adjusted on the field boundary
-
-     +   A sign (+ or -) always be placed before a number
-
-     MIN_FIELD_WIDTH  minimum field width
-
-     *   (ndigits-1) 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
+/* 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)\
+{\
+    volatile GFC_REAL_ ## x tmp, one = 1.0;\
+    tmp = * (GFC_REAL_ ## x *)source;\
+    if (ISFINITE (y,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;\
+         }\
+       nbefore = e%3;\
+       if (nbefore < 0)\
+         nbefore = 3 + nbefore;\
+      }\
+    else\
+      nprinted = -1;\
+}\
 
-#define DTOA \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-         "e", ndigits - 1, tmp);
+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.  */
 
-#define DTOAL \
-snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-         "Le", ndigits - 1, tmp);
+  switch (len)
+    {
+    case 4:
+      EN_PREC(4,)
+      break;
 
-#else
+    case 8:
+      EN_PREC(8,)
+      break;
 
-#define DTOA \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-        "e", ndigits - 1, tmp);
+#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
+    default:
+      internal_error (NULL, "bad real kind");
+    }
 
-#define DTOAL \
-sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
-        "Le", ndigits - 1, tmp);
+  if (nprinted == -1)
+    return -1;
 
-#endif
+  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;
+}
+  
 
 #define WRITE_FLOAT(x,y)\
 {\
        GFC_REAL_ ## x tmp;\
        tmp = * (GFC_REAL_ ## x *)source;\
-       sign_bit = __builtin_signbit (tmp);\
-       if (!isfinite (tmp))\
+       sign_bit = SIGNBIT (y,tmp);\
+       if (!ISFINITE (y,tmp))\
          { \
-           write_infnan (dtp, f, isnan (tmp), sign_bit);\
+           write_infnan (dtp, f, ISNAN (y,tmp), sign_bit);\
            return;\
          }\
        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 \
+       if (f->format == FMT_G)\
          output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
-                                   zero_flag, ndigits, edigits);\
+                                   zero_flag, comp_d);\
+       else\
+         {\
+           if (f->format == FMT_F)\
+             nprinted = FDTOA(y,precision,tmp);                \
+           else\
+             nprinted = DTOA(y,precision,tmp);                                 \
+           output_float (dtp, f, buffer, size, nprinted, precision,\
+                         sign_bit, zero_flag);\
+         }\
 }\
 
 /* Output a real number according to its format.  */
 
 static void
-write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
+write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
+            int len, int comp_d)
 {
-
-#if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
-# define MIN_FIELD_WIDTH 46
-#else
-# define MIN_FIELD_WIDTH 31
-#endif
-#define STR(x) STR1(x)
-#define STR1(x) #x
-
-  /* This must be large enough to accurately hold any value.  */
-  char buffer[MIN_FIELD_WIDTH+1];
-  int sign_bit, ndigits, edigits;
+  int sign_bit, nprinted;
+  int precision;  /* Precision for snprintf call.  */
   bool zero_flag;
-  size_t size;
-
-  size = MIN_FIELD_WIDTH+1;
 
-  /* 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;
-
-  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;
-    }
+  if (f->format != FMT_EN)
+    precision = determine_precision (dtp, f, len);
   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;
-    }
+    precision = determine_en_precision (dtp, f, source, len);
+
+  /* 4932 is the maximum exponent of long double and quad precision, 3
+     extra characters for the sign, the decimal point, and the
+     trailing null, and finally some extra digits depending on the
+     requested precision.  */
+  const size_t size = 4932 + 3 + precision;
+#define BUF_STACK_SZ 5000
+  char buf_stack[BUF_STACK_SZ];
+  char *buffer;
+  if (size > BUF_STACK_SZ)
+     buffer = xmalloc (size);
+  else
+     buffer = buf_stack;
 
   switch (len)
     {
@@ -1041,10 +1302,16 @@ write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
 #endif
 #ifdef HAVE_GFC_REAL_16
     case 16:
+# ifdef GFC_REAL_16_IS_FLOAT128
+      WRITE_FLOAT(16,Q)
+# else
       WRITE_FLOAT(16,L)
+# endif
       break;
 #endif
     default:
       internal_error (NULL, "bad real kind");
     }
+  if (size > BUF_STACK_SZ)
+     free (buffer);
 }