Remove redundant checks in graphite_can_represent_scev.
[official-gcc/graphite-test-results.git] / libgfortran / io / write_float.def
blob45c2a17a50d805e7f9c7e45d2aab6d7f8e5ba1b5
1 /* Copyright (C) 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
2 Contributed by Andy Vaught
3 Write float code factoring to this file by Jerry DeLisle
4 F2003 I/O support contributed by Jerry DeLisle
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 /* Output a real number according to its format which is FMT_G free. */
64 static void
65 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
66 int sign_bit, bool zero_flag, int ndigits, int edigits)
68 char *out;
69 char *digits;
70 int e;
71 char expchar, rchar;
72 format_token ft;
73 int w;
74 int d;
75 /* Number of digits before the decimal point. */
76 int nbefore;
77 /* Number of zeros after the decimal point. */
78 int nzero;
79 /* Number of digits after the decimal point. */
80 int nafter;
81 /* Number of zeros after the decimal point, whatever the precision. */
82 int nzero_real;
83 int leadzero;
84 int nblanks;
85 int i;
86 sign_t sign;
88 ft = f->format;
89 w = f->u.real.w;
90 d = f->u.real.d;
92 rchar = '5';
93 nzero_real = -1;
95 /* We should always know the field width and precision. */
96 if (d < 0)
97 internal_error (&dtp->common, "Unspecified precision");
99 sign = calculate_sign (dtp, sign_bit);
101 /* The following code checks the given string has punctuation in the correct
102 places. Uncomment if needed for debugging.
103 if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
104 || buffer[ndigits + 2] != 'e'))
105 internal_error (&dtp->common, "printf is broken"); */
107 /* Read the exponent back in. */
108 e = atoi (&buffer[ndigits + 3]) + 1;
110 /* Make sure zero comes out as 0.0e0. */
111 if (zero_flag)
113 e = 0;
114 if (compile_options.sign_zero == 1)
115 sign = calculate_sign (dtp, sign_bit);
116 else
117 sign = calculate_sign (dtp, 0);
119 /* Handle special cases. */
120 if (w == 0)
121 w = d + 2;
123 /* For this one we choose to not output a decimal point.
124 F95 10.5.1.2.1 */
125 if (w == 1 && ft == FMT_F)
127 out = write_block (dtp, w);
128 if (out == NULL)
129 return;
130 *out = '0';
131 return;
136 /* Normalize the fractional component. */
137 buffer[2] = buffer[1];
138 digits = &buffer[2];
140 /* Figure out where to place the decimal point. */
141 switch (ft)
143 case FMT_F:
144 if (d == 0 && e <= 0 && dtp->u.p.scale_factor == 0)
146 memmove (digits + 1, digits, ndigits - 1);
147 digits[0] = '0';
148 e++;
151 nbefore = e + dtp->u.p.scale_factor;
152 if (nbefore < 0)
154 nzero = -nbefore;
155 nzero_real = nzero;
156 if (nzero > d)
157 nzero = d;
158 nafter = d - nzero;
159 nbefore = 0;
161 else
163 nzero = 0;
164 nafter = d;
166 expchar = 0;
167 break;
169 case FMT_E:
170 case FMT_D:
171 i = dtp->u.p.scale_factor;
172 if (d <= 0 && i == 0)
174 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
175 "greater than zero in format specifier 'E' or 'D'");
176 return;
178 if (i <= -d || i >= d + 2)
180 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
181 "out of range in format specifier 'E' or 'D'");
182 return;
185 if (!zero_flag)
186 e -= i;
187 if (i < 0)
189 nbefore = 0;
190 nzero = -i;
191 nafter = d + i;
193 else if (i > 0)
195 nbefore = i;
196 nzero = 0;
197 nafter = (d - i) + 1;
199 else /* i == 0 */
201 nbefore = 0;
202 nzero = 0;
203 nafter = d;
206 if (ft == FMT_E)
207 expchar = 'E';
208 else
209 expchar = 'D';
210 break;
212 case FMT_EN:
213 /* The exponent must be a multiple of three, with 1-3 digits before
214 the decimal point. */
215 if (!zero_flag)
216 e--;
217 if (e >= 0)
218 nbefore = e % 3;
219 else
221 nbefore = (-e) % 3;
222 if (nbefore != 0)
223 nbefore = 3 - nbefore;
225 e -= nbefore;
226 nbefore++;
227 nzero = 0;
228 nafter = d;
229 expchar = 'E';
230 break;
232 case FMT_ES:
233 if (!zero_flag)
234 e--;
235 nbefore = 1;
236 nzero = 0;
237 nafter = d;
238 expchar = 'E';
239 break;
241 default:
242 /* Should never happen. */
243 internal_error (&dtp->common, "Unexpected format token");
246 /* Round the value. The value being rounded is an unsigned magnitude.
247 The ROUND_COMPATIBLE is rounding away from zero when there is a tie. */
248 switch (dtp->u.p.current_unit->round_status)
250 case ROUND_ZERO: /* Do nothing and truncation occurs. */
251 goto skip;
252 case ROUND_UP:
253 if (sign_bit)
254 goto skip;
255 rchar = '0';
256 break;
257 case ROUND_DOWN:
258 if (!sign_bit)
259 goto skip;
260 rchar = '0';
261 break;
262 case ROUND_NEAREST:
263 /* Round compatible unless there is a tie. A tie is a 5 with
264 all trailing zero's. */
265 i = nafter + nbefore;
266 if (digits[i] == '5')
268 for(i++ ; i < ndigits; i++)
270 if (digits[i] != '0')
271 goto do_rnd;
273 /* It is a tie so round to even. */
274 switch (digits[nafter + nbefore - 1])
276 case '1':
277 case '3':
278 case '5':
279 case '7':
280 case '9':
281 /* If odd, round away from zero to even. */
282 break;
283 default:
284 /* If even, skip rounding, truncate to even. */
285 goto skip;
288 /* Fall through. */
289 case ROUND_PROCDEFINED:
290 case ROUND_UNSPECIFIED:
291 case ROUND_COMPATIBLE:
292 rchar = '5';
293 /* Just fall through and do the actual rounding. */
296 do_rnd:
298 if (nbefore + nafter == 0)
300 ndigits = 0;
301 if (nzero_real == d && digits[0] >= rchar)
303 /* We rounded to zero but shouldn't have */
304 nzero--;
305 nafter = 1;
306 digits[0] = '1';
307 ndigits = 1;
310 else if (nbefore + nafter < ndigits)
312 ndigits = nbefore + nafter;
313 i = ndigits;
314 if (digits[i] >= rchar)
316 /* Propagate the carry. */
317 for (i--; i >= 0; i--)
319 if (digits[i] != '9')
321 digits[i]++;
322 break;
324 digits[i] = '0';
327 if (i < 0)
329 /* The carry overflowed. Fortunately we have some spare
330 space at the start of the buffer. We may discard some
331 digits, but this is ok because we already know they are
332 zero. */
333 digits--;
334 digits[0] = '1';
335 if (ft == FMT_F)
337 if (nzero > 0)
339 nzero--;
340 nafter++;
342 else
343 nbefore++;
345 else if (ft == FMT_EN)
347 nbefore++;
348 if (nbefore == 4)
350 nbefore = 1;
351 e += 3;
354 else
355 e++;
360 skip:
362 /* Calculate the format of the exponent field. */
363 if (expchar)
365 edigits = 1;
366 for (i = abs (e); i >= 10; i /= 10)
367 edigits++;
369 if (f->u.real.e < 0)
371 /* Width not specified. Must be no more than 3 digits. */
372 if (e > 999 || e < -999)
373 edigits = -1;
374 else
376 edigits = 4;
377 if (e > 99 || e < -99)
378 expchar = ' ';
381 else
383 /* Exponent width specified, check it is wide enough. */
384 if (edigits > f->u.real.e)
385 edigits = -1;
386 else
387 edigits = f->u.real.e + 2;
390 else
391 edigits = 0;
393 /* Zero values always output as positive, even if the value was negative
394 before rounding. */
395 for (i = 0; i < ndigits; i++)
397 if (digits[i] != '0')
398 break;
400 if (i == ndigits)
402 /* The output is zero, so set the sign according to the sign bit unless
403 -fno-sign-zero was specified. */
404 if (compile_options.sign_zero == 1)
405 sign = calculate_sign (dtp, sign_bit);
406 else
407 sign = calculate_sign (dtp, 0);
410 /* Pick a field size if none was specified. */
411 if (w <= 0)
412 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
414 /* Work out how much padding is needed. */
415 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
416 if (sign != S_NONE)
417 nblanks--;
419 if (dtp->u.p.g0_no_blanks)
421 w -= nblanks;
422 nblanks = 0;
425 /* Create the ouput buffer. */
426 out = write_block (dtp, w);
427 if (out == NULL)
428 return;
430 /* Check the value fits in the specified field width. */
431 if (nblanks < 0 || edigits == -1)
433 star_fill (out, w);
434 return;
437 /* See if we have space for a zero before the decimal point. */
438 if (nbefore == 0 && nblanks > 0)
440 leadzero = 1;
441 nblanks--;
443 else
444 leadzero = 0;
446 /* Pad to full field width. */
448 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
450 memset (out, ' ', nblanks);
451 out += nblanks;
454 /* Output the initial sign (if any). */
455 if (sign == S_PLUS)
456 *(out++) = '+';
457 else if (sign == S_MINUS)
458 *(out++) = '-';
460 /* Output an optional leading zero. */
461 if (leadzero)
462 *(out++) = '0';
464 /* Output the part before the decimal point, padding with zeros. */
465 if (nbefore > 0)
467 if (nbefore > ndigits)
469 i = ndigits;
470 memcpy (out, digits, i);
471 ndigits = 0;
472 while (i < nbefore)
473 out[i++] = '0';
475 else
477 i = nbefore;
478 memcpy (out, digits, i);
479 ndigits -= i;
482 digits += i;
483 out += nbefore;
486 /* Output the decimal point. */
487 *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
489 /* Output leading zeros after the decimal point. */
490 if (nzero > 0)
492 for (i = 0; i < nzero; i++)
493 *(out++) = '0';
496 /* Output digits after the decimal point, padding with zeros. */
497 if (nafter > 0)
499 if (nafter > ndigits)
500 i = ndigits;
501 else
502 i = nafter;
504 memcpy (out, digits, i);
505 while (i < nafter)
506 out[i++] = '0';
508 digits += i;
509 ndigits -= i;
510 out += nafter;
513 /* Output the exponent. */
514 if (expchar)
516 if (expchar != ' ')
518 *(out++) = expchar;
519 edigits--;
521 #if HAVE_SNPRINTF
522 snprintf (buffer, size, "%+0*d", edigits, e);
523 #else
524 sprintf (buffer, "%+0*d", edigits, e);
525 #endif
526 memcpy (out, buffer, edigits);
529 if (dtp->u.p.no_leading_blank)
531 out += edigits;
532 memset( out , ' ' , nblanks );
533 dtp->u.p.no_leading_blank = 0;
536 #undef STR
537 #undef STR1
538 #undef MIN_FIELD_WIDTH
542 /* Write "Infinite" or "Nan" as appropriate for the given format. */
544 static void
545 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
547 char * p, fin;
548 int nb = 0;
550 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
552 nb = f->u.real.w;
554 /* If the field width is zero, the processor must select a width
555 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
557 if (nb == 0) nb = 4;
558 p = write_block (dtp, nb);
559 if (p == NULL)
560 return;
561 if (nb < 3)
563 memset (p, '*',nb);
564 return;
567 memset(p, ' ', nb);
568 if (!isnan_flag)
570 if (sign_bit)
573 /* If the sign is negative and the width is 3, there is
574 insufficient room to output '-Inf', so output asterisks */
576 if (nb == 3)
578 memset (p, '*',nb);
579 return;
582 /* The negative sign is mandatory */
584 fin = '-';
586 else
588 /* The positive sign is optional, but we output it for
589 consistency */
590 fin = '+';
592 if (nb > 8)
594 /* We have room, so output 'Infinity' */
595 memcpy(p + nb - 8, "Infinity", 8);
596 else
598 /* For the case of width equals 8, there is not enough room
599 for the sign and 'Infinity' so we go with 'Inf' */
600 memcpy(p + nb - 3, "Inf", 3);
602 if (nb < 9 && nb > 3)
603 p[nb - 4] = fin; /* Put the sign in front of Inf */
604 else if (nb > 8)
605 p[nb - 9] = fin; /* Put the sign in front of Infinity */
607 else
608 memcpy(p + nb - 3, "NaN", 3);
609 return;
614 /* Returns the value of 10**d. */
616 #define CALCULATE_EXP(x) \
617 inline static GFC_REAL_ ## x \
618 calculate_exp_ ## x (int d)\
620 int i;\
621 GFC_REAL_ ## x r = 1.0;\
622 for (i = 0; i< (d >= 0 ? d : -d); i++)\
623 r *= 10;\
624 r = (d >= 0) ? r : 1.0 / r;\
625 return r;\
628 CALCULATE_EXP(4)
630 CALCULATE_EXP(8)
632 #ifdef HAVE_GFC_REAL_10
633 CALCULATE_EXP(10)
634 #endif
636 #ifdef HAVE_GFC_REAL_16
637 CALCULATE_EXP(16)
638 #endif
639 #undef CALCULATE_EXP
641 /* Generate corresponding I/O format for FMT_G and output.
642 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
643 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
645 Data Magnitude Equivalent Conversion
646 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
647 m = 0 F(w-n).(d-1), n' '
648 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
649 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
650 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
651 ................ ..........
652 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
653 m >= 10**d-0.5 Ew.d[Ee]
655 notes: for Gw.d , n' ' means 4 blanks
656 for Gw.dEe, n' ' means e+2 blanks */
658 #define OUTPUT_FLOAT_FMT_G(x) \
659 static void \
660 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
661 GFC_REAL_ ## x m, char *buffer, size_t size, \
662 int sign_bit, bool zero_flag, int ndigits, int edigits) \
664 int e = f->u.real.e;\
665 int d = f->u.real.d;\
666 int w = f->u.real.w;\
667 fnode *newf;\
668 GFC_REAL_ ## x rexp_d;\
669 int low, high, mid;\
670 int ubound, lbound;\
671 char *p;\
672 int save_scale_factor, nb = 0;\
674 save_scale_factor = dtp->u.p.scale_factor;\
675 newf = (fnode *) get_mem (sizeof (fnode));\
677 rexp_d = calculate_exp_ ## x (-d);\
678 if ((m > 0.0 && m < 0.1 - 0.05 * rexp_d) || (rexp_d * (m + 0.5) >= 1.0) ||\
679 ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
681 newf->format = FMT_E;\
682 newf->u.real.w = w;\
683 newf->u.real.d = d;\
684 newf->u.real.e = e;\
685 nb = 0;\
686 goto finish;\
689 mid = 0;\
690 low = 0;\
691 high = d + 1;\
692 lbound = 0;\
693 ubound = d + 1;\
695 while (low <= high)\
697 GFC_REAL_ ## x temp;\
698 mid = (low + high) / 2;\
700 temp = (calculate_exp_ ## x (mid - 1) * (1 - 0.5 * rexp_d));\
702 if (m < temp)\
704 ubound = mid;\
705 if (ubound == lbound + 1)\
706 break;\
707 high = mid - 1;\
709 else if (m > temp)\
711 lbound = mid;\
712 if (ubound == lbound + 1)\
714 mid ++;\
715 break;\
717 low = mid + 1;\
719 else\
721 mid++;\
722 break;\
726 if (e < 0)\
727 nb = 4;\
728 else\
729 nb = e + 2;\
731 newf->format = FMT_F;\
732 newf->u.real.w = f->u.real.w - nb;\
734 if (m == 0.0)\
735 newf->u.real.d = d - 1;\
736 else\
737 newf->u.real.d = - (mid - d - 1);\
739 dtp->u.p.scale_factor = 0;\
741 finish:\
742 output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
743 edigits);\
744 dtp->u.p.scale_factor = save_scale_factor;\
746 free (newf);\
748 if (nb > 0 && !dtp->u.p.g0_no_blanks)\
750 p = write_block (dtp, nb);\
751 if (p == NULL)\
752 return;\
753 memset (p, ' ', nb);\
757 OUTPUT_FLOAT_FMT_G(4)
759 OUTPUT_FLOAT_FMT_G(8)
761 #ifdef HAVE_GFC_REAL_10
762 OUTPUT_FLOAT_FMT_G(10)
763 #endif
765 #ifdef HAVE_GFC_REAL_16
766 OUTPUT_FLOAT_FMT_G(16)
767 #endif
769 #undef OUTPUT_FLOAT_FMT_G
772 /* Define a macro to build code for write_float. */
774 /* Note: Before output_float is called, sprintf is used to print to buffer the
775 number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
776 (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
777 before the decimal point.
779 # The result will always contain a decimal point, even if no
780 digits follow it
782 - The converted value is to be left adjusted on the field boundary
784 + A sign (+ or -) always be placed before a number
786 MIN_FIELD_WIDTH minimum field width
788 * (ndigits-1) is used as the precision
790 e format: [-]d.ddde±dd where there is one digit before the
791 decimal-point character and the number of digits after it is
792 equal to the precision. The exponent always contains at least two
793 digits; if the value is zero, the exponent is 00. */
795 #ifdef HAVE_SNPRINTF
797 #define DTOA \
798 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
799 "e", ndigits - 1, tmp);
801 #define DTOAL \
802 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
803 "Le", ndigits - 1, tmp);
805 #else
807 #define DTOA \
808 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
809 "e", ndigits - 1, tmp);
811 #define DTOAL \
812 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
813 "Le", ndigits - 1, tmp);
815 #endif
817 #define WRITE_FLOAT(x,y)\
819 GFC_REAL_ ## x tmp;\
820 tmp = * (GFC_REAL_ ## x *)source;\
821 sign_bit = __builtin_signbit (tmp);\
822 if (!isfinite (tmp))\
824 write_infnan (dtp, f, isnan (tmp), sign_bit);\
825 return;\
827 tmp = sign_bit ? -tmp : tmp;\
828 zero_flag = (tmp == 0.0);\
830 DTOA ## y\
832 if (f->format != FMT_G)\
833 output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
834 edigits);\
835 else \
836 output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
837 zero_flag, ndigits, edigits);\
840 /* Output a real number according to its format. */
842 static void
843 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
846 #if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
847 # define MIN_FIELD_WIDTH 46
848 #else
849 # define MIN_FIELD_WIDTH 31
850 #endif
851 #define STR(x) STR1(x)
852 #define STR1(x) #x
854 /* This must be large enough to accurately hold any value. */
855 char buffer[MIN_FIELD_WIDTH+1];
856 int sign_bit, ndigits, edigits;
857 bool zero_flag;
858 size_t size;
860 size = MIN_FIELD_WIDTH+1;
862 /* printf pads blanks for us on the exponent so we just need it big enough
863 to handle the largest number of exponent digits expected. */
864 edigits=4;
866 if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
867 || ((f->format == FMT_D || f->format == FMT_E)
868 && dtp->u.p.scale_factor != 0))
870 /* Always convert at full precision to avoid double rounding. */
871 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
873 else
875 /* The number of digits is known, so let printf do the rounding. */
876 if (f->format == FMT_ES)
877 ndigits = f->u.real.d + 1;
878 else
879 ndigits = f->u.real.d;
880 if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
881 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
884 switch (len)
886 case 4:
887 WRITE_FLOAT(4,)
888 break;
890 case 8:
891 WRITE_FLOAT(8,)
892 break;
894 #ifdef HAVE_GFC_REAL_10
895 case 10:
896 WRITE_FLOAT(10,L)
897 break;
898 #endif
899 #ifdef HAVE_GFC_REAL_16
900 case 16:
901 WRITE_FLOAT(16,L)
902 break;
903 #endif
904 default:
905 internal_error (NULL, "bad real kind");