Index: gcc/testsuite/gfortran.dg/char4_iunit_1.f03
===================================================================
--- gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(revision 174673)
+++ gcc/testsuite/gfortran.dg/char4_iunit_1.f03	(working copy)
@@ -24,7 +24,7 @@ program char4_iunit_1
   write(string, *) .true., .false. , .true.
   if (string .ne. 4_" T F T                                    ") call abort
   write(string, *) 1.2345e-06, 4.2846e+10_8
-  if (string .ne. 4_"   1.23450002E-06   42846000000.000000      ") call abort
+  if (string .ne. 4_"  1.234500019E-06   42846000000.000000      ") call abort
   write(string, *) nan, inf
   if (string .ne. 4_"              NaN         Infinity    ") call abort
   write(string, '(10x,f3.1,3x,f9.1)') nan, inf
Index: gcc/testsuite/gfortran.dg/fmt_g0_6.f08
===================================================================
--- gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(revision 174673)
+++ gcc/testsuite/gfortran.dg/fmt_g0_6.f08	(working copy)
@@ -9,7 +9,7 @@ program test_g0fr
     
     call check_all(0.0_RT, 15, 2, 0)
     call check_all(0.991_RT, 15, 2, 0)
-    call check_all(0.995_RT, 15, 2, 0)
+    call check_all(0.99500000001_RT, 15, 2, 0)
     call check_all(0.996_RT, 15, 2, 0)
     call check_all(0.999_RT, 15, 2, 0)
 contains
Index: libgfortran/io/write.c
===================================================================
--- libgfortran/io/write.c	(revision 174673)
+++ libgfortran/io/write.c	(working copy)
@@ -1155,35 +1155,35 @@ write_z (st_parameter_dt *dtp, const fnode *f, con
 void
 write_d (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_e (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_f (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_en (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
 void
 write_es (st_parameter_dt *dtp, const fnode *f, const char *p, int len)
 {
-  write_float (dtp, f, p, len, 0);
+  write_float (dtp, f, p, len);
 }
 
 
@@ -1476,7 +1476,7 @@ write_real (st_parameter_dt *dtp, const char *sour
   int org_scale = dtp->u.p.scale_factor;
   dtp->u.p.scale_factor = 1;
   set_fnode_default (dtp, &f, length);
-  write_float (dtp, &f, source , length, 1);
+  write_float (dtp, &f, source , length);
   dtp->u.p.scale_factor = org_scale;
 }
 
@@ -1487,19 +1487,11 @@ void
 write_real_g0 (st_parameter_dt *dtp, const char *source, int length, int d)
 {
   fnode f;
-  int comp_d; 
   set_fnode_default (dtp, &f, length);
   if (d > 0)
     f.u.real.d = d;
-
-  /* Compensate for extra digits when using scale factor, d is not
-     specified, and the magnitude is such that E editing is used.  */
-  if (dtp->u.p.scale_factor > 0 && d == 0)
-    comp_d = 1;
-  else
-    comp_d = 0;
   dtp->u.p.g0_no_blanks = 1;
-  write_float (dtp, &f, source , length, comp_d);
+  write_float (dtp, &f, source , length);
   dtp->u.p.g0_no_blanks = 0;
 }
 
Index: libgfortran/io/write_float.def
===================================================================
--- libgfortran/io/write_float.def	(revision 174673)
+++ libgfortran/io/write_float.def	(working copy)
@@ -59,8 +59,22 @@ calculate_sign (st_parameter_dt *dtp, int negative
 }
 
 
-/* Output a real number according to its format which is FMT_G free.  */
+/* Output a real number according to its format.
 
+   What this does:
+
+	- Calculates the number of digits to emit before and after the
+	  decimal point.
+	- Performs the required rounding.
+	- Calculates the number of leading blanks to emit.
+	- Calculates the number of trailing blanks to emit.
+	- Calculates the exponent format if any.
+	- Determines if the sign should be displayed.
+	- Allocates the write buffer according to the width.
+	- Emits the formatted number according to the above.
+
+   buffer contains the digit string when we enter this function.  */
+
 static try
 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
 	      int sign_bit, bool zero_flag, int ndigits, int edigits)
@@ -78,8 +92,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
   int nafter;
   /* Number of zeros after the decimal point, whatever the precision.  */
   int nzero_real;
+  /* Flag for optional leading zero, if there is room.  */
   int leadzero;
-  int nblanks;
+  /* Number of leading and trailing blanks.  */
+  int lblanks, tblanks; 
   sign_t sign;
 
   ft = f->format;
@@ -87,7 +103,6 @@ output_float (st_parameter_dt *dtp, const fnode *f
   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.  */
@@ -142,19 +157,50 @@ output_float (st_parameter_dt *dtp, const fnode *f
       expchar = 0;
       break;
 
+    case FMT_G:
+      nbefore = 0;
+      nzero = 0;
+      nafter = d;
+      expchar = 0;
+
+      if (e == d)
+	{
+	  if (e > 0)
+	    {
+	      nbefore = e;
+	      nafter = d > nbefore ? d - nbefore : 0;
+	    }
+	  else
+	    expchar = 'E';
+	  break;
+	}
+      else if (e < d && e >= 0)
+	{
+	  nbefore = e;
+	  nafter = d > nbefore ? d - nbefore : 0;
+	  break;
+	}
+      else if (e >= -d && e < 0 && f->u.real.e == -1)
+	{
+	  nafter = d;
+	  expchar = 'E';
+	}
+
+    /* Fall through.  */
+
     case FMT_E:
     case FMT_D:
-      i = dtp->u.p.scale_factor;
       if (d <= 0 && p == 0)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
-			  "greater than zero in format specifier 'E' or 'D'");
+			  "greater than zero in format specifier 'D', 'E', "
+			  "or 'G'");
 	  return FAILURE;
 	}
       if (p <= -d || p >= d + 2)
 	{
 	  generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
-			  "out of range in format specifier 'E' or 'D'");
+			  "out of range in format specifier 'E', 'D', or 'G'");
 	  return FAILURE;
 	}
 
@@ -179,10 +225,11 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  nafter = d;
 	}
 
