OSDN Git Service

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