OSDN Git Service

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