OSDN Git Service

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