OSDN Git Service

2008-09-22 Jerry DeLisle <jvdelisle@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       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 }