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