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