-      if (ft == FMT_E)
+      if (ft == FMT_E || ft == FMT_G)
 	expchar = 'E';
       else
 	expchar = 'D';
+      
       break;
 
     case FMT_EN:
@@ -224,6 +271,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 
   /* Round the value.  The value being rounded is an unsigned magnitude.
      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
+  rchar = '5';
   switch (dtp->u.p.current_unit->round_status)
     {
       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
@@ -275,6 +323,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   rchar = '0';
   if (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++)
     {
@@ -284,7 +333,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
   goto skip;
     
   do_rnd:
- 
+
   if (nbefore + nafter == 0)
     {
       ndigits = 0;
@@ -321,6 +370,7 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	         zero.  */
 	      digits--;
 	      digits[0] = '1';
+
 	      if (ft == FMT_F)
 		{
 		  if (nzero > 0)
@@ -331,6 +381,27 @@ output_float (st_parameter_dt *dtp, const fnode *f
 		  else
 		    nbefore++;
 		}
+	      else if (ft == FMT_G)
+		{
+		  if (e < d && e >= 0)
+		    {
+		      nbefore++;
+		      nafter--;
+		    }
+		  else if (e == d)
+		    {
+		      nbefore -= d;
+		      nafter += d;
+		      e++;
+		      expchar = 'E';
+		    }
+		  else
+		    {
+		      e++;
+		      if (e == 0)
+		        expchar = 0;
+		    }
+		}
 	      else if (ft == FMT_EN)
 		{
 		  nbefore++;
@@ -346,11 +417,68 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	}
     }
 
+  /* 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. 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.  */
+  for (i = 0; i < ndigits; i++)
+    {
+      if (digits[i] != '0')
+	break;
+    }
+  if (i == ndigits)
+    {
+      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)
+        sign = calculate_sign (dtp, sign_bit);
+      else
+	sign = calculate_sign (dtp, 0);
+    }
+
   skip:
 
