re PR libfortran/36857 (Non-English locale breaks gfortran float formatting ("printf...
[gcc.git] / libgfortran / io / write_float.def
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
5
6 This file is part of the GNU Fortran 95 runtime library (libgfortran).
7
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.
12
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.)
21
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.
26
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. */
31
32 #include "config.h"
33
34 typedef enum
35 { S_NONE, S_MINUS, S_PLUS }
36 sign_t;
37
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. */
40
41 static sign_t
42 calculate_sign (st_parameter_dt *dtp, int negative_flag)
43 {
44 sign_t s = S_NONE;
45
46 if (negative_flag)
47 s = S_MINUS;
48 else
49 switch (dtp->u.p.sign_status)
50 {
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 s = options.optional_plus ? S_PLUS : S_NONE;
59 break;
60 }
61
62 return s;
63 }
64
65
66 /* Output a real number according to its format which is FMT_G free. */
67
68 static void
69 output_float (st_parameter_dt *dtp, const fnode *f, char *buffer, size_t size,
70 int sign_bit, bool zero_flag, int ndigits, int edigits)
71 {
72 char *out;
73 char *digits;
74 int e;
75 char expchar;
76 format_token ft;
77 int w;
78 int d;
79 /* Number of digits before the decimal point. */
80 int nbefore;
81 /* Number of zeros after the decimal point. */
82 int nzero;
83 /* Number of digits after the decimal point. */
84 int nafter;
85 /* Number of zeros after the decimal point, whatever the precision. */
86 int nzero_real;
87 int leadzero;
88 int nblanks;
89 int i;
90 sign_t sign;
91
92 ft = f->format;
93 w = f->u.real.w;
94 d = f->u.real.d;
95
96 nzero_real = -1;
97
98 /* We should always know the field width and precision. */
99 if (d < 0)
100 internal_error (&dtp->common, "Unspecified precision");
101
102 sign = calculate_sign (dtp, sign_bit);
103
104 /* The following code checks the given string has punctuation in the correct
105 places. Uncomment if needed for debugging.
106 if (d != 0 && ((buffer[2] != '.' && buffer[2] != ',')
107 || buffer[ndigits + 2] != 'e'))
108 internal_error (&dtp->common, "printf is broken"); */
109
110 /* Read the exponent back in. */
111 e = atoi (&buffer[ndigits + 3]) + 1;
112
113 /* Make sure zero comes out as 0.0e0. */
114 if (zero_flag)
115 {
116 e = 0;
117 if (compile_options.sign_zero == 1)
118 sign = calculate_sign (dtp, sign_bit);
119 else
120 sign = calculate_sign (dtp, 0);
121 }
122
123 /* Normalize the fractional component. */
124 buffer[2] = buffer[1];
125 digits = &buffer[2];
126
127 /* Figure out where to place the decimal point. */
128 switch (ft)
129 {
130 case FMT_F:
131 nbefore = e + dtp->u.p.scale_factor;
132 if (nbefore < 0)
133 {
134 nzero = -nbefore;
135 nzero_real = nzero;
136 if (nzero > d)
137 nzero = d;
138 nafter = d - nzero;
139 nbefore = 0;
140 }
141 else
142 {
143 nzero = 0;
144 nafter = d;
145 }
146 expchar = 0;
147 break;
148
149 case FMT_E:
150 case FMT_D:
151 i = dtp->u.p.scale_factor;
152 if (d <= 0 && i == 0)
153 {
154 generate_error (&dtp->common, LIBERROR_FORMAT, "Precision not "
155 "greater than zero in format specifier 'E' or 'D'");
156 return;
157 }
158 if (i <= -d || i >= d + 2)
159 {
160 generate_error (&dtp->common, LIBERROR_FORMAT, "Scale factor "
161 "out of range in format specifier 'E' or 'D'");
162 return;
163 }
164
165 if (!zero_flag)
166 e -= i;
167 if (i < 0)
168 {
169 nbefore = 0;
170 nzero = -i;
171 nafter = d + i;
172 }
173 else if (i > 0)
174 {
175 nbefore = i;
176 nzero = 0;
177 nafter = (d - i) + 1;
178 }
179 else /* i == 0 */
180 {
181 nbefore = 0;
182 nzero = 0;
183 nafter = d;
184 }
185
186 if (ft == FMT_E)
187 expchar = 'E';
188 else
189 expchar = 'D';
190 break;
191
192 case FMT_EN:
193 /* The exponent must be a multiple of three, with 1-3 digits before
194 the decimal point. */
195 if (!zero_flag)
196 e--;
197 if (e >= 0)
198 nbefore = e % 3;
199 else
200 {
201 nbefore = (-e) % 3;
202 if (nbefore != 0)
203 nbefore = 3 - nbefore;
204 }
205 e -= nbefore;
206 nbefore++;
207 nzero = 0;
208 nafter = d;
209 expchar = 'E';
210 break;
211
212 case FMT_ES:
213 if (!zero_flag)
214 e--;
215 nbefore = 1;
216 nzero = 0;
217 nafter = d;
218 expchar = 'E';
219 break;
220
221 default:
222 /* Should never happen. */
223 internal_error (&dtp->common, "Unexpected format token");
224 }
225
226 /* Round the value. */
227 if (nbefore + nafter == 0)
228 {
229 ndigits = 0;
230 if (nzero_real == d && digits[0] >= '5')
231 {
232 /* We rounded to zero but shouldn't have */
233 nzero--;
234 nafter = 1;
235 digits[0] = '1';
236 ndigits = 1;
237 }
238 }
239 else if (nbefore + nafter < ndigits)
240 {
241 ndigits = nbefore + nafter;
242 i = ndigits;
243 if (digits[i] >= '5')
244 {
245 /* Propagate the carry. */
246 for (i--; i >= 0; i--)
247 {
248 if (digits[i] != '9')
249 {
250 digits[i]++;
251 break;
252 }
253 digits[i] = '0';
254 }
255
256 if (i < 0)
257 {
258 /* The carry overflowed. Fortunately we have some spare space
259 at the start of the buffer. We may discard some digits, but
260 this is ok because we already know they are zero. */
261 digits--;
262 digits[0] = '1';
263 if (ft == FMT_F)
264 {
265 if (nzero > 0)
266 {
267 nzero--;
268 nafter++;
269 }
270 else
271 nbefore++;
272 }
273 else if (ft == FMT_EN)
274 {
275 nbefore++;
276 if (nbefore == 4)
277 {
278 nbefore = 1;
279 e += 3;
280 }
281 }
282 else
283 e++;
284 }
285 }
286 }
287
288 /* Calculate the format of the exponent field. */
289 if (expchar)
290 {
291 edigits = 1;
292 for (i = abs (e); i >= 10; i /= 10)
293 edigits++;
294
295 if (f->u.real.e < 0)
296 {
297 /* Width not specified. Must be no more than 3 digits. */
298 if (e > 999 || e < -999)
299 edigits = -1;
300 else
301 {
302 edigits = 4;
303 if (e > 99 || e < -99)
304 expchar = ' ';
305 }
306 }
307 else
308 {
309 /* Exponent width specified, check it is wide enough. */
310 if (edigits > f->u.real.e)
311 edigits = -1;
312 else
313 edigits = f->u.real.e + 2;
314 }
315 }
316 else
317 edigits = 0;
318
319 /* Pick a field size if none was specified. */
320 if (w <= 0)
321 w = nbefore + nzero + nafter + (sign != S_NONE ? 2 : 1);
322
323 /* Create the ouput buffer. */
324 out = write_block (dtp, w);
325 if (out == NULL)
326 return;
327
328 /* Zero values always output as positive, even if the value was negative
329 before rounding. */
330 for (i = 0; i < ndigits; i++)
331 {
332 if (digits[i] != '0')
333 break;
334 }
335 if (i == ndigits)
336 {
337 /* The output is zero, so set the sign according to the sign bit unless
338 -fno-sign-zero was specified. */
339 if (compile_options.sign_zero == 1)
340 sign = calculate_sign (dtp, sign_bit);
341 else
342 sign = calculate_sign (dtp, 0);
343 }
344
345 /* Work out how much padding is needed. */
346 nblanks = w - (nbefore + nzero + nafter + edigits + 1);
347 if (sign != S_NONE)
348 nblanks--;
349
350 /* Check the value fits in the specified field width. */
351 if (nblanks < 0 || edigits == -1)
352 {
353 star_fill (out, w);
354 return;
355 }
356
357 /* See if we have space for a zero before the decimal point. */
358 if (nbefore == 0 && nblanks > 0)
359 {
360 leadzero = 1;
361 nblanks--;
362 }
363 else
364 leadzero = 0;
365
366 /* Pad to full field width. */
367
368 if ( ( nblanks > 0 ) && !dtp->u.p.no_leading_blank)
369 {
370 memset (out, ' ', nblanks);
371 out += nblanks;
372 }
373
374 /* Output the initial sign (if any). */
375 if (sign == S_PLUS)
376 *(out++) = '+';
377 else if (sign == S_MINUS)
378 *(out++) = '-';
379
380 /* Output an optional leading zero. */
381 if (leadzero)
382 *(out++) = '0';
383
384 /* Output the part before the decimal point, padding with zeros. */
385 if (nbefore > 0)
386 {
387 if (nbefore > ndigits)
388 {
389 i = ndigits;
390 memcpy (out, digits, i);
391 ndigits = 0;
392 while (i < nbefore)
393 out[i++] = '0';
394 }
395 else
396 {
397 i = nbefore;
398 memcpy (out, digits, i);
399 ndigits -= i;
400 }
401
402 digits += i;
403 out += nbefore;
404 }
405 /* Output the decimal point. */
406 *(out++) = dtp->u.p.decimal_status == DECIMAL_POINT ? '.' : ',';
407
408 /* Output leading zeros after the decimal point. */
409 if (nzero > 0)
410 {
411 for (i = 0; i < nzero; i++)
412 *(out++) = '0';
413 }
414
415 /* Output digits after the decimal point, padding with zeros. */
416 if (nafter > 0)
417 {
418 if (nafter > ndigits)
419 i = ndigits;
420 else
421 i = nafter;
422
423 memcpy (out, digits, i);
424 while (i < nafter)
425 out[i++] = '0';
426
427 digits += i;
428 ndigits -= i;
429 out += nafter;
430 }
431
432 /* Output the exponent. */
433 if (expchar)
434 {
435 if (expchar != ' ')
436 {
437 *(out++) = expchar;
438 edigits--;
439 }
440 #if HAVE_SNPRINTF
441 snprintf (buffer, size, "%+0*d", edigits, e);
442 #else
443 sprintf (buffer, "%+0*d", edigits, e);
444 #endif
445 memcpy (out, buffer, edigits);
446 }
447 if (dtp->u.p.no_leading_blank)
448 {
449 out += edigits;
450 memset( out , ' ' , nblanks );
451 dtp->u.p.no_leading_blank = 0;
452 }
453 #undef STR
454 #undef STR1
455 #undef MIN_FIELD_WIDTH
456 }
457
458
459 /* Write "Infinite" or "Nan" as appropriate for the given format. */
460
461 static void
462 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
463 {
464 char * p, fin;
465 int nb = 0;
466
467 if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
468 {
469 nb = f->u.real.w;
470
471 /* If the field width is zero, the processor must select a width
472 not zero. 4 is chosen to allow output of '-Inf' or '+Inf' */
473
474 if (nb == 0) nb = 4;
475 p = write_block (dtp, nb);
476 if (p == NULL)
477 return;
478 if (nb < 3)
479 {
480 memset (p, '*',nb);
481 return;
482 }
483
484 memset(p, ' ', nb);
485 if (!isnan_flag)
486 {
487 if (sign_bit)
488 {
489
490 /* If the sign is negative and the width is 3, there is
491 insufficient room to output '-Inf', so output asterisks */
492
493 if (nb == 3)
494 {
495 memset (p, '*',nb);
496 return;
497 }
498
499 /* The negative sign is mandatory */
500
501 fin = '-';
502 }
503 else
504
505 /* The positive sign is optional, but we output it for
506 consistency */
507 fin = '+';
508
509 if (nb > 8)
510
511 /* We have room, so output 'Infinity' */
512 memcpy(p + nb - 8, "Infinity", 8);
513 else
514
515 /* For the case of width equals 8, there is not enough room
516 for the sign and 'Infinity' so we go with 'Inf' */
517 memcpy(p + nb - 3, "Inf", 3);
518
519 if (nb < 9 && nb > 3)
520 p[nb - 4] = fin; /* Put the sign in front of Inf */
521 else if (nb > 8)
522 p[nb - 9] = fin; /* Put the sign in front of Infinity */
523 }
524 else
525 memcpy(p + nb - 3, "NaN", 3);
526 return;
527 }
528 }
529
530
531 /* Returns the value of 10**d. */
532
533 #define CALCULATE_EXP(x) \
534 inline static GFC_REAL_ ## x \
535 calculate_exp_ ## x (int d)\
536 {\
537 int i;\
538 GFC_REAL_ ## x r = 1.0;\
539 for (i = 0; i< (d >= 0 ? d : -d); i++)\
540 r *= 10;\
541 r = (d >= 0) ? r : 1.0 / r;\
542 return r;\
543 }
544
545 CALCULATE_EXP(4)
546
547 CALCULATE_EXP(8)
548
549 #ifdef HAVE_GFC_REAL_10
550 CALCULATE_EXP(10)
551 #endif
552
553 #ifdef HAVE_GFC_REAL_16
554 CALCULATE_EXP(16)
555 #endif
556 #undef CALCULATE_EXP
557
558 /* Generate corresponding I/O format for FMT_G and output.
559 The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
560 LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
561
562 Data Magnitude Equivalent Conversion
563 0< m < 0.1-0.5*10**(-d-1) Ew.d[Ee]
564 m = 0 F(w-n).(d-1), n' '
565 0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d) F(w-n).d, n' '
566 1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1) F(w-n).(d-1), n' '
567 10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2) F(w-n).(d-2), n' '
568 ................ ..........
569 10**(d-1)-0.5*10**(-1)<= m <10**d-0.5 F(w-n).0,n(' ')
570 m >= 10**d-0.5 Ew.d[Ee]
571
572 notes: for Gw.d , n' ' means 4 blanks
573 for Gw.dEe, n' ' means e+2 blanks */
574
575 #define OUTPUT_FLOAT_FMT_G(x) \
576 static void \
577 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
578 GFC_REAL_ ## x m, char *buffer, size_t size, \
579 int sign_bit, bool zero_flag, int ndigits, int edigits) \
580 { \
581 int e = f->u.real.e;\
582 int d = f->u.real.d;\
583 int w = f->u.real.w;\
584 fnode *newf;\
585 GFC_REAL_ ## x exp_d;\
586 int low, high, mid;\
587 int ubound, lbound;\
588 char *p;\
589 int save_scale_factor, nb = 0;\
590 \
591 save_scale_factor = dtp->u.p.scale_factor;\
592 newf = get_mem (sizeof (fnode));\
593 \
594 exp_d = calculate_exp_ ## x (d);\
595 if ((m > 0.0 && m < 0.1 - 0.05 / exp_d) || (m >= exp_d - 0.5 ) ||\
596 ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
597 { \
598 newf->format = FMT_E;\
599 newf->u.real.w = w;\
600 newf->u.real.d = d;\
601 newf->u.real.e = e;\
602 nb = 0;\
603 goto finish;\
604 }\
605 \
606 mid = 0;\
607 low = 0;\
608 high = d + 1;\
609 lbound = 0;\
610 ubound = d + 1;\
611 \
612 while (low <= high)\
613 { \
614 GFC_REAL_ ## x temp;\
615 mid = (low + high) / 2;\
616 \
617 temp = 0.1 * calculate_exp_ ## x (mid) - 0.5\
618 * calculate_exp_ ## x (mid - d - 1);\
619 \
620 if (m < temp)\
621 { \
622 ubound = mid;\
623 if (ubound == lbound + 1)\
624 break;\
625 high = mid - 1;\
626 }\
627 else if (m > temp)\
628 { \
629 lbound = mid;\
630 if (ubound == lbound + 1)\
631 { \
632 mid ++;\
633 break;\
634 }\
635 low = mid + 1;\
636 }\
637 else\
638 break;\
639 }\
640 \
641 if (e < 0)\
642 nb = 4;\
643 else\
644 nb = e + 2;\
645 \
646 newf->format = FMT_F;\
647 newf->u.real.w = f->u.real.w - nb;\
648 \
649 if (m == 0.0)\
650 newf->u.real.d = d - 1;\
651 else\
652 newf->u.real.d = - (mid - d - 1);\
653 \
654 dtp->u.p.scale_factor = 0;\
655 \
656 finish:\
657 output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
658 edigits);\
659 dtp->u.p.scale_factor = save_scale_factor;\
660 \
661 free_mem(newf);\
662 \
663 if (nb > 0)\
664 { \
665 p = write_block (dtp, nb);\
666 if (p == NULL)\
667 return;\
668 memset (p, ' ', nb);\
669 }\
670 }\
671
672 OUTPUT_FLOAT_FMT_G(4)
673
674 OUTPUT_FLOAT_FMT_G(8)
675
676 #ifdef HAVE_GFC_REAL_10
677 OUTPUT_FLOAT_FMT_G(10)
678 #endif
679
680 #ifdef HAVE_GFC_REAL_16
681 OUTPUT_FLOAT_FMT_G(16)
682 #endif
683
684 #undef OUTPUT_FLOAT_FMT_G
685
686
687 /* Define a macro to build code for write_float. */
688
689 /* Note: Before output_float is called, sprintf is used to print to buffer the
690 number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
691 (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
692 before the decimal point.
693
694 # The result will always contain a decimal point, even if no
695 digits follow it
696
697 - The converted value is to be left adjusted on the field boundary
698
699 + A sign (+ or -) always be placed before a number
700
701 MIN_FIELD_WIDTH minimum field width
702
703 * (ndigits-1) is used as the precision
704
705 e format: [-]d.ddde±dd where there is one digit before the
706 decimal-point character and the number of digits after it is
707 equal to the precision. The exponent always contains at least two
708 digits; if the value is zero, the exponent is 00. */
709
710 #ifdef HAVE_SNPRINTF
711
712 #define DTOA \
713 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
714 "e", ndigits - 1, tmp);
715
716 #define DTOAL \
717 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
718 "Le", ndigits - 1, tmp);
719
720 #else
721
722 #define DTOA \
723 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
724 "e", ndigits - 1, tmp);
725
726 #define DTOAL \
727 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
728 "Le", ndigits - 1, tmp);
729
730 #endif
731
732 #define WRITE_FLOAT(x,y)\
733 {\
734 GFC_REAL_ ## x tmp;\
735 tmp = * (GFC_REAL_ ## x *)source;\
736 sign_bit = signbit (tmp);\
737 if (!isfinite (tmp))\
738 { \
739 write_infnan (dtp, f, isnan (tmp), sign_bit);\
740 return;\
741 }\
742 tmp = sign_bit ? -tmp : tmp;\
743 if (f->u.real.d == 0 && f->format == FMT_F)\
744 {\
745 if (tmp < 0.5)\
746 tmp = 0.0;\
747 else if (tmp < 1.0)\
748 tmp = tmp + 0.5;\
749 }\
750 zero_flag = (tmp == 0.0);\
751 \
752 DTOA ## y\
753 \
754 if (f->format != FMT_G)\
755 output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
756 edigits);\
757 else \
758 output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
759 zero_flag, ndigits, edigits);\
760 }\
761
762 /* Output a real number according to its format. */
763
764 static void
765 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
766 {
767
768 #if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
769 # define MIN_FIELD_WIDTH 46
770 #else
771 # define MIN_FIELD_WIDTH 31
772 #endif
773 #define STR(x) STR1(x)
774 #define STR1(x) #x
775
776 /* This must be large enough to accurately hold any value. */
777 char buffer[MIN_FIELD_WIDTH+1];
778 int sign_bit, ndigits, edigits;
779 bool zero_flag;
780 size_t size;
781
782 size = MIN_FIELD_WIDTH+1;
783
784 /* printf pads blanks for us on the exponent so we just need it big enough
785 to handle the largest number of exponent digits expected. */
786 edigits=4;
787
788 if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G
789 || ((f->format == FMT_D || f->format == FMT_E)
790 && dtp->u.p.scale_factor != 0))
791 {
792 /* Always convert at full precision to avoid double rounding. */
793 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
794 }
795 else
796 {
797 /* The number of digits is known, so let printf do the rounding. */
798 if (f->format == FMT_ES)
799 ndigits = f->u.real.d + 1;
800 else
801 ndigits = f->u.real.d;
802 if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
803 ndigits = MIN_FIELD_WIDTH - 4 - edigits;
804 }
805
806 switch (len)
807 {
808 case 4:
809 WRITE_FLOAT(4,)
810 break;
811
812 case 8:
813 WRITE_FLOAT(8,)
814 break;
815
816 #ifdef HAVE_GFC_REAL_10
817 case 10:
818 WRITE_FLOAT(10,L)
819 break;
820 #endif
821 #ifdef HAVE_GFC_REAL_16
822 case 16:
823 WRITE_FLOAT(16,L)
824 break;
825 #endif
826 default:
827 internal_error (NULL, "bad real kind");
828 }
829 }