Merged revisions 143552,143554,143557,143560,143562,143564-143567,143570-143573,14357...
[official-gcc.git] / libgfortran / io / write_float.def
blob17bd7844560125513b5472d6729132d2955a6a06
1 /* Copyright (C) 2007, 2008 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 95 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 2, or (at your option)
11 any later version.
13 In addition to the permissions in the GNU General Public License, the
14 Free Software Foundation gives you unlimited permission to link the
15 compiled version of this file into combinations with other programs,
16 and to distribute those combinations without any restriction coming
17 from the use of this file. (The General Public License restrictions
18 do apply in other respects; for example, they cover modification of
19 the file, and distribution when not linked into a combine
20 executable.)
22 Libgfortran is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 GNU General Public License for more details.
27 You should have received a copy of the GNU General Public License
28 along with Libgfortran; see the file COPYING. If not, write to
29 the Free Software Foundation, 51 Franklin Street, Fifth Floor,
30 Boston, MA 02110-1301, USA. */
32 #include "config.h"
34 typedef enum
35 { S_NONE, S_MINUS, S_PLUS }
36 sign_t;
38 /* Given a flag that indicates if a value is negative or not, return a
39 sign_t that gives the sign that we need to produce. */
41 static sign_t
42 calculate_sign (st_parameter_dt *dtp, int negative_flag)
44 sign_t s = S_NONE;
46 if (negative_flag)
47 s = S_MINUS;
48 else
49 switch (dtp->u.p.sign_status)
51 case SIGN_SP: /* Show sign. */
52 s = S_PLUS;
53 break;
54 case SIGN_SS: /* Suppress sign. */
55 s = S_NONE;
56 break;
57 case SIGN_S: /* Processor defined. */
58 case SIGN_UNSPECIFIED:
59 s = options.optional_plus ? S_PLUS : S_NONE;
60 break;
63 return s;
67 /* Output a real number according to its format which is FMT_G free. */
69 static void
70 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
71 int sign_bit, bool zero_flag, int ndigits, int edigits)
73 char *out;
74 char *digits;
75 int e;
76 char expchar;
77 format_token ft;
78 int w;
79 int d;
80 /* Number of digits before the decimal point. */
81 int nbefore;
82 /* Number of zeros after the decimal point. */
83 int nzero;
84 /* Number of digits after the decimal point. */
85 int nafter;
86 /* Number of zeros after the decimal point, whatever the precision. */
87 int nzero_real;
88 int leadzero;
89 int nblanks;
90 int i;
91 sign_t sign;
93 ft = f->format;
94 w = f->u.real.w;
95 d = f->u.real.d;
97 nzero_real = -1;
99 /* We should always know the field width and precision. */
100 if (d < 0)
101 internal_error (&dtp->common, "Unspecified precision");
103 sign = calculate_sign (dtp, sign_bit);
105 /* The following code checks the given string has punctuation in the correct
106 places. Uncomment if needed for debugging.
107 if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
108 || buffer[ndigits + 2] != 'e'))
109 internal_error (&dtp->common, "printf is broken"); */
111 /* Read the exponent back in. */
112 e = atoi (&buffer[ndigits + 3]) + 1;
114 /* Make sure zero comes out as 0.0e0. */
115 if (zero_flag)
117 e = 0;
118 if (compile_options.sign_zero == 1)
119 sign = calculate_sign (dtp, sign_bit);
120 else
121 sign = calculate_sign (dtp, 0);
123 /* Handle special cases. */
124 if (w == 0)
125 w = d + 2;
127 /* For this one we choose to not output a decimal point.
128 F95 10.5.1.2.1 */
129 if (w == 1 && ft == FMT_F)
131 out = write_block (dtp, w);
132 if (out == NULL)
133 return;
134 *out = '0';
135 return;
140 /* Normalize the fractional component. */
141 buffer[2] = buffer[1];
142 digits = &buffer[2];
144 /* Figure out where to place the decimal point. */
145 switch (ft)
147 case FMT_F:
148 nbefore = e + dtp->u.p.scale_factor;
149 if (nbefore < 0)
151 nzero = -nbefore;
152 nzero_real = nzero;
153 if (nzero > d)
154 nzero = d;
155 nafter = d - nzero;
156 nbefore = 0;
158 else
160 nzero = 0;
161 nafter = d;
163 expchar = 0;
164 break;
166 case FMT_E:
167 case FMT_D:
168 i = dtp->u.p.scale_factor;
169 if (d <= 0 && i == 0)
171 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
172 "greater than zero in format specifier 'E' or 'D'");
173 return;
175 if (i <= -d || i >= d + 2)
177 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
178 "out of range in format specifier 'E' or 'D'");
179 return;
182 if (!zero_flag)
183 e -= i;
184 if (i < 0)
186 nbefore = 0;
187 nzero = -i;
188 nafter = d + i;
190 else if (i > 0)
192 nbefore = i;
193 nzero = 0;
194 nafter = (d - i) + 1;
196 else /* i == 0 */
198 nbefore = 0;
199 nzero = 0;
200 nafter = d;
203 if (ft == FMT_E)
204 expchar = 'E';
205 else
206 expchar = 'D';
207 break;
209 case FMT_EN:
210 /* The exponent must be a multiple of three, with 1-3 digits before
211 the decimal point. */
212 if (!zero_flag)
213 e--;
214 if (e >= 0)
215 nbefore = e % 3;
216 else
218 nbefore = (-e) % 3;
219 if (nbefore != 0)
220 nbefore = 3 - nbefore;
222 e -= nbefore;
223 nbefore++;
224 nzero = 0;
225 nafter = d;
226 expchar = 'E';
227 break;
229 case FMT_ES:
230 if (!zero_flag)
231 e--;
232 nbefore = 1;
233 nzero = 0;
234 nafter = d;
235 expchar = 'E';
236 break;
238 default:
239 /* Should never happen. */
240 internal_error (&dtp->common, "Unexpected format token");
243 /* Round the value. */
244 if (nbefore + nafter == 0)
246 ndigits = 0;
247 if (nzero_real == d && digits[0] >= '5')
249 /* We rounded to zero but shouldn't have */
250 nzero--;
251 nafter = 1;
252 digits[0] = '1';
253 ndigits = 1;
256 else if (nbefore + nafter < ndigits)
258 ndigits = nbefore + nafter;
259 i = ndigits;
260 if (digits[i] >= '5')
262 /* Propagate the carry. */
263 for (i--; i >= 0; i--)
265 if (digits[i] != '9')
267 digits[i]++;
268 break;
270 digits[i] = '0';
273 if (i < 0)
275 /* The carry overflowed. Fortunately we have some spare space
276 at the start of the buffer. We may discard some digits, but
277 this is ok because we already know they are zero. */
278 digits--;
279 digits[0] = '1';
280 if (ft == FMT_F)
282 if (nzero > 0)
284 nzero--;
285 nafter++;
287 else
288 nbefore++;
290 else if (ft == FMT_EN)
292 nbefore++;
293 if (nbefore == 4)
295 nbefore = 1;
296 e += 3;
299 else
300 e++;
305 /* Calculate the format of the exponent field. */
306 if (expchar)
308 edigits = 1;
309 for (i = abs (e); i >= 10; i /= 10)
310 edigits++;
312 if (f->u.real.e < 0)
314 /* Width not specified. Must be no more than 3 digits. */
315 if (e > 999 || e < -999)
316 edigits = -1;
317 else
319 edigits = 4;
320 if (e > 99 || e < -99)
321 expchar = ' ';
324 else
326 /* Exponent width specified, check it is wide enough. */
327 if (edigits > f->u.real.e)
328 edigits = -1;
329 else
330 edigits = f->u.real.e + 2;
333 else
334 edigits = 0;
336 /* Zero values always output as positive, even if the value was negative
337 before rounding. */
338 for (i = 0; i < ndigits; i++)
340 if (digits[i] != '0')
341 break;
343 if (i == ndigits)
345 /* The output is zero, so set the sign according to the sign bit unless
346 -fno-sign-zero was specified. */
347 if (compile_options.sign_zero == 1)
348 sign = calculate_sign (dtp, sign_bit);
349 else
350 sign = calculate_sign (dtp, 0);
353 /* Pick a field size if none was specified. */
354 if (w <= 0)
355 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
357 /* Work out how much padding is needed. */
358 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
359 if (sign != S_NONE)
360 nblanks--;
362 if (dtp->u.p.g0_no_blanks)
364 w -= nblanks;
365 nblanks = 0;
368 /* Create the ouput buffer. */
369 out = write_block (dtp, w);
370 if (out == NULL)
371 return;
373 /* Check the value fits in the specified field width. */
374 if (nblanks < 0 || edigits == -1)
376 star_fill (out, w);
377 return;
380 /* See if we have space for a zero before the decimal point. */
381 if (nbefore == 0 && nblanks > 0)
383 leadzero = 1;
384 nblanks--;
386 else
387 leadzero = 0;
389 /* Pad to full field width. */
391 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
393 memset (out, ' ', nblanks);
394 out += nblanks;
397 /* Output the initial sign (if any). */
398 if (sign == S_PLUS)
399 *(out++) = '+';
400 else if (sign == S_MINUS)
401 *(out++) = '-';
403 /* Output an optional leading zero. */
404 if (leadzero)
405 *(out++) = '0';
407 /* Output the part before the decimal point, padding with zeros. */
408 if (nbefore > 0)
410 if (nbefore > ndigits)
412 i = ndigits;
413 memcpy (out, digits, i);
414 ndigits = 0;
415 while (i < nbefore)
416 out[i++] = '0';
418 else
420 i = nbefore;
421 memcpy (out, digits, i);
422 ndigits -= i;
425 digits += i;
426 out += nbefore;
429 /* Output the decimal point. */
430 *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
432 /* Output leading zeros after the decimal point. */
433 if (nzero > 0)
435 for (i = 0; i < nzero; i++)
436 *(out++) = '0';
439 /* Output digits after the decimal point, padding with zeros. */
440 if (nafter > 0)
442 if (nafter > ndigits)
443 i = ndigits;
444 else
445 i = nafter;
447 memcpy (out, digits, i);
448 while (i < nafter)
449 out[i++] = '0';
451 digits += i;
452 ndigits -= i;
453 out += nafter;
456 /* Output the exponent. */
457 if (expchar)
459 if (expchar != ' ')
461 *(out++) = expchar;
462 edigits--;
464 #if HAVE_SNPRINTF
465 snprintf (buffer, size, "%+0*d", edigits, e);
466 #else
467 sprintf (buffer, "%+0*d", edigits, e);
468 #endif
469 memcpy (out, buffer, edigits);
472 if (dtp->u.p.no_leading_blank)
474 out += edigits;
475 memset( out , ' ' , nblanks );
476 dtp->u.p.no_leading_blank = 0;
479 #undef STR
480 #undef STR1
481 #undef MIN_FIELD_WIDTH
485 /* Write "Infinite" or "Nan" as appropriate for the given format. */
487 static void
488 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
490 char * p, fin;
491 int nb = 0;
493 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
495 nb = f->u.real.w;
497 /* If the field width is zero, the processor must select a width
498 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
500 if (nb == 0) nb = 4;
501 p = write_block (dtp, nb);
502 if (p == NULL)
503 return;
504 if (nb < 3)
506 memset (p, '*',nb);
507 return;
510 memset(p, ' ', nb);
511 if (!isnan_flag)
513 if (sign_bit)
516 /* If the sign is negative and the width is 3, there is
517 insufficient room to output '-Inf', so output asterisks */
519 if (nb == 3)
521 memset (p, '*',nb);
522 return;
525 /* The negative sign is mandatory */
527 fin = '-';
529 else
531 /* The positive sign is optional, but we output it for
532 consistency */
533 fin = '+';
535 if (nb > 8)
537 /* We have room, so output 'Infinity' */
538 memcpy(p + nb - 8, "Infinity", 8);
539 else
541 /* For the case of width equals 8, there is not enough room
542 for the sign and 'Infinity' so we go with 'Inf' */
543 memcpy(p + nb - 3, "Inf", 3);
545 if (nb < 9 && nb > 3)
546 p[nb - 4] = fin; /* Put the sign in front of Inf */
547 else if (nb > 8)
548 p[nb - 9] = fin; /* Put the sign in front of Infinity */
550 else
551 memcpy(p + nb - 3, "NaN", 3);
552 return;
557 /* Returns the value of 10**d. */
559 #define CALCULATE_EXP(x) \
560 inline static GFC_REAL_ ## x \
561 calculate_exp_ ## x (int d)\
563 int i;\
564 GFC_REAL_ ## x r = 1.0;\
565 for (i = 0; i< (d >= 0 ? d : -d); i++)\
566 r *= 10;\
567 r = (d >= 0) ? r : 1.0 / r;\
568 return r;\
571 CALCULATE_EXP(4)
573 CALCULATE_EXP(8)
575 #ifdef HAVE_GFC_REAL_10
576 CALCULATE_EXP(10)
577 #endif
579 #ifdef HAVE_GFC_REAL_16
580 CALCULATE_EXP(16)
581 #endif
582 #undef CALCULATE_EXP
584 /* Generate corresponding I/O format for FMT_G and output.
585 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
586 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
588 Data Magnitude Equivalent Conversion
589 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
590 m = 0 F(w-n).(d-1), n' '
591 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
592 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
593 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
594 ................ ..........
595 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
596 m >= 10**d-0.5 Ew.d[Ee]
598 notes: for Gw.d , n' ' means 4 blanks
599 for Gw.dEe, n' ' means e+2 blanks */
601 #define OUTPUT_FLOAT_FMT_G(x) \
602 static void \
603 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
604 GFC_REAL_ ## x m, char *buffer, size_t size, \
605 int sign_bit, bool zero_flag, int ndigits, int edigits) \
607 int e = f->u.real.e;\
608 int d = f->u.real.d;\
609 int w = f->u.real.w;\
610 fnode *newf;\
611 GFC_REAL_ ## x exp_d;\
612 int low, high, mid;\
613 int ubound, lbound;\
614 char *p;\
615 int save_scale_factor, nb = 0;\
617 save_scale_factor = dtp->u.p.scale_factor;\
618 newf = (fnode *) get_mem (sizeof (fnode));\
620 exp_d = calculate_exp_ ## x (d);\
621 if ((m > 0.0 && m < 0.1 - 0.05 / exp_d) || (m >= exp_d - 0.5 ) ||\
622 ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
624 newf->format = FMT_E;\
625 newf->u.real.w = w;\
626 newf->u.real.d = d;\
627 newf->u.real.e = e;\
628 nb = 0;\
629 goto finish;\
632 mid = 0;\
633 low = 0;\
634 high = d + 1;\
635 lbound = 0;\
636 ubound = d + 1;\
638 while (low <= high)\
640 GFC_REAL_ ## x temp;\
641 mid = (low + high) / 2;\
643 temp = (calculate_exp_ ## x (mid) - \
644 5 * calculate_exp_ ## x (mid - d - 1)) / 10;\
646 if (m < temp)\
648 ubound = mid;\
649 if (ubound == lbound + 1)\
650 break;\
651 high = mid - 1;\
653 else if (m > temp)\
655 lbound = mid;\
656 if (ubound == lbound + 1)\
658 mid ++;\
659 break;\
661 low = mid + 1;\
663 else\
665 mid++;\
666 break;\
670 if (e < 0)\
671 nb = 4;\
672 else\
673 nb = e + 2;\
675 newf->format = FMT_F;\
676 newf->u.real.w = f->u.real.w - nb;\
678 if (m == 0.0)\
679 newf->u.real.d = d - 1;\
680 else\
681 newf->u.real.d = - (mid - d - 1);\
683 dtp->u.p.scale_factor = 0;\
685 finish:\
686 output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
687 edigits);\
688 dtp->u.p.scale_factor = save_scale_factor;\
690 free_mem(newf);\
692 if (nb > 0 && !dtp->u.p.g0_no_blanks)\
694 p = write_block (dtp, nb);\
695 if (p == NULL)\
696 return;\
697 memset (p, ' ', nb);\
701 OUTPUT_FLOAT_FMT_G(4)
703 OUTPUT_FLOAT_FMT_G(8)
705 #ifdef HAVE_GFC_REAL_10
706 OUTPUT_FLOAT_FMT_G(10)
707 #endif
709 #ifdef HAVE_GFC_REAL_16
710 OUTPUT_FLOAT_FMT_G(16)
711 #endif
713 #undef OUTPUT_FLOAT_FMT_G
716 /* Define a macro to build code for write_float. */
718 /* Note: Before output_float is called, sprintf is used to print to buffer the
719 number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
720 (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
721 before the decimal point.
723 # The result will always contain a decimal point, even if no
724 digits follow it
726 - The converted value is to be left adjusted on the field boundary
728 + A sign (+ or -) always be placed before a number
730 MIN_FIELD_WIDTH minimum field width
732 * (ndigits-1) is used as the precision
734 e format: [-]d.ddde±dd where there is one digit before the
735 decimal-point character and the number of digits after it is
736 equal to the precision. The exponent always contains at least two
737 digits; if the value is zero, the exponent is 00. */
739 #ifdef HAVE_SNPRINTF
741 #define DTOA \
742 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
743 "e", ndigits - 1, tmp);
745 #define DTOAL \
746 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
747 "Le", ndigits - 1, tmp);
749 #else
751 #define DTOA \
752 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
753 "e", ndigits - 1, tmp);
755 #define DTOAL \
756 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
757 "Le", ndigits - 1, tmp);
759 #endif
761 #define WRITE_FLOAT(x,y)\
763 GFC_REAL_ ## x tmp;\
764 tmp = * (GFC_REAL_ ## x *)source;\
765 sign_bit = signbit (tmp);\
766 if (!isfinite (tmp))\
768 write_infnan (dtp, f, isnan (tmp), sign_bit);\
769 return;\
771 tmp = sign_bit ? -tmp : tmp;\
772 if (f->u.real.d == 0 && f->format == FMT_F\
773 && dtp->u.p.scale_factor == 0)\
775 if (tmp < 0.5)\
776 tmp = 0.0;\
777 else if (tmp < 1.0)\
778 tmp = 1.0;\
780 zero_flag = (tmp == 0.0);\
782 DTOA ## y\
784 if (f->format != FMT_G)\
785 output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
786 edigits);\
787 else \
788 output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
789 zero_flag, ndigits, edigits);\
792 /* Output a real number according to its format. */
794 static void
795 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
798 #if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
799 # define MIN_FIELD_WIDTH 46
800 #else
801 # define MIN_FIELD_WIDTH 31
802 #endif
803 #define STR(x) STR1(x)
804 #define STR1(x) #x
806 /* This must be large enough to accurately hold any value. */
807 char buffer[MIN_FIELD_WIDTH+1];
808 int sign_bit, ndigits, edigits;
809 bool zero_flag;
810 size_t size;
812 size = MIN_FIELD_WIDTH+1;
814 /* printf pads blanks for us on the exponent so we just need it big enough
815 to handle the largest number of exponent digits expected. */
816 edigits=4;
818 if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
819 || ((f->format == FMT_D || f->format == FMT_E)
820 && dtp->u.p.scale_factor != 0))
822 /* Always convert at full precision to avoid double rounding. */
823 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
825 else
827 /* The number of digits is known, so let printf do the rounding. */
828 if (f->format == FMT_ES)
829 ndigits = f->u.real.d + 1;
830 else
831 ndigits = f->u.real.d;
832 if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
833 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
836 switch (len)
838 case 4:
839 WRITE_FLOAT(4,)
840 break;
842 case 8:
843 WRITE_FLOAT(8,)
844 break;
846 #ifdef HAVE_GFC_REAL_10
847 case 10:
848 WRITE_FLOAT(10,L)
849 break;
850 #endif
851 #ifdef HAVE_GFC_REAL_16
852 case 16:
853 WRITE_FLOAT(16,L)
854 break;
855 #endif
856 default:
857 internal_error (NULL, "bad real kind");