+  /* Pick a field size if none was specified, taking into account small
+     values that may have been rounded to zero.  */
+  if (w <= 0)
+    {
+      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;
+	}
+    }
+  
+  /* Figure out the leading and trailing blanks.  */
+  lblanks = tblanks = 0;
+
+  /* G formatting is a special case of several things.  */
+  if (ft == FMT_G)
+    {
+      if (zero_flag && nbefore == 0 && nafter > 0)
+	{
+	  nbefore++;
+	  nafter--;
+	}
+      tblanks = edigits = f->u.real.e == -1 ? 4 : f->u.real.e + 2;
+      lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+		     (sign != S_NONE ? 1 : 0));
+    }
+  else
+    {
+      lblanks = w - (nbefore + nzero + nafter + 1 + e);
+      if (sign != S_NONE)
+	lblanks--;
+    }
+
   /* Calculate the format of the exponent field.  */
   if (expchar)
     {
+      tblanks = 0;
       edigits = 1;
       for (i = abs (e); i >= 10; i /= 10)
 	edigits++;
@@ -379,76 +507,57 @@ output_float (st_parameter_dt *dtp, const fnode *f
   else
     edigits = 0;
 
-  /* 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.  */
-  for (i = 0; i < ndigits; i++)
+  if (ft == FMT_G)
     {
-      if (digits[i] != '0')
-	break;
-    }
-
-  /* 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)
-    {
-      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)
-        sign = calculate_sign (dtp, sign_bit);
-      else
-	sign = calculate_sign (dtp, 0);
-    }
-
-  /* Pick a field size if none was specified, taking into account small
-     values that may have been rounded to zero.  */
-  if (w <= 0)
-    {
-      if (zero_flag)
-	w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
-      else
+      if (dtp->u.p.g0_no_blanks)
 	{
-	  w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
-	  w = w == 1 ? 2 : w;
+	  w = w - tblanks - lblanks;
+	  lblanks = tblanks = 0;
 	}
+      else if (expchar == 0)
+	{
+	  lblanks = w - (nbefore + nzero + nafter + 1 + tblanks +
+			 (sign != S_NONE ? 1 : 0));
+	  if (dtp->u.p.no_leading_blank)
+	    {
+	      tblanks += lblanks;
+	      lblanks = 0; 
+	    }
+	}
     }
-  
-  /* Work out how much padding is needed.  */
-  nblanks = w - (nbefore + nzero + nafter + edigits + 1);
-  if (sign != S_NONE)
-    nblanks--;
+  else
+    lblanks = w - (nbefore + nzero + nafter + 1 + edigits +
+		 (sign != S_NONE ? 1 : 0));
 
-  if (dtp->u.p.g0_no_blanks)
-    {
-      w -= nblanks;
-      nblanks = 0;
-    }
-
   /* Create the ouput buffer.  */
-  out = write_block (dtp, w);
+  if ((ft == FMT_G && dtp->u.p.g0_no_blanks)
+      || (ft == FMT_F && f->u.real.w == 0))
+    out = write_block (dtp, w);
+  else
+    out = write_block (dtp, f->u.real.w);
+
   if (out == NULL)
     return FAILURE;
 
   /* Check the value fits in the specified field width.  */
-  if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
+  if (lblanks < 0 || tblanks < 0 || edigits == -1 || w == 1
+      || (w == 2 && sign != S_NONE))
     {
       if (unlikely (is_char4_unit (dtp)))
 	{
 	  gfc_char4_t *out4 = (gfc_char4_t *) out;
-	  memset4 (out4, '*', w);
+	  memset4 (out4, '*', f->u.real.w);
 	  return FAILURE;
 	}
-      star_fill (out, w);
+      star_fill (out, f->u.real.w);
       return FAILURE;
     }
 
   /* See if we have space for a zero before the decimal point.  */
-  if (nbefore == 0 && nblanks > 0)
+  if (nbefore == 0 && lblanks > 0)
     {
       leadzero = 1;
-      nblanks--;
+      lblanks--;
     }
   else
     leadzero = 0;
@@ -461,10 +570,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       gfc_char4_t *out4 = (gfc_char4_t *) out;
       /* Pad to full field width.  */
 
-      if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+      if ((lblanks > 0) && !dtp->u.p.no_leading_blank)
 	{
-	  memset4 (out4, ' ', nblanks);
-	  out4 += nblanks;
+	  memset4 (out4, ' ', lblanks);
+	  out4 += lblanks;
 	}
 
       /* Output the initial sign (if any).  */
@@ -539,21 +648,17 @@ output_float (st_parameter_dt *dtp, const fnode *f
 	  memcpy4 (out4, buffer, edigits);
 	}
 
-      if (dtp->u.p.no_leading_blank)
-	{
-	  out4 += edigits;
-	  memset4 (out4, ' ' , nblanks);
-	  dtp->u.p.no_leading_blank = 0;
-	}
+      /* Emit trailing blanks if needed.  */
+      memset4 (out4, ' ' , tblanks);
       return SUCCESS;
     } /* End of character(kind=4) internal unit code.  */
 
-  /* Pad to full field width.  */
+  /* Pad leading to full field width.  */
 
-  if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
+  if (lblanks > 0 && !dtp->u.p.no_leading_blank)
     {
-      memset (out, ' ', nblanks);
-      out += nblanks;
+      memset (out, ' ', lblanks);
+      out += lblanks;
     }
 
   /* Output the initial sign (if any).  */
@@ -627,13 +732,10 @@ output_float (st_parameter_dt *dtp, const fnode *f
       memcpy (out, buffer, edigits);
     }
 
-  if (dtp->u.p.no_leading_blank)
-    {
-      out += edigits;
-      memset( out , ' ' , nblanks );
-      dtp->u.p.no_leading_blank = 0;
-    }
-
+  /* Emit trailing blanks if needed.  */
+  if (tblanks)
+    memset( out , ' ' , tblanks );
+  
 #undef STR
 #undef STR1
 #undef MIN_FIELD_WIDTH
@@ -770,185 +872,6 @@ write_infnan (st_parameter_dt *dtp, const fnode *f
     }
 }
 
