OSDN Git Service

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