OSDN Git Service

3709b206620080016d70a844f0780ca5a2ff8878
[pf3gnuchains/gcc-fork.git] / libgfortran / io / write_float.def
1 /* Copyright (C) 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
2    Contributed by Andy Vaught
3    Write float code factoring to this file by Jerry DeLisle   
4    F2003 I/O support contributed by Jerry DeLisle
5
6 This file is part of the GNU Fortran runtime library (libgfortran).
7
8 Libgfortran is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 3, or (at your option)
11 any later version.
12
13 Libgfortran is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 GNU General Public License for more details.
17
18 Under Section 7 of GPL version 3, you are granted additional
19 permissions described in the GCC Runtime Library Exception, version
20 3.1, as published by the Free Software Foundation.
21
22 You should have received a copy of the GNU General Public License and
23 a copy of the GCC Runtime Library Exception along with this program;
24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
25 <http://www.gnu.org/licenses/>.  */
26
27 #include "config.h"
28
29 typedef enum
30 { S_NONE, S_MINUS, S_PLUS }
31 sign_t;
32
33 /* Given a flag that indicates if a value is negative or not, return a
34    sign_t that gives the sign that we need to produce.  */
35
36 static sign_t
37 calculate_sign (st_parameter_dt *dtp, int negative_flag)
38 {
39   sign_t s = S_NONE;
40
41   if (negative_flag)
42     s = S_MINUS;
43   else
44     switch (dtp->u.p.sign_status)
45       {
46       case SIGN_SP:     /* Show sign. */
47         s = S_PLUS;
48         break;
49       case SIGN_SS:     /* Suppress sign. */
50         s = S_NONE;
51         break;
52       case SIGN_S:      /* Processor defined. */
53       case SIGN_UNSPECIFIED:
54         s = options.optional_plus ? S_PLUS : S_NONE;
55         break;
56       }
57
58   return s;
59 }
60
61
62 /* Output a real number according to its format which is FMT_G free.  */
63
64 static try
65 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size, 
66               int sign_bit, bool zero_flag, int ndigits, int edigits)
67 {
68   char *out;
69   char *digits;
70   int e;
71   char expchar, rchar;
72   format_token ft;
73   int w;
74   int d;
75   /* Number of digits before the decimal point.  */
76   int nbefore;
77   /* Number of zeros after the decimal point.  */
78   int nzero;
79   /* Number of digits after the decimal point.  */
80   int nafter;
81   /* Number of zeros after the decimal point, whatever the precision.  */
82   int nzero_real;
83   int leadzero;
84   int nblanks;
85   int i;
86   sign_t sign;
87
88   ft = f->format;
89   w = f->u.real.w;
90   d = f->u.real.d;
91
92   rchar = '5';
93   nzero_real = -1;
94
95   /* We should always know the field width and precision.  */
96   if (d < 0)
97     internal_error (&dtp->common, "Unspecified precision");
98
99   sign = calculate_sign (dtp, sign_bit);
100   
101   /* The following code checks the given string has punctuation in the correct
102      places.  Uncomment if needed for debugging.
103      if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
104                     || buffer[ndigits + 2] != 'e'))
105        internal_error (&dtp->common, "printf is broken");  */
106
107   /* Read the exponent back in.  */
108   e = atoi (&buffer[ndigits + 3]) + 1;
109
110   /* Make sure zero comes out as 0.0e0.   */
111   if (zero_flag)
112     e = 0;
113
114   /* Normalize the fractional component.  */
115   buffer[2] = buffer[1];
116   digits = &buffer[2];
117
118   /* Figure out where to place the decimal point.  */
119   switch (ft)
120     {
121     case FMT_F:
122       if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0)
123         {
124           memmove (digits + 1, digits, ndigits - 1);
125           digits[0] = '0';
126           e++;
127         }
128
129       nbefore = e + dtp->u.p.scale_factor;
130       if (nbefore < 0)
131         {
132           nzero = -nbefore;
133           nzero_real = nzero;
134           if (nzero > d)
135             nzero = d;
136           nafter = d - nzero;
137           nbefore = 0;
138         }
139       else
140         {
141           nzero = 0;
142           nafter = d;
143         }
144       expchar = 0;
145       break;
146
147     case FMT_E:
148     case FMT_D:
149       i = dtp->u.p.scale_factor;
150       if (d <= 0 && i == 0)
151         {
152           generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
153                           "greater than zero in format specifier 'E' or 'D'");
154           return FAILURE;
155         }
156       if (i <= -d || i >= d + 2)
157         {
158           generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
159                           "out of range in format specifier 'E' or 'D'");
160           return FAILURE;
161         }
162
163       if (!zero_flag)
164         e -= i;
165       if (i < 0)
166         {
167           nbefore = 0;
168           nzero = -i;
169           nafter = d + i;
170         }
171       else if (i > 0)
172         {
173           nbefore = i;
174           nzero = 0;
175           nafter = (d - i) + 1;
176         }
177       else /* i == 0 */
178         {
179           nbefore = 0;
180           nzero = 0;
181           nafter = d;
182         }
183
184       if (ft == FMT_E)
185         expchar = 'E';
186       else
187         expchar = 'D';
188       break;
189
190     case FMT_EN:
191       /* The exponent must be a multiple of three, with 1-3 digits before
192          the decimal point.  */
193       if (!zero_flag)
194         e--;
195       if (e >= 0)
196         nbefore = e % 3;
197       else
198         {
199           nbefore = (-e) % 3;
200           if (nbefore != 0)
201             nbefore = 3 - nbefore;
202         }
203       e -= nbefore;
204       nbefore++;
205       nzero = 0;
206       nafter = d;
207       expchar = 'E';
208       break;
209
210     case FMT_ES:
211       if (!zero_flag)
212         e--;
213       nbefore = 1;
214       nzero = 0;
215       nafter = d;
216       expchar = 'E';
217       break;
218
219     default:
220       /* Should never happen.  */
221       internal_error (&dtp->common, "Unexpected format token");
222     }
223
224   /* Round the value.  The value being rounded is an unsigned magnitude.
225      The ROUND_COMPATIBLE is rounding away from zero when there is a tie.  */
226   switch (dtp->u.p.current_unit->round_status)
227     {
228       case ROUND_ZERO: /* Do nothing and truncation occurs.  */
229         goto skip;
230       case ROUND_UP:
231         if (sign_bit)
232           goto skip;
233         rchar = '0';
234         break;
235       case ROUND_DOWN:
236         if (!sign_bit)
237           goto skip;
238         rchar = '0';
239         break;
240       case ROUND_NEAREST:
241         /* Round compatible unless there is a tie. A tie is a 5 with
242            all trailing zero's.  */
243         i = nafter + nbefore;
244         if (digits[i] == '5')
245           {
246             for(i++ ; i < ndigits; i++)
247               {
248                 if (digits[i] != '0')
249                   goto do_rnd;
250               }
251             /* It is a  tie so round to even.  */
252             switch (digits[nafter + nbefore - 1])
253               {
254                 case '1':
255                 case '3':
256                 case '5':
257                 case '7':
258                 case '9':
259                   /* If odd, round away from zero to even.  */
260                   break;
261                 default:
262                   /* If even, skip rounding, truncate to even.  */
263                   goto skip;
264               }
265           }
266          /* Fall through.  */ 
267       case ROUND_PROCDEFINED:
268       case ROUND_UNSPECIFIED:
269       case ROUND_COMPATIBLE:
270         rchar = '5';
271         /* Just fall through and do the actual rounding.  */
272     }
273     
274   do_rnd:
275  
276   if (nbefore + nafter == 0)
277     {
278       ndigits = 0;
279       if (nzero_real == d && digits[0] >= rchar)
280         {
281           /* We rounded to zero but shouldn't have */
282           nzero--;
283           nafter = 1;
284           digits[0] = '1';
285           ndigits = 1;
286         }
287     }
288   else if (nbefore + nafter < ndigits)
289     {
290       ndigits = nbefore + nafter;
291       i = ndigits;
292       if (digits[i] >= rchar)
293         {
294           /* Propagate the carry.  */
295           for (i--; i >= 0; i--)
296             {
297               if (digits[i] != '9')
298                 {
299                   digits[i]++;
300                   break;
301                 }
302               digits[i] = '0';
303             }
304
305           if (i < 0)
306             {
307               /* The carry overflowed.  Fortunately we have some spare
308                  space at the start of the buffer.  We may discard some
309                  digits, but this is ok because we already know they are
310                  zero.  */
311               digits--;
312               digits[0] = '1';
313               if (ft == FMT_F)
314                 {
315                   if (nzero > 0)
316                     {
317                       nzero--;
318                       nafter++;
319                     }
320                   else
321                     nbefore++;
322                 }
323               else if (ft == FMT_EN)
324                 {
325                   nbefore++;
326                   if (nbefore == 4)
327                     {
328                       nbefore = 1;
329                       e += 3;
330                     }
331                 }
332               else
333                 e++;
334             }
335         }
336     }
337
338   skip:
339
340   /* Calculate the format of the exponent field.  */
341   if (expchar)
342     {
343       edigits = 1;
344       for (i = abs (e); i >= 10; i /= 10)
345         edigits++;
346
347       if (f->u.real.e < 0)
348         {
349           /* Width not specified.  Must be no more than 3 digits.  */
350           if (e > 999 || e < -999)
351             edigits = -1;
352           else
353             {
354               edigits = 4;
355               if (e > 99 || e < -99)
356                 expchar = ' ';
357             }
358         }
359       else
360         {
361           /* Exponent width specified, check it is wide enough.  */
362           if (edigits > f->u.real.e)
363             edigits = -1;
364           else
365             edigits = f->u.real.e + 2;
366         }
367     }
368   else
369     edigits = 0;
370
371   /* Scan the digits string and count the number of zeros.  If we make it
372      all the way through the loop, we know the value is zero after the
373      rounding completed above.  */
374   for (i = 0; i < ndigits; i++)
375     {
376       if (digits[i] != '0')
377         break;
378     }
379
380   /* To format properly, we need to know if the rounded result is zero and if
381      so, we set the zero_flag which may have been already set for
382      actual zero.  */
383   if (i == ndigits)
384     {
385       zero_flag = true;
386       /* The output is zero, so set the sign according to the sign bit unless
387          -fno-sign-zero was specified.  */
388       if (compile_options.sign_zero == 1)
389         sign = calculate_sign (dtp, sign_bit);
390       else
391         sign = calculate_sign (dtp, 0);
392     }
393
394   /* Pick a field size if none was specified, taking into account small
395      values that may have been rounded to zero.  */
396   if (w <= 0)
397     {
398       if (zero_flag)
399         w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
400       else
401         {
402           w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
403           w = w == 1 ? 2 : w;
404         }
405     }
406   
407   /* Work out how much padding is needed.  */
408   nblanks = w - (nbefore + nzero + nafter + edigits + 1);
409   if (sign != S_NONE)
410     nblanks--;
411
412   if (dtp->u.p.g0_no_blanks)
413     {
414       w -= nblanks;
415       nblanks = 0;
416     }
417
418   /* Create the ouput buffer.  */
419   out = write_block (dtp, w);
420   if (out == NULL)
421     return FAILURE;
422
423   /* Check the value fits in the specified field width.  */
424   if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
425     {
426       if (unlikely (is_char4_unit (dtp)))
427         {
428           gfc_char4_t *out4 = (gfc_char4_t *) out;
429           memset4 (out4, '*', w);
430           return FAILURE;
431         }
432       star_fill (out, w);
433       return FAILURE;
434     }
435
436   /* See if we have space for a zero before the decimal point.  */
437   if (nbefore == 0 && nblanks > 0)
438     {
439       leadzero = 1;
440       nblanks--;
441     }
442   else
443     leadzero = 0;
444
445   /* For internal character(kind=4) units, we duplicate the code used for
446      regular output slightly modified.  This needs to be maintained
447      consistent with the regular code that follows this block.  */
448   if (unlikely (is_char4_unit (dtp)))
449     {
450       gfc_char4_t *out4 = (gfc_char4_t *) out;
451       /* Pad to full field width.  */
452
453       if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
454         {
455           memset4 (out4, ' ', nblanks);
456           out4 += nblanks;
457         }
458
459       /* Output the initial sign (if any).  */
460       if (sign == S_PLUS)
461         *(out4++) = '+';
462       else if (sign == S_MINUS)
463         *(out4++) = '-';
464
465       /* Output an optional leading zero.  */
466       if (leadzero)
467         *(out4++) = '0';
468
469       /* Output the part before the decimal point, padding with zeros.  */
470       if (nbefore > 0)
471         {
472           if (nbefore > ndigits)
473             {
474               i = ndigits;
475               memcpy4 (out4, digits, i);
476               ndigits = 0;
477               while (i < nbefore)
478                 out4[i++] = '0';
479             }
480           else
481             {
482               i = nbefore;
483               memcpy4 (out4, digits, i);
484               ndigits -= i;
485             }
486
487           digits += i;
488           out4 += nbefore;
489         }
490
491       /* Output the decimal point.  */
492       *(out4++) = dtp->u.p.current_unit->decimal_status
493                     == DECIMAL_POINT ? '.' : ',';
494
495       /* Output leading zeros after the decimal point.  */
496       if (nzero > 0)
497         {
498           for (i = 0; i < nzero; i++)
499             *(out4++) = '0';
500         }
501
502       /* Output digits after the decimal point, padding with zeros.  */
503       if (nafter > 0)
504         {
505           if (nafter > ndigits)
506             i = ndigits;
507           else
508             i = nafter;
509
510           memcpy4 (out4, digits, i);
511           while (i < nafter)
512             out4[i++] = '0';
513
514           digits += i;
515           ndigits -= i;
516           out4 += nafter;
517         }
518
519       /* Output the exponent.  */
520       if (expchar)
521         {
522           if (expchar != ' ')
523             {
524               *(out4++) = expchar;
525               edigits--;
526             }
527 #if HAVE_SNPRINTF
528           snprintf (buffer, size, "%+0*d", edigits, e);
529 #else
530           sprintf (buffer, "%+0*d", edigits, e);
531 #endif
532           memcpy4 (out4, buffer, edigits);
533         }
534
535       if (dtp->u.p.no_leading_blank)
536         {
537           out4 += edigits;
538           memset4 (out4, ' ' , nblanks);
539           dtp->u.p.no_leading_blank = 0;
540         }
541       return SUCCESS;
542     } /* End of character(kind=4) internal unit code.  */
543
544   /* Pad to full field width.  */
545
546   if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
547     {
548       memset (out, ' ', nblanks);
549       out += nblanks;
550     }
551
552   /* Output the initial sign (if any).  */
553   if (sign == S_PLUS)
554     *(out++) = '+';
555   else if (sign == S_MINUS)
556     *(out++) = '-';
557
558   /* Output an optional leading zero.  */
559   if (leadzero)
560     *(out++) = '0';
561
562   /* Output the part before the decimal point, padding with zeros.  */
563   if (nbefore > 0)
564     {
565       if (nbefore > ndigits)
566         {
567           i = ndigits;
568           memcpy (out, digits, i);
569           ndigits = 0;
570           while (i < nbefore)
571             out[i++] = '0';
572         }
573       else
574         {
575           i = nbefore;
576           memcpy (out, digits, i);
577           ndigits -= i;
578         }
579
580       digits += i;
581       out += nbefore;
582     }
583
584   /* Output the decimal point.  */
585   *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
586
587   /* Output leading zeros after the decimal point.  */
588   if (nzero > 0)
589     {
590       for (i = 0; i < nzero; i++)
591         *(out++) = '0';
592     }
593
594   /* Output digits after the decimal point, padding with zeros.  */
595   if (nafter > 0)
596     {
597       if (nafter > ndigits)
598         i = ndigits;
599       else
600         i = nafter;
601
602       memcpy (out, digits, i);
603       while (i < nafter)
604         out[i++] = '0';
605
606       digits += i;
607       ndigits -= i;
608       out += nafter;
609     }
610
611   /* Output the exponent.  */
612   if (expchar)
613     {
614       if (expchar != ' ')
615         {
616           *(out++) = expchar;
617           edigits--;
618         }
619 #if HAVE_SNPRINTF
620       snprintf (buffer, size, "%+0*d", edigits, e);
621 #else
622       sprintf (buffer, "%+0*d", edigits, e);
623 #endif
624       memcpy (out, buffer, edigits);
625     }
626
627   if (dtp->u.p.no_leading_blank)
628     {
629       out += edigits;
630       memset( out , ' ' , nblanks );
631       dtp->u.p.no_leading_blank = 0;
632     }
633
634 #undef STR
635 #undef STR1
636 #undef MIN_FIELD_WIDTH
637   return SUCCESS;
638 }
639
640
641 /* Write "Infinite" or "Nan" as appropriate for the given format.  */
642
643 static void
644 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
645 {
646   char * p, fin;
647   int nb = 0;
648   sign_t sign;
649   int mark;
650
651   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
652     {
653       sign = calculate_sign (dtp, sign_bit);
654       mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
655
656       nb =  f->u.real.w;
657   
658       /* If the field width is zero, the processor must select a width 
659          not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
660      
661       if (nb == 0)
662         {
663           if (isnan_flag)
664             nb = 3;
665           else
666             nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
667         }
668       p = write_block (dtp, nb);
669       if (p == NULL)
670         return;
671       if (nb < 3)
672         {
673           if (unlikely (is_char4_unit (dtp)))
674             {
675               gfc_char4_t *p4 = (gfc_char4_t *) p;
676               memset4 (p4, '*', nb);
677             }
678           else
679             memset (p, '*', nb);
680           return;
681         }
682
683       if (unlikely (is_char4_unit (dtp)))
684         {
685           gfc_char4_t *p4 = (gfc_char4_t *) p;
686           memset4 (p4, ' ', nb);
687         }
688       else
689         memset(p, ' ', nb);
690
691       if (!isnan_flag)
692         {
693           if (sign_bit)
694             {
695               /* If the sign is negative and the width is 3, there is
696                  insufficient room to output '-Inf', so output asterisks */
697               if (nb == 3)
698                 {
699                   if (unlikely (is_char4_unit (dtp)))
700                     {
701                       gfc_char4_t *p4 = (gfc_char4_t *) p;
702                       memset4 (p4, '*', nb);
703                     }
704                   else
705                     memset (p, '*', nb);
706                   return;
707                 }
708               /* The negative sign is mandatory */
709               fin = '-';
710             }    
711           else
712             /* The positive sign is optional, but we output it for
713                consistency */
714             fin = '+';
715             
716           if (unlikely (is_char4_unit (dtp)))
717             {
718               gfc_char4_t *p4 = (gfc_char4_t *) p;
719
720               if (nb > mark)
721                 /* We have room, so output 'Infinity' */
722                 memcpy4 (p4 + nb - 8, "Infinity", 8);
723               else
724                 /* For the case of width equals mark, there is not enough room
725                    for the sign and 'Infinity' so we go with 'Inf' */
726                 memcpy4 (p4 + nb - 3, "Inf", 3);
727
728               if (sign == S_PLUS || sign == S_MINUS)
729                 {
730                   if (nb < 9 && nb > 3)
731                     /* Put the sign in front of Inf */
732                     p4[nb - 4] = (gfc_char4_t) fin;
733                   else if (nb > 8)
734                     /* Put the sign in front of Infinity */
735                     p4[nb - 9] = (gfc_char4_t) fin;
736                 }
737               return;
738             }
739
740           if (nb > mark)
741             /* We have room, so output 'Infinity' */
742             memcpy(p + nb - 8, "Infinity", 8);
743           else
744             /* For the case of width equals 8, there is not enough room
745                for the sign and 'Infinity' so we go with 'Inf' */
746             memcpy(p + nb - 3, "Inf", 3);
747
748           if (sign == S_PLUS || sign == S_MINUS)
749             {
750               if (nb < 9 && nb > 3)
751                 p[nb - 4] = fin;  /* Put the sign in front of Inf */
752               else if (nb > 8)
753                 p[nb - 9] = fin;  /* Put the sign in front of Infinity */
754             }
755         }
756       else
757         {
758           if (unlikely (is_char4_unit (dtp)))
759             {
760               gfc_char4_t *p4 = (gfc_char4_t *) p;
761               memcpy4 (p4 + nb - 3, "NaN", 3);
762             }
763           else
764             memcpy(p + nb - 3, "NaN", 3);
765         }
766       return;
767     }
768 }
769
770
771 /* Returns the value of 10**d.  */
772
773 #define CALCULATE_EXP(x) \
774 inline static GFC_REAL_ ## x \
775 calculate_exp_ ## x  (int d)\
776 {\
777   int i;\
778   GFC_REAL_ ## x r = 1.0;\
779   for (i = 0; i< (d >= 0 ? d : -d); i++)\
780     r *= 10;\
781   r = (d >= 0) ? r : 1.0 / r;\
782   return r;\
783 }
784
785 CALCULATE_EXP(4)
786
787 CALCULATE_EXP(8)
788
789 #ifdef HAVE_GFC_REAL_10
790 CALCULATE_EXP(10)
791 #endif
792
793 #ifdef HAVE_GFC_REAL_16
794 CALCULATE_EXP(16)
795 #endif
796 #undef CALCULATE_EXP
797
798 /* Generate corresponding I/O format for FMT_G and output.
799    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
800    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
801
802    Data Magnitude                              Equivalent Conversion
803    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
804    m = 0                                       F(w-n).(d-1), n' '
805    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
806    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
807    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
808    ................                           ..........
809    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
810    m >= 10**d-0.5                              Ew.d[Ee]
811
812    notes: for Gw.d ,  n' ' means 4 blanks
813           for Gw.dEe, n' ' means e+2 blanks  */
814
815 #define OUTPUT_FLOAT_FMT_G(x) \
816 static void \
817 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
818                       GFC_REAL_ ## x m, char *buffer, size_t size, \
819                       int sign_bit, bool zero_flag, int ndigits, int edigits) \
820 { \
821   int e = f->u.real.e;\
822   int d = f->u.real.d;\
823   int w = f->u.real.w;\
824   fnode *newf;\
825   GFC_REAL_ ## x rexp_d;\
826   int low, high, mid;\
827   int ubound, lbound;\
828   char *p, pad = ' ';\
829   int save_scale_factor, nb = 0;\
830   try result;\
831 \
832   save_scale_factor = dtp->u.p.scale_factor;\
833   newf = (fnode *) get_mem (sizeof (fnode));\
834 \
835   rexp_d = calculate_exp_ ## x (-d);\
836   if ((m > 0.0 && m < 0.1 - 0.05 * rexp_d) || (rexp_d * (m + 0.5) >= 1.0) ||\
837       ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
838     { \
839       newf->format = FMT_E;\
840       newf->u.real.w = w;\
841       newf->u.real.d = d;\
842       newf->u.real.e = e;\
843       nb = 0;\
844       goto finish;\
845     }\
846 \
847   mid = 0;\
848   low = 0;\
849   high = d + 1;\
850   lbound = 0;\
851   ubound = d + 1;\
852 \
853   while (low <= high)\
854     { \
855       GFC_REAL_ ## x temp;\
856       mid = (low + high) / 2;\
857 \
858       temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
859 \
860       if (m < temp)\
861         { \
862           ubound = mid;\
863           if (ubound == lbound + 1)\
864             break;\
865           high = mid - 1;\
866         }\
867       else if (m > temp)\
868         { \
869           lbound = mid;\
870           if (ubound == lbound + 1)\
871             { \
872               mid ++;\
873               break;\
874             }\
875           low = mid + 1;\
876         }\
877       else\
878         {\
879           mid++;\
880           break;\
881         }\
882     }\
883 \
884   if (e > 4)\
885     e = 4;\
886   if (e < 0)\
887     nb = 4;\
888   else\
889     nb = e + 2;\
890 \
891   nb = nb >= w ? 0 : nb;\
892   newf->format = FMT_F;\
893   newf->u.real.w = f->u.real.w - nb;\
894 \
895   if (m == 0.0)\
896     newf->u.real.d = d - 1;\
897   else\
898     newf->u.real.d = - (mid - d - 1);\
899 \
900   dtp->u.p.scale_factor = 0;\
901 \
902  finish:\
903   result = output_float (dtp, newf, buffer, size, sign_bit, zero_flag, \
904                          ndigits, edigits);\
905   dtp->u.p.scale_factor = save_scale_factor;\
906 \
907   free (newf);\
908 \
909   if (nb > 0 && !dtp->u.p.g0_no_blanks)\
910     {\
911       p = write_block (dtp, nb);\
912       if (p == NULL)\
913         return;\
914       if (result == FAILURE)\
915         pad = '*';\
916       if (unlikely (is_char4_unit (dtp)))\
917         {\
918           gfc_char4_t *p4 = (gfc_char4_t *) p;\
919           memset4 (p4, pad, nb);\
920         }\
921       else\
922         memset (p, pad, nb);\
923     }\
924 }\
925
926 OUTPUT_FLOAT_FMT_G(4)
927
928 OUTPUT_FLOAT_FMT_G(8)
929
930 #ifdef HAVE_GFC_REAL_10
931 OUTPUT_FLOAT_FMT_G(10)
932 #endif
933
934 #ifdef HAVE_GFC_REAL_16
935 OUTPUT_FLOAT_FMT_G(16)
936 #endif
937
938 #undef OUTPUT_FLOAT_FMT_G
939
940
941 /* Define a macro to build code for write_float.  */
942
943   /* Note: Before output_float is called, sprintf is used to print to buffer the
944      number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
945      (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
946      before the decimal point.
947
948      #   The result will always contain a decimal point, even if no
949          digits follow it
950
951      -   The converted value is to be left adjusted on the field boundary
952
953      +   A sign (+ or -) always be placed before a number
954
955      MIN_FIELD_WIDTH  minimum field width
956
957      *   (ndigits-1) is used as the precision
958
959      e format: [-]d.ddde±dd where there is one digit before the
960        decimal-point character and the number of digits after it is
961        equal to the precision. The exponent always contains at least two
962        digits; if the value is zero, the exponent is 00.  */
963
964 #ifdef HAVE_SNPRINTF
965
966 #define DTOA \
967 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
968           "e", ndigits - 1, tmp);
969
970 #define DTOAL \
971 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
972           "Le", ndigits - 1, tmp);
973
974 #else
975
976 #define DTOA \
977 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
978          "e", ndigits - 1, tmp);
979
980 #define DTOAL \
981 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
982          "Le", ndigits - 1, tmp);
983
984 #endif
985
986 #if defined(GFC_REAL_16_IS_FLOAT128)
987 #define DTOAQ \
988 __qmath_(quadmath_snprintf) (buffer, sizeof buffer, \
989                              "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
990                              "Qe", ndigits - 1, tmp);
991 #endif
992
993 #define WRITE_FLOAT(x,y)\
994 {\
995         GFC_REAL_ ## x tmp;\
996         tmp = * (GFC_REAL_ ## x *)source;\
997         sign_bit = signbit (tmp);\
998         if (!isfinite (tmp))\
999           { \
1000             write_infnan (dtp, f, isnan (tmp), sign_bit);\
1001             return;\
1002           }\
1003         tmp = sign_bit ? -tmp : tmp;\
1004         zero_flag = (tmp == 0.0);\
1005 \
1006         DTOA ## y\
1007 \
1008         if (f->format != FMT_G)\
1009           output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
1010                         edigits);\
1011         else \
1012           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1013                                     zero_flag, ndigits, edigits);\
1014 }\
1015
1016 /* Output a real number according to its format.  */
1017
1018 static void
1019 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
1020 {
1021
1022 #if defined(HAVE_GFC_REAL_16) || __LDBL_DIG__ > 18
1023 # define MIN_FIELD_WIDTH 46
1024 #else
1025 # define MIN_FIELD_WIDTH 31
1026 #endif
1027 #define STR(x) STR1(x)
1028 #define STR1(x) #x
1029
1030   /* This must be large enough to accurately hold any value.  */
1031   char buffer[MIN_FIELD_WIDTH+1];
1032   int sign_bit, ndigits, edigits;
1033   bool zero_flag;
1034   size_t size;
1035
1036   size = MIN_FIELD_WIDTH+1;
1037
1038   /* printf pads blanks for us on the exponent so we just need it big enough
1039      to handle the largest number of exponent digits expected.  */
1040   edigits=4;
1041
1042   if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G 
1043       || ((f->format == FMT_D || f->format == FMT_E)
1044       && dtp->u.p.scale_factor != 0))
1045     {
1046       /* Always convert at full precision to avoid double rounding.  */
1047       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1048     }
1049   else
1050     {
1051       /* The number of digits is known, so let printf do the rounding.  */
1052       if (f->format == FMT_ES)
1053         ndigits = f->u.real.d + 1;
1054       else
1055         ndigits = f->u.real.d;
1056       if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
1057         ndigits = MIN_FIELD_WIDTH - 4 - edigits;
1058     }
1059
1060   switch (len)
1061     {
1062     case 4:
1063       WRITE_FLOAT(4,)
1064       break;
1065
1066     case 8:
1067       WRITE_FLOAT(8,)
1068       break;
1069
1070 #ifdef HAVE_GFC_REAL_10
1071     case 10:
1072       WRITE_FLOAT(10,L)
1073       break;
1074 #endif
1075 #ifdef HAVE_GFC_REAL_16
1076     case 16:
1077 # ifdef GFC_REAL_16_IS_FLOAT128
1078       WRITE_FLOAT(16,Q)
1079 # else
1080       WRITE_FLOAT(16,L)
1081 # endif
1082       break;
1083 #endif
1084     default:
1085       internal_error (NULL, "bad real kind");
1086     }
1087 }