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