-
-/* Returns the value of 10**d.  */
-
-#define CALCULATE_EXP(x) \
-inline static GFC_REAL_ ## x \
-calculate_exp_ ## x  (int d)\
-{\
-  int i;\
-  GFC_REAL_ ## x r = 1.0;\
-  for (i = 0; i< (d >= 0 ? d : -d); i++)\
-    r *= 10;\
-  r = (d >= 0) ? r : 1.0 / r;\
-  return r;\
-}
-
-CALCULATE_EXP(4)
-
-CALCULATE_EXP(8)
-
-#ifdef HAVE_GFC_REAL_10
-CALCULATE_EXP(10)
-#endif
-
-#ifdef HAVE_GFC_REAL_16
-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]
-
-   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 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 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, r = 0.5;\
-  int low, high, mid;\
-  int ubound, lbound;\
-  char *p, pad = ' ';\
-  int save_scale_factor, nb = 0;\
-  try result;\
-\
-  save_scale_factor = dtp->u.p.scale_factor;\
-  newf = (fnode *) get_mem (sizeof (fnode));\
-\
-  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;\
-    }\
-\
-  rexp_d = calculate_exp_ ## x (-d);\
-  if ((m > 0.0 && ((m < 0.1 - 0.1 * r * rexp_d) || (rexp_d * (m + r) >= 1.0)))\
-      || ((m == 0.0) && !(compile_options.allow_std\
-			  & (GFC_STD_F2003 | GFC_STD_F2008))))\
-    { \
-      newf->format = FMT_E;\
-      newf->u.real.w = w;\
-      newf->u.real.d = d - comp_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)\
-    { \
-      volatile GFC_REAL_ ## x temp;\
-      mid = (low + high) / 2;\
-\
-      temp = (calculate_exp_ ## x (mid - 1) * (1 - r * 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;\
-	}\
-    }\
-\
-  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;\
-\
- finish:\
-  result = 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 (result == FAILURE)\
-        pad = '*';\
-      if (unlikely (is_char4_unit (dtp)))\
-	{\
-	  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(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, snprintf is used to print to buffer the
@@ -1003,19 +926,15 @@ __qmath_(quadmath_snprintf) (buffer, sizeof buffer
 \
 	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, comp_d);\
+	output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
+		      edigits);\
 }\
 
 /* Output a real number according to its format.  */
 
 static void
 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
-            int len, int comp_d)
+            int len)
 {
 
 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18