OSDN Git Service

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