* config/rs6000/t-spe (MULTILIB_EXCEPTIONS): Allow isel without SPE.
[official-gcc.git] / libgfortran / io / write_float.def
blob99f6ff8e83f851f4feb97d4c2175e816b9f29624
1 /* Copyright (C) 2007-2014 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
6 This file is part of the GNU Fortran runtime library (libgfortran).
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.
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.
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.
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/>. */
27 #include "config.h"
29 typedef enum
30 { S_NONE, S_MINUS, S_PLUS }
31 sign_t;
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. */
36 static sign_t
37 calculate_sign (st_parameter_dt *dtp, int negative_flag)
39 sign_t s = S_NONE;
41 if (negative_flag)
42 s = S_MINUS;
43 else
44 switch (dtp->u.p.sign_status)
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;
58 return s;
62 /* Determine the precision except for EN format. For G format,
63 determines an upper bound to be used for sizing the buffer. */
65 static int
66 determine_precision (st_parameter_dt * dtp, const fnode * f, int len)
68 int precision = f->u.real.d;
70 switch (f->format)
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;
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;
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)
102 precision += 2 * len + 4;
103 if (precision < 0)
104 precision = 0;
107 return precision;
111 /* Output a real number according to its format which is FMT_G free. */
113 static bool
114 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
115 int nprinted, int precision, int sign_bit, bool zero_flag)
117 char *out;
118 char *digits;
119 int e, w, d, p, i;
120 char expchar, rchar;
121 format_token ft;
122 /* Number of digits before the decimal point. */
123 int nbefore;
124 /* Number of zeros after the decimal point. */
125 int nzero;
126 /* Number of digits after the decimal point. */
127 int nafter;
128 int leadzero;
129 int nblanks;
130 int ndigits, edigits;
131 sign_t sign;
133 ft = f->format;
134 w = f->u.real.w;
135 d = f->u.real.d;
136 p = dtp->u.p.scale_factor;
138 rchar = '5';
140 /* We should always know the field width and precision. */
141 if (d < 0)
142 internal_error (&dtp->common, "Unspecified precision");
144 sign = calculate_sign (dtp, sign_bit);
146 /* Calculate total number of digits. */
147 if (ft == FMT_F)
148 ndigits = nprinted - 2;
149 else
150 ndigits = precision + 1;
152 /* Read the exponent back in. */
153 if (ft != FMT_F)
154 e = atoi (&buffer[ndigits + 3]) + 1;
155 else
156 e = 0;
158 /* Make sure zero comes out as 0.0e0. */
159 if (zero_flag)
160 e = 0;
162 /* Normalize the fractional component. */
163 if (ft != FMT_F)
165 buffer[2] = buffer[1];
166 digits = &buffer[2];
168 else
169 digits = &buffer[1];
171 /* Figure out where to place the decimal point. */
172 switch (ft)
174 case FMT_F:
175 nbefore = ndigits - precision;
176 /* Make sure the decimal point is a '.'; depending on the
177 locale, this might not be the case otherwise. */
178 digits[nbefore] = '.';
179 if (p != 0)
181 if (p > 0)
184 memmove (digits + nbefore, digits + nbefore + 1, p);
185 digits[nbefore + p] = '.';
186 nbefore += p;
187 nafter = d - p;
188 if (nafter < 0)
189 nafter = 0;
190 nafter = d;
191 nzero = 0;
193 else /* p < 0 */
195 if (nbefore + p >= 0)
197 nzero = 0;
198 memmove (digits + nbefore + p + 1, digits + nbefore + p, -p);
199 nbefore += p;
200 digits[nbefore] = '.';
201 nafter = d;
203 else
205 nzero = -(nbefore + p);
206 memmove (digits + 1, digits, nbefore);
207 digits++;
208 nafter = d + nbefore;
209 nbefore = 0;
211 if (nzero > d)
212 nzero = d;
215 else
217 nzero = 0;
218 nafter = d;
221 while (digits[0] == '0' && nbefore > 0)
223 digits++;
224 nbefore--;
225 ndigits--;
228 expchar = 0;
229 /* If we need to do rounding ourselves, get rid of the dot by
230 moving the fractional part. */
231 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
232 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
233 memmove (digits + nbefore, digits + nbefore + 1, ndigits - nbefore);
234 break;
236 case FMT_E:
237 case FMT_D:
238 i = dtp->u.p.scale_factor;
239 if (d <= 0 && p == 0)
241 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
242 "greater than zero in format specifier 'E' or 'D'");
243 return false;
245 if (p <= -d || p >= d + 2)
247 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
248 "out of range in format specifier 'E' or 'D'");
249 return false;
252 if (!zero_flag)
253 e -= p;
254 if (p < 0)
256 nbefore = 0;
257 nzero = -p;
258 nafter = d + p;
260 else if (p > 0)
262 nbefore = p;
263 nzero = 0;
264 nafter = (d - p) + 1;
266 else /* p == 0 */
268 nbefore = 0;
269 nzero = 0;
270 nafter = d;
273 if (ft == FMT_E)
274 expchar = 'E';
275 else
276 expchar = 'D';
277 break;
279 case FMT_EN:
280 /* The exponent must be a multiple of three, with 1-3 digits before
281 the decimal point. */
282 if (!zero_flag)
283 e--;
284 if (e >= 0)
285 nbefore = e % 3;
286 else
288 nbefore = (-e) % 3;
289 if (nbefore != 0)
290 nbefore = 3 - nbefore;
292 e -= nbefore;
293 nbefore++;
294 nzero = 0;
295 nafter = d;
296 expchar = 'E';
297 break;
299 case FMT_ES:
300 if (!zero_flag)
301 e--;
302 nbefore = 1;
303 nzero = 0;
304 nafter = d;
305 expchar = 'E';
306 break;
308 default:
309 /* Should never happen. */
310 internal_error (&dtp->common, "Unexpected format token");
313 if (zero_flag)
314 goto skip;
316 /* Round the value. The value being rounded is an unsigned magnitude. */
317 switch (dtp->u.p.current_unit->round_status)
319 /* For processor defined and unspecified rounding we use
320 snprintf to print the exact number of digits needed, and thus
321 let snprintf handle the rounding. On system claiming support
322 for IEEE 754, this ought to be round to nearest, ties to
323 even, corresponding to the Fortran ROUND='NEAREST'. */
324 case ROUND_PROCDEFINED:
325 case ROUND_UNSPECIFIED:
326 case ROUND_ZERO: /* Do nothing and truncation occurs. */
327 goto skip;
328 case ROUND_UP:
329 if (sign_bit)
330 goto skip;
331 goto updown;
332 case ROUND_DOWN:
333 if (!sign_bit)
334 goto skip;
335 goto updown;
336 case ROUND_NEAREST:
337 /* Round compatible unless there is a tie. A tie is a 5 with
338 all trailing zero's. */
339 i = nafter + nbefore;
340 if (digits[i] == '5')
342 for(i++ ; i < ndigits; i++)
344 if (digits[i] != '0')
345 goto do_rnd;
347 /* It is a tie so round to even. */
348 switch (digits[nafter + nbefore - 1])
350 case '1':
351 case '3':
352 case '5':
353 case '7':
354 case '9':
355 /* If odd, round away from zero to even. */
356 break;
357 default:
358 /* If even, skip rounding, truncate to even. */
359 goto skip;
362 /* Fall through. */
363 /* The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */
364 case ROUND_COMPATIBLE:
365 rchar = '5';
366 goto do_rnd;
369 updown:
371 rchar = '0';
372 if (ft != FMT_F && w > 0 && d == 0 && p == 0)
373 nbefore = 1;
374 /* Scan for trailing zeros to see if we really need to round it. */
375 for(i = nbefore + nafter; i < ndigits; i++)
377 if (digits[i] != '0')
378 goto do_rnd;
380 goto skip;
382 do_rnd:
384 if (nbefore + nafter == 0)
385 /* Handle the case Fw.0 and value < 1.0 */
387 ndigits = 0;
388 if (digits[0] >= rchar)
390 /* We rounded to zero but shouldn't have */
391 nbefore = 1;
392 digits--;
393 digits[0] = '1';
394 ndigits = 1;
397 else if (nbefore + nafter < ndigits)
399 i = ndigits = nbefore + nafter;
400 if (digits[i] >= rchar)
402 /* Propagate the carry. */
403 for (i--; i >= 0; i--)
405 if (digits[i] != '9')
407 digits[i]++;
408 break;
410 digits[i] = '0';
413 if (i < 0)
415 /* The carry overflowed. Fortunately we have some spare
416 space at the start of the buffer. We may discard some
417 digits, but this is ok because we already know they are
418 zero. */
419 digits--;
420 digits[0] = '1';
421 if (ft == FMT_F)
423 if (nzero > 0)
425 nzero--;
426 nafter++;
428 else
429 nbefore++;
431 else if (ft == FMT_EN)
433 nbefore++;
434 if (nbefore == 4)
436 nbefore = 1;
437 e += 3;
440 else
441 e++;
446 skip:
448 /* Calculate the format of the exponent field. */
449 if (expchar)
451 edigits = 1;
452 for (i = abs (e); i >= 10; i /= 10)
453 edigits++;
455 if (f->u.real.e < 0)
457 /* Width not specified. Must be no more than 3 digits. */
458 if (e > 999 || e < -999)
459 edigits = -1;
460 else
462 edigits = 4;
463 if (e > 99 || e < -99)
464 expchar = ' ';
467 else
469 /* Exponent width specified, check it is wide enough. */
470 if (edigits > f->u.real.e)
471 edigits = -1;
472 else
473 edigits = f->u.real.e + 2;
476 else
477 edigits = 0;
479 /* Scan the digits string and count the number of zeros. If we make it
480 all the way through the loop, we know the value is zero after the
481 rounding completed above. */
482 int hasdot = 0;
483 for (i = 0; i < ndigits + hasdot; i++)
485 if (digits[i] == '.')
486 hasdot = 1;
487 else if (digits[i] != '0')
488 break;
491 /* To format properly, we need to know if the rounded result is zero and if
492 so, we set the zero_flag which may have been already set for
493 actual zero. */
494 if (i == ndigits + hasdot)
496 zero_flag = true;
497 /* The output is zero, so set the sign according to the sign bit unless
498 -fno-sign-zero was specified. */
499 if (compile_options.sign_zero == 1)
500 sign = calculate_sign (dtp, sign_bit);
501 else
502 sign = calculate_sign (dtp, 0);
505 /* Pick a field size if none was specified, taking into account small
506 values that may have been rounded to zero. */
507 if (w <= 0)
509 if (zero_flag)
510 w = d + (sign != S_NONE ? 2 : 1) + (d == 0 ? 1 : 0);
511 else
513 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
514 w = w == 1 ? 2 : w;
518 /* Work out how much padding is needed. */
519 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
520 if (sign != S_NONE)
521 nblanks--;
523 if (dtp->u.p.g0_no_blanks)
525 w -= nblanks;
526 nblanks = 0;
529 /* Create the ouput buffer. */
530 out = write_block (dtp, w);
531 if (out == NULL)
532 return false;
534 /* Check the value fits in the specified field width. */
535 if (nblanks < 0 || edigits == -1 || w == 1 || (w == 2 && sign != S_NONE))
537 if (unlikely (is_char4_unit (dtp)))
539 gfc_char4_t *out4 = (gfc_char4_t *) out;
540 memset4 (out4, '*', w);
541 return false;
543 star_fill (out, w);
544 return false;
547 /* See if we have space for a zero before the decimal point. */
548 if (nbefore == 0 && nblanks > 0)
550 leadzero = 1;
551 nblanks--;
553 else
554 leadzero = 0;
556 /* For internal character(kind=4) units, we duplicate the code used for
557 regular output slightly modified. This needs to be maintained
558 consistent with the regular code that follows this block. */
559 if (unlikely (is_char4_unit (dtp)))
561 gfc_char4_t *out4 = (gfc_char4_t *) out;
562 /* Pad to full field width. */
564 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
566 memset4 (out4, ' ', nblanks);
567 out4 += nblanks;
570 /* Output the initial sign (if any). */
571 if (sign == S_PLUS)
572 *(out4++) = '+';
573 else if (sign == S_MINUS)
574 *(out4++) = '-';
576 /* Output an optional leading zero. */
577 if (leadzero)
578 *(out4++) = '0';
580 /* Output the part before the decimal point, padding with zeros. */
581 if (nbefore > 0)
583 if (nbefore > ndigits)
585 i = ndigits;
586 memcpy4 (out4, digits, i);
587 ndigits = 0;
588 while (i < nbefore)
589 out4[i++] = '0';
591 else
593 i = nbefore;
594 memcpy4 (out4, digits, i);
595 ndigits -= i;
598 digits += i;
599 out4 += nbefore;
602 /* Output the decimal point. */
603 *(out4++) = dtp->u.p.current_unit->decimal_status
604 == DECIMAL_POINT ? '.' : ',';
605 if (ft == FMT_F
606 && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
607 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
608 digits++;
610 /* Output leading zeros after the decimal point. */
611 if (nzero > 0)
613 for (i = 0; i < nzero; i++)
614 *(out4++) = '0';
617 /* Output digits after the decimal point, padding with zeros. */
618 if (nafter > 0)
620 if (nafter > ndigits)
621 i = ndigits;
622 else
623 i = nafter;
625 memcpy4 (out4, digits, i);
626 while (i < nafter)
627 out4[i++] = '0';
629 digits += i;
630 ndigits -= i;
631 out4 += nafter;
634 /* Output the exponent. */
635 if (expchar)
637 if (expchar != ' ')
639 *(out4++) = expchar;
640 edigits--;
642 snprintf (buffer, size, "%+0*d", edigits, e);
643 memcpy4 (out4, buffer, edigits);
646 if (dtp->u.p.no_leading_blank)
648 out4 += edigits;
649 memset4 (out4, ' ' , nblanks);
650 dtp->u.p.no_leading_blank = 0;
652 return true;
653 } /* End of character(kind=4) internal unit code. */
655 /* Pad to full field width. */
657 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
659 memset (out, ' ', nblanks);
660 out += nblanks;
663 /* Output the initial sign (if any). */
664 if (sign == S_PLUS)
665 *(out++) = '+';
666 else if (sign == S_MINUS)
667 *(out++) = '-';
669 /* Output an optional leading zero. */
670 if (leadzero)
671 *(out++) = '0';
673 /* Output the part before the decimal point, padding with zeros. */
674 if (nbefore > 0)
676 if (nbefore > ndigits)
678 i = ndigits;
679 memcpy (out, digits, i);
680 ndigits = 0;
681 while (i < nbefore)
682 out[i++] = '0';
684 else
686 i = nbefore;
687 memcpy (out, digits, i);
688 ndigits -= i;
691 digits += i;
692 out += nbefore;
695 /* Output the decimal point. */
696 *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
697 if (ft == FMT_F
698 && (dtp->u.p.current_unit->round_status == ROUND_UNSPECIFIED
699 || dtp->u.p.current_unit->round_status == ROUND_PROCDEFINED))
700 digits++;
702 /* Output leading zeros after the decimal point. */
703 if (nzero > 0)
705 for (i = 0; i < nzero; i++)
706 *(out++) = '0';
709 /* Output digits after the decimal point, padding with zeros. */
710 if (nafter > 0)
712 if (nafter > ndigits)
713 i = ndigits;
714 else
715 i = nafter;
717 memcpy (out, digits, i);
718 while (i < nafter)
719 out[i++] = '0';
721 digits += i;
722 ndigits -= i;
723 out += nafter;
726 /* Output the exponent. */
727 if (expchar)
729 if (expchar != ' ')
731 *(out++) = expchar;
732 edigits--;
734 snprintf (buffer, size, "%+0*d", edigits, e);
735 memcpy (out, buffer, edigits);
738 if (dtp->u.p.no_leading_blank)
740 out += edigits;
741 memset( out , ' ' , nblanks );
742 dtp->u.p.no_leading_blank = 0;
745 return true;
749 /* Write "Infinite" or "Nan" as appropriate for the given format. */
751 static void
752 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
754 char * p, fin;
755 int nb = 0;
756 sign_t sign;
757 int mark;
759 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
761 sign = calculate_sign (dtp, sign_bit);
762 mark = (sign == S_PLUS || sign == S_MINUS) ? 8 : 7;
764 nb = f->u.real.w;
766 /* If the field width is zero, the processor must select a width
767 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
769 if ((nb == 0) || dtp->u.p.g0_no_blanks)
771 if (isnan_flag)
772 nb = 3;
773 else
774 nb = (sign == S_PLUS || sign == S_MINUS) ? 4 : 3;
776 p = write_block (dtp, nb);
777 if (p == NULL)
778 return;
779 if (nb < 3)
781 if (unlikely (is_char4_unit (dtp)))
783 gfc_char4_t *p4 = (gfc_char4_t *) p;
784 memset4 (p4, '*', nb);
786 else
787 memset (p, '*', nb);
788 return;
791 if (unlikely (is_char4_unit (dtp)))
793 gfc_char4_t *p4 = (gfc_char4_t *) p;
794 memset4 (p4, ' ', nb);
796 else
797 memset(p, ' ', nb);
799 if (!isnan_flag)
801 if (sign_bit)
803 /* If the sign is negative and the width is 3, there is
804 insufficient room to output '-Inf', so output asterisks */
805 if (nb == 3)
807 if (unlikely (is_char4_unit (dtp)))
809 gfc_char4_t *p4 = (gfc_char4_t *) p;
810 memset4 (p4, '*', nb);
812 else
813 memset (p, '*', nb);
814 return;
816 /* The negative sign is mandatory */
817 fin = '-';
819 else
820 /* The positive sign is optional, but we output it for
821 consistency */
822 fin = '+';
824 if (unlikely (is_char4_unit (dtp)))
826 gfc_char4_t *p4 = (gfc_char4_t *) p;
828 if (nb > mark)
829 /* We have room, so output 'Infinity' */
830 memcpy4 (p4 + nb - 8, "Infinity", 8);
831 else
832 /* For the case of width equals mark, there is not enough room
833 for the sign and 'Infinity' so we go with 'Inf' */
834 memcpy4 (p4 + nb - 3, "Inf", 3);
836 if (sign == S_PLUS || sign == S_MINUS)
838 if (nb < 9 && nb > 3)
839 /* Put the sign in front of Inf */
840 p4[nb - 4] = (gfc_char4_t) fin;
841 else if (nb > 8)
842 /* Put the sign in front of Infinity */
843 p4[nb - 9] = (gfc_char4_t) fin;
845 return;
848 if (nb > mark)
849 /* We have room, so output 'Infinity' */
850 memcpy(p + nb - 8, "Infinity", 8);
851 else
852 /* For the case of width equals 8, there is not enough room
853 for the sign and 'Infinity' so we go with 'Inf' */
854 memcpy(p + nb - 3, "Inf", 3);
856 if (sign == S_PLUS || sign == S_MINUS)
858 if (nb < 9 && nb > 3)
859 p[nb - 4] = fin; /* Put the sign in front of Inf */
860 else if (nb > 8)
861 p[nb - 9] = fin; /* Put the sign in front of Infinity */
864 else
866 if (unlikely (is_char4_unit (dtp)))
868 gfc_char4_t *p4 = (gfc_char4_t *) p;
869 memcpy4 (p4 + nb - 3, "NaN", 3);
871 else
872 memcpy(p + nb - 3, "NaN", 3);
874 return;
879 /* Returns the value of 10**d. */
881 #define CALCULATE_EXP(x) \
882 static GFC_REAL_ ## x \
883 calculate_exp_ ## x (int d)\
885 int i;\
886 GFC_REAL_ ## x r = 1.0;\
887 for (i = 0; i< (d >= 0 ? d : -d); i++)\
888 r *= 10;\
889 r = (d >= 0) ? r : 1.0 / r;\
890 return r;\
893 CALCULATE_EXP(4)
895 CALCULATE_EXP(8)
897 #ifdef HAVE_GFC_REAL_10
898 CALCULATE_EXP(10)
899 #endif
901 #ifdef HAVE_GFC_REAL_16
902 CALCULATE_EXP(16)
903 #endif
904 #undef CALCULATE_EXP
907 /* Define a macro to build code for write_float. */
909 /* Note: Before output_float is called, snprintf is used to print to buffer the
910 number in the format +D.DDDDe+ddd.
912 # The result will always contain a decimal point, even if no
913 digits follow it
915 - The converted value is to be left adjusted on the field boundary
917 + A sign (+ or -) always be placed before a number
919 * prec is used as the precision
921 e format: [-]d.ddde±dd where there is one digit before the
922 decimal-point character and the number of digits after it is
923 equal to the precision. The exponent always contains at least two
924 digits; if the value is zero, the exponent is 00. */
927 #define TOKENPASTE(x, y) TOKENPASTE2(x, y)
928 #define TOKENPASTE2(x, y) x ## y
930 #define DTOA(suff,prec,val) TOKENPASTE(DTOA2,suff)(prec,val)
932 #define DTOA2(prec,val) \
933 snprintf (buffer, size, "%+-#.*e", (prec), (val))
935 #define DTOA2L(prec,val) \
936 snprintf (buffer, size, "%+-#.*Le", (prec), (val))
939 #if defined(GFC_REAL_16_IS_FLOAT128)
940 #define DTOA2Q(prec,val) \
941 __qmath_(quadmath_snprintf) (buffer, size, "%+-#.*Qe", (prec), (val))
942 #endif
944 #define FDTOA(suff,prec,val) TOKENPASTE(FDTOA2,suff)(prec,val)
946 /* For F format, we print to the buffer with f format. */
947 #define FDTOA2(prec,val) \
948 snprintf (buffer, size, "%+-#.*f", (prec), (val))
950 #define FDTOA2L(prec,val) \
951 snprintf (buffer, size, "%+-#.*Lf", (prec), (val))
954 #if defined(GFC_REAL_16_IS_FLOAT128)
955 #define FDTOA2Q(prec,val) \
956 __qmath_(quadmath_snprintf) (buffer, size, "%+-#.*Qf", \
957 (prec), (val))
958 #endif
961 #if defined(GFC_REAL_16_IS_FLOAT128)
962 #define ISFINITE2Q(val) finiteq(val)
963 #endif
964 #define ISFINITE2(val) isfinite(val)
965 #define ISFINITE2L(val) isfinite(val)
967 #define ISFINITE(suff,val) TOKENPASTE(ISFINITE2,suff)(val)
970 #if defined(GFC_REAL_16_IS_FLOAT128)
971 #define SIGNBIT2Q(val) signbitq(val)
972 #endif
973 #define SIGNBIT2(val) signbit(val)
974 #define SIGNBIT2L(val) signbit(val)
976 #define SIGNBIT(suff,val) TOKENPASTE(SIGNBIT2,suff)(val)
979 #if defined(GFC_REAL_16_IS_FLOAT128)
980 #define ISNAN2Q(val) isnanq(val)
981 #endif
982 #define ISNAN2(val) isnan(val)
983 #define ISNAN2L(val) isnan(val)
985 #define ISNAN(suff,val) TOKENPASTE(ISNAN2,suff)(val)
989 /* Generate corresponding I/O format for FMT_G and output.
990 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
991 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
993 Data Magnitude Equivalent Conversion
994 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
995 m = 0 F(w-n).(d-1), n' '
996 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
997 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
998 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
999 ................ ..........
1000 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
1001 m >= 10**d-0.5 Ew.d[Ee]
1003 notes: for Gw.d , n' ' means 4 blanks
1004 for Gw.dEe, n' ' means e+2 blanks
1005 for rounding modes adjustment, r, See Fortran F2008 10.7.5.2.2
1006 the asm volatile is required for 32-bit x86 platforms. */
1008 #define OUTPUT_FLOAT_FMT_G(x,y) \
1009 static void \
1010 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
1011 GFC_REAL_ ## x m, char *buffer, size_t size, \
1012 int sign_bit, bool zero_flag, int comp_d) \
1014 int e = f->u.real.e;\
1015 int d = f->u.real.d;\
1016 int w = f->u.real.w;\
1017 fnode newf;\
1018 GFC_REAL_ ## x exp_d, r = 0.5, r_sc;\
1019 int low, high, mid;\
1020 int ubound, lbound;\
1021 char *p, pad = ' ';\
1022 int save_scale_factor, nb = 0;\
1023 bool result;\
1024 int nprinted, precision;\
1025 volatile GFC_REAL_ ## x temp;\
1027 save_scale_factor = dtp->u.p.scale_factor;\
1029 switch (dtp->u.p.current_unit->round_status)\
1031 case ROUND_ZERO:\
1032 r = sign_bit ? 1.0 : 0.0;\
1033 break;\
1034 case ROUND_UP:\
1035 r = 1.0;\
1036 break;\
1037 case ROUND_DOWN:\
1038 r = 0.0;\
1039 break;\
1040 default:\
1041 break;\
1044 exp_d = calculate_exp_ ## x (d);\
1045 r_sc = (1 - r / exp_d);\
1046 temp = 0.1 * r_sc;\
1047 if ((m > 0.0 && ((m < temp) || (r >= (exp_d - m))))\
1048 || ((m == 0.0) && !(compile_options.allow_std\
1049 & (GFC_STD_F2003 | GFC_STD_F2008)))\
1050 || d == 0)\
1052 newf.format = FMT_E;\
1053 newf.u.real.w = w;\
1054 newf.u.real.d = d - comp_d;\
1055 newf.u.real.e = e;\
1056 nb = 0;\
1057 precision = determine_precision (dtp, &newf, x);\
1058 nprinted = DTOA(y,precision,m); \
1059 goto finish;\
1062 mid = 0;\
1063 low = 0;\
1064 high = d + 1;\
1065 lbound = 0;\
1066 ubound = d + 1;\
1068 while (low <= high)\
1070 mid = (low + high) / 2;\
1072 temp = (calculate_exp_ ## x (mid - 1) * r_sc);\
1074 if (m < temp)\
1076 ubound = mid;\
1077 if (ubound == lbound + 1)\
1078 break;\
1079 high = mid - 1;\
1081 else if (m > temp)\
1083 lbound = mid;\
1084 if (ubound == lbound + 1)\
1086 mid ++;\
1087 break;\
1089 low = mid + 1;\
1091 else\
1093 mid++;\
1094 break;\
1098 nb = e <= 0 ? 4 : e + 2;\
1099 nb = nb >= w ? w - 1 : nb;\
1100 newf.format = FMT_F;\
1101 newf.u.real.w = w - nb;\
1102 newf.u.real.d = m == 0.0 ? d - 1 : -(mid - d - 1) ;\
1103 dtp->u.p.scale_factor = 0;\
1104 precision = determine_precision (dtp, &newf, x); \
1105 nprinted = FDTOA(y,precision,m); \
1107 finish:\
1108 result = output_float (dtp, &newf, buffer, size, nprinted, precision,\
1109 sign_bit, zero_flag);\
1110 dtp->u.p.scale_factor = save_scale_factor;\
1113 if (nb > 0 && !dtp->u.p.g0_no_blanks)\
1115 p = write_block (dtp, nb);\
1116 if (p == NULL)\
1117 return;\
1118 if (!result)\
1119 pad = '*';\
1120 if (unlikely (is_char4_unit (dtp)))\
1122 gfc_char4_t *p4 = (gfc_char4_t *) p;\
1123 memset4 (p4, pad, nb);\
1125 else \
1126 memset (p, pad, nb);\
1130 OUTPUT_FLOAT_FMT_G(4,)
1132 OUTPUT_FLOAT_FMT_G(8,)
1134 #ifdef HAVE_GFC_REAL_10
1135 OUTPUT_FLOAT_FMT_G(10,L)
1136 #endif
1138 #ifdef HAVE_GFC_REAL_16
1139 # ifdef GFC_REAL_16_IS_FLOAT128
1140 OUTPUT_FLOAT_FMT_G(16,Q)
1141 #else
1142 OUTPUT_FLOAT_FMT_G(16,L)
1143 #endif
1144 #endif
1146 #undef OUTPUT_FLOAT_FMT_G
1149 /* EN format is tricky since the number of significant digits depends
1150 on the magnitude. Solve it by first printing a temporary value and
1151 figure out the number of significant digits from the printed
1152 exponent. Values y, 0.95*10.0**e <= y <10.0**e, are rounded to
1153 10.0**e even when the final result will not be rounded to 10.0**e.
1154 For these values the exponent returned by atoi has to be decremented
1155 by one. The values y in the ranges
1156 (1000.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*(n+1))
1157 (100.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+2)
1158 (10.0-0.5*10.0**(-d))*10.0**(3*n) <= y < 10.0*(3*n+1)
1159 are correctly rounded respectively to 1.0...0*10.0*(3*(n+1)),
1160 100.0...0*10.0*(3*n), and 10.0...0*10.0*(3*n), where 0...0
1161 represents d zeroes, by the lines 279 to 297. */
1163 #define EN_PREC(x,y)\
1165 volatile GFC_REAL_ ## x tmp, one = 1.0;\
1166 tmp = * (GFC_REAL_ ## x *)source;\
1167 if (ISFINITE (y,tmp))\
1169 nprinted = DTOA(y,0,tmp);\
1170 int e = atoi (&buffer[4]);\
1171 if (buffer[1] == '1')\
1173 tmp = (calculate_exp_ ## x (-e)) * tmp;\
1174 tmp = one - (tmp < 0 ? -tmp : tmp); \
1175 if (tmp > 0)\
1176 e = e - 1;\
1178 nbefore = e%3;\
1179 if (nbefore < 0)\
1180 nbefore = 3 + nbefore;\
1182 else\
1183 nprinted = -1;\
1186 static int
1187 determine_en_precision (st_parameter_dt *dtp, const fnode *f,
1188 const char *source, int len)
1190 int nprinted;
1191 char buffer[10];
1192 const size_t size = 10;
1193 int nbefore; /* digits before decimal point - 1. */
1195 switch (len)
1197 case 4:
1198 EN_PREC(4,)
1199 break;
1201 case 8:
1202 EN_PREC(8,)
1203 break;
1205 #ifdef HAVE_GFC_REAL_10
1206 case 10:
1207 EN_PREC(10,L)
1208 break;
1209 #endif
1210 #ifdef HAVE_GFC_REAL_16
1211 case 16:
1212 # ifdef GFC_REAL_16_IS_FLOAT128
1213 EN_PREC(16,Q)
1214 # else
1215 EN_PREC(16,L)
1216 # endif
1217 break;
1218 #endif
1219 default:
1220 internal_error (NULL, "bad real kind");
1223 if (nprinted == -1)
1224 return -1;
1226 int prec = f->u.real.d + nbefore;
1227 if (dtp->u.p.current_unit->round_status != ROUND_UNSPECIFIED
1228 && dtp->u.p.current_unit->round_status != ROUND_PROCDEFINED)
1229 prec += 2 * len + 4;
1230 return prec;
1234 #define WRITE_FLOAT(x,y)\
1236 GFC_REAL_ ## x tmp;\
1237 tmp = * (GFC_REAL_ ## x *)source;\
1238 sign_bit = SIGNBIT (y,tmp);\
1239 if (!ISFINITE (y,tmp))\
1241 write_infnan (dtp, f, ISNAN (y,tmp), sign_bit);\
1242 return;\
1244 tmp = sign_bit ? -tmp : tmp;\
1245 zero_flag = (tmp == 0.0);\
1246 if (f->format == FMT_G)\
1247 output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
1248 zero_flag, comp_d);\
1249 else\
1251 if (f->format == FMT_F)\
1252 nprinted = FDTOA(y,precision,tmp); \
1253 else\
1254 nprinted = DTOA(y,precision,tmp); \
1255 output_float (dtp, f, buffer, size, nprinted, precision,\
1256 sign_bit, zero_flag);\
1260 /* Output a real number according to its format. */
1262 static void
1263 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, \
1264 int len, int comp_d)
1266 int sign_bit, nprinted;
1267 int precision; /* Precision for snprintf call. */
1268 bool zero_flag;
1270 if (f->format != FMT_EN)
1271 precision = determine_precision (dtp, f, len);
1272 else
1273 precision = determine_en_precision (dtp, f, source, len);
1275 /* 4932 is the maximum exponent of long double and quad precision, 3
1276 extra characters for the sign, the decimal point, and the
1277 trailing null, and finally some extra digits depending on the
1278 requested precision. */
1279 const size_t size = 4932 + 3 + precision;
1280 char buffer[size];
1282 switch (len)
1284 case 4:
1285 WRITE_FLOAT(4,)
1286 break;
1288 case 8:
1289 WRITE_FLOAT(8,)
1290 break;
1292 #ifdef HAVE_GFC_REAL_10
1293 case 10:
1294 WRITE_FLOAT(10,L)
1295 break;
1296 #endif
1297 #ifdef HAVE_GFC_REAL_16
1298 case 16:
1299 # ifdef GFC_REAL_16_IS_FLOAT128
1300 WRITE_FLOAT(16,Q)
1301 # else
1302 WRITE_FLOAT(16,L)
1303 # endif
1304 break;
1305 #endif
1306 default:
1307 internal_error (NULL, "bad real kind");