xref: /netbsd-src/external/gpl3/gcc/dist/libgfortran/io/write_float.def (revision b1e838363e3c6fc78a55519254d99869742dd33c) 
1/* Copyright (C) 2007-2022 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
803#ifdef HAVE_GFC_REAL_17
804CALCULATE_EXP(17)
805#endif
806#undef CALCULATE_EXP
807
808
809/* Define macros to build code for format_float.  */
810
811  /* Note: Before output_float is called, snprintf is used to print to buffer the
812     number in the format +D.DDDDe+ddd.
813
814     #   The result will always contain a decimal point, even if no
815	 digits follow it
816
817     -   The converted value is to be left adjusted on the field boundary
818
819     +   A sign (+ or -) always be placed before a number
820
821     *   prec is used as the precision
822
823     e format: [-]d.ddde±dd where there is one digit before the
824       decimal-point character and the number of digits after it is
825       equal to the precision. The exponent always contains at least two
826       digits; if the value is zero, the exponent is 00.  */
827
828
829#define TOKENPASTE(x, y) TOKENPASTE2(x, y)
830#define TOKENPASTE2(x, y) x ## y
831
832#define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
833
834#define DTOA2(prec,val) \
835snprintf (buffer, size, "%+-#.*e", (prec), (val))
836
837#define DTOA2L(prec,val) \
838snprintf (buffer, size, "%+-#.*Le", (prec), (val))
839
840
841#if defined(HAVE_GFC_REAL_17)
842# if defined(POWER_IEEE128)
843#  define DTOA2Q(prec,val) \
844__snprintfieee128 (buffer, size, "%+-#.*Le", (prec), (val))
845# else
846#  define DTOA2Q(prec,val) \
847quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val))
848# endif
849#elif defined(GFC_REAL_16_IS_FLOAT128)
850# define DTOA2Q(prec,val) \
851quadmath_snprintf (buffer, size, "%+-#.*Qe", (prec), (val))
852#endif
853
854#define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
855
856/* For F format, we print to the buffer with f format.  */
857#define FDTOA2(prec,val) \
858snprintf (buffer, size, "%+-#.*f", (prec), (val))
859
860#define FDTOA2L(prec,val) \
861snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
862
863
864#if defined(HAVE_GFC_REAL_17)
865# if defined(POWER_IEEE128)
866#  define FDTOA2Q(prec,val) \
867__snprintfieee128 (buffer, size, "%+-#.*Lf", (prec), (val))
868# else
869# define FDTOA2Q(prec,val) \
870quadmath_snprintf (buffer, size, "%+-#.*Qf", (prec), (val))
871# endif
872#elif defined(GFC_REAL_16_IS_FLOAT128)
873# define FDTOA2Q(prec,val) \
874quadmath_snprintf (buffer, size, "%+-#.*Qf", (prec), (val))
875#endif
876
877
878/* EN format is tricky since the number of significant digits depends
879   on the magnitude.  Solve it by first printing a temporary value and
880   figure out the number of significant digits from the printed
881   exponent.  Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
882   10.0**e even when the final result will not be rounded to 10.0**e.
883   For these values the exponent returned by atoi has to be decremented
884   by one. The values y in the ranges
885       (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))
886        (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
887         (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
888   are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
889   100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
890   represents d zeroes, by the lines 279 to 297. */
891#define EN_PREC(x,y)\
892{\
893    volatile GFC_REAL_ ## x tmp, one = 1.0;\
894    tmp = * (GFC_REAL_ ## x *)source;\
895    if (isfinite (tmp))\
896      {\
897	nprinted = DTOA(y,0,tmp);\
898	int e = atoi (&buffer[4]);\
899	if (buffer[1] == '1')\
900	  {\
901	    tmp = (calculate_exp_ ## x (-e)) * tmp;\
902	    tmp = one - (tmp < 0 ? -tmp : tmp);\
903	    if (tmp > 0)\
904	      e = e - 1;\
905	  }\
906	nbefore = e%3;\
907	if (nbefore < 0)\
908	  nbefore = 3 + nbefore;\
909      }\
910    else\
911      nprinted = -1;\
912}\
913
914static int
915determine_en_precision (st_parameter_dt *dtp, const fnode *f,
916			const char *source, int len)
917{
918  int nprinted;
919  char buffer[10];
920  const size_t size = 10;
921  int nbefore; /* digits before decimal point - 1.  */
922
923  switch (len)
924    {
925    case 4:
926      EN_PREC(4,)
927      break;
928
929    case 8:
930      EN_PREC(8,)
931      break;
932
933#ifdef HAVE_GFC_REAL_10
934    case 10:
935      EN_PREC(10,L)
936      break;
937#endif
938#ifdef HAVE_GFC_REAL_16
939    case 16:
940# ifdef GFC_REAL_16_IS_FLOAT128
941      EN_PREC(16,Q)
942# else
943      EN_PREC(16,L)
944# endif
945      break;
946#endif
947#ifdef HAVE_GFC_REAL_17
948    case 17:
949      EN_PREC(17,Q)
950#endif
951      break;
952    default:
953      internal_error (NULL, "bad real kind");
954    }
955
956  if (nprinted == -1)
957    return -1;
958
959  int prec = f->u.real.d + nbefore;
960  if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
961      && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
962    prec += 2 * len + 4;
963  return prec;
964}
965
966
967/* Generate corresponding I/O format. and output.
968   The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
969   LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
970
971   Data Magnitude                              Equivalent Conversion
972   0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
973   m = 0                                       F(w-n).(d-1), n' '
974   0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
975   1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
976   10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
977   ................                           ..........
978   10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
979   m >= 10**d-0.5                              Ew.d[Ee]
980
981   notes: for Gw.d ,  n' ' means 4 blanks
982	  for Gw.dEe, n' ' means e+2 blanks
983	  for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
984	  the asm volatile is required for 32-bit x86 platforms.  */
985#define FORMAT_FLOAT(x,y)\
986{\
987  int npad = 0;\
988  GFC_REAL_ ## x m;\
989  m = * (GFC_REAL_ ## x *)source;\
990  sign_bit = signbit (m);\
991  if (!isfinite (m))\
992    { \
993      build_infnan_string (dtp, f, isnan (m), sign_bit, result, res_len);\
994      return;\
995    }\
996  m = sign_bit ? -m : m;\
997  zero_flag = (m == 0.0);\
998  if (f->format == FMT_G)\
999    {\
1000      int e = f->u.real.e;\
1001      int d = f->u.real.d;\
1002      int w = f->u.real.w;\
1003      fnode newf;\
1004      GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
1005      int low, high, mid;\
1006      int ubound, lbound;\
1007      int save_scale_factor;\
1008      volatile GFC_REAL_ ## x temp;\
1009      save_scale_factor = dtp->u.p.scale_factor;\
1010      if (w == DEFAULT_WIDTH)\
1011	{\
1012	  w = default_width;\
1013	  d = precision;\
1014	}\
1015      /* The switch between FMT_E and FMT_F is based on the absolute value.  \
1016         Set r=0 for rounding toward zero and r = 1 otherwise.  \
1017	 If (exp_d - m) == 1 there is no rounding needed.  */\
1018      switch (dtp->u.p.current_unit->round_status)\
1019	{\
1020	  case ROUND_ZERO:\
1021	    r = 0.0;\
1022	    break;\
1023	  case ROUND_UP:\
1024	    r = sign_bit ? 0.0 : 1.0;\
1025	    break;\
1026	  case ROUND_DOWN:\
1027	    r = sign_bit ? 1.0 : 0.0;\
1028	    break;\
1029	  default:\
1030	    break;\
1031	}\
1032      exp_d = calculate_exp_ ## x (d);\
1033      r_sc = (1 - r / exp_d);\
1034      temp = 0.1 * r_sc;\
1035      if ((m > 0.0 && ((m < temp) || (r < 1 && r >= (exp_d - m))\
1036				  || (r == 1 && 1 > (exp_d - m))))\
1037	  || ((m == 0.0) && !(compile_options.allow_std\
1038			      & (GFC_STD_F2003 | GFC_STD_F2008)))\
1039	  ||  d == 0)\
1040	{ \
1041	  newf.format = FMT_E;\
1042	  newf.u.real.w = w;\
1043	  newf.u.real.d = d - comp_d;\
1044	  newf.u.real.e = e;\
1045	  npad = 0;\
1046	  precision = determine_precision (dtp, &newf, x);\
1047	  nprinted = DTOA(y,precision,m);\
1048	}\
1049      else \
1050	{\
1051	  mid = 0;\
1052	  low = 0;\
1053	  high = d + 1;\
1054	  lbound = 0;\
1055	  ubound = d + 1;\
1056	  while (low <= high)\
1057	    {\
1058	      mid = (low + high) / 2;\
1059	      temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
1060	      if (m < temp)\
1061		{ \
1062		  ubound = mid;\
1063		  if (ubound == lbound + 1)\
1064		    break;\
1065		  high = mid - 1;\
1066		}\
1067	      else if (m > temp)\
1068		{ \
1069		  lbound = mid;\
1070		  if (ubound == lbound + 1)\
1071		    { \
1072		      mid ++;\
1073		      break;\
1074		    }\
1075		  low = mid + 1;\
1076		}\
1077	      else\
1078		{\
1079		  mid++;\
1080		  break;\
1081		}\
1082	    }\
1083	  npad = e <= 0 ? 4 : e + 2;\
1084	  npad = npad >= w ? w - 1 : npad;\
1085	  npad = dtp->u.p.g0_no_blanks ? 0 : npad;\
1086	  newf.format = FMT_F;\
1087	  newf.u.real.w = w - npad;\
1088	  newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
1089	  dtp->u.p.scale_factor = 0;\
1090	  precision = determine_precision (dtp, &newf, x);\
1091	  nprinted = FDTOA(y,precision,m);\
1092	}\
1093      build_float_string (dtp, &newf, buffer, size, nprinted, precision,\
1094				   sign_bit, zero_flag, npad, default_width,\
1095				   result, res_len);\
1096      dtp->u.p.scale_factor = save_scale_factor;\
1097    }\
1098  else\
1099    {\
1100      if (f->format == FMT_F)\
1101	nprinted = FDTOA(y,precision,m);\
1102      else\
1103	nprinted = DTOA(y,precision,m);\
1104      build_float_string (dtp, f, buffer, size, nprinted, precision,\
1105				   sign_bit, zero_flag, npad, default_width,\
1106				   result, res_len);\
1107    }\
1108}\
1109
1110/* Output a real number according to its format.  */
1111
1112
1113static void
1114get_float_string (st_parameter_dt *dtp, const fnode *f, const char *source,
1115		  int kind, int comp_d, char *buffer, int precision,
1116		  size_t size, char *result, size_t *res_len)
1117{
1118  int sign_bit, nprinted;
1119  bool zero_flag;
1120  int default_width = 0;
1121
1122  if (f->u.real.w == DEFAULT_WIDTH)
1123    /* This codepath can only be reached with -fdec-format-defaults. The default
1124     * values are based on those used in the Oracle Fortran compiler.
1125     */
1126    {
1127      default_width = default_width_for_float (kind);
1128      precision = default_precision_for_float (kind);
1129    }
1130
1131  switch (kind)
1132    {
1133    case 4:
1134      FORMAT_FLOAT(4,)
1135      break;
1136
1137    case 8:
1138      FORMAT_FLOAT(8,)
1139      break;
1140
1141#ifdef HAVE_GFC_REAL_10
1142    case 10:
1143      FORMAT_FLOAT(10,L)
1144      break;
1145#endif
1146#ifdef HAVE_GFC_REAL_16
1147    case 16:
1148# ifdef GFC_REAL_16_IS_FLOAT128
1149      FORMAT_FLOAT(16,Q)
1150# else
1151      FORMAT_FLOAT(16,L)
1152# endif
1153      break;
1154#endif
1155#ifdef HAVE_GFC_REAL_17
1156    case 17:
1157      FORMAT_FLOAT(17,Q)
1158      break;
1159#endif
1160    default:
1161      internal_error (NULL, "bad real kind");
1162    }
1163  return;
1164}
1165