OSDN Git Service

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