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