OSDN Git Service

73a6ed14a1b290b74043198294a5dec2ff5deb25
[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   *(out++) = dtp->u.p.current_unit->decimal_status == DECIMAL_POINT ? '.' : ',';
408
409   /* Output leading zeros after the decimal point.  */
410   if (nzero > 0)
411     {
412       for (i = 0; i < nzero; i++)
413         *(out++) = '0';
414     }
415
416   /* Output digits after the decimal point, padding with zeros.  */
417   if (nafter > 0)
418     {
419       if (nafter > ndigits)
420         i = ndigits;
421       else
422         i = nafter;
423
424       memcpy (out, digits, i);
425       while (i < nafter)
426         out[i++] = '0';
427
428       digits += i;
429       ndigits -= i;
430       out += nafter;
431     }
432
433   /* Output the exponent.  */
434   if (expchar)
435     {
436       if (expchar != ' ')
437         {
438           *(out++) = expchar;
439           edigits--;
440         }
441 #if HAVE_SNPRINTF
442       snprintf (buffer, size, "%+0*d", edigits, e);
443 #else
444       sprintf (buffer, "%+0*d", edigits, e);
445 #endif
446       memcpy (out, buffer, edigits);
447     }
448   if (dtp->u.p.no_leading_blank)
449     {
450       out += edigits;
451       memset( out , ' ' , nblanks );
452       dtp->u.p.no_leading_blank = 0;
453     }
454 #undef STR
455 #undef STR1
456 #undef MIN_FIELD_WIDTH
457 }
458
459
460 /* Write "Infinite" or "Nan" as appropriate for the given format.  */
461
462 static void
463 write_infnan (st_parameter_dt *dtp, const fnode *f, int isnan_flag, int sign_bit)
464 {
465   char * p, fin;
466   int nb = 0;
467
468   if (f->format != FMT_B && f->format != FMT_O && f->format != FMT_Z)
469     {
470           nb =  f->u.real.w;
471           
472           /* If the field width is zero, the processor must select a width 
473              not zero.  4 is chosen to allow output of '-Inf' or '+Inf' */
474              
475           if (nb == 0) nb = 4;
476           p = write_block (dtp, nb);
477           if (p == NULL)
478             return;
479           if (nb < 3)
480             {
481               memset (p, '*',nb);
482               return;
483             }
484
485           memset(p, ' ', nb);
486           if (!isnan_flag)
487             {
488               if (sign_bit)
489                 {
490                 
491                   /* If the sign is negative and the width is 3, there is
492                      insufficient room to output '-Inf', so output asterisks */
493                      
494                   if (nb == 3)
495                     {
496                       memset (p, '*',nb);
497                       return;
498                     }
499                     
500                   /* The negative sign is mandatory */
501                     
502                   fin = '-';
503                 }    
504               else
505               
506                   /* The positive sign is optional, but we output it for
507                      consistency */
508                   fin = '+';
509
510               if (nb > 8)
511               
512                 /* We have room, so output 'Infinity' */
513                 memcpy(p + nb - 8, "Infinity", 8);
514               else
515               
516                 /* For the case of width equals 8, there is not enough room
517                    for the sign and 'Infinity' so we go with 'Inf' */
518                 memcpy(p + nb - 3, "Inf", 3);
519
520               if (nb < 9 && nb > 3)
521                 p[nb - 4] = fin;  /* Put the sign in front of Inf */
522               else if (nb > 8)
523                 p[nb - 9] = fin;  /* Put the sign in front of Infinity */
524             }
525           else
526             memcpy(p + nb - 3, "NaN", 3);
527           return;
528         }
529     }
530
531
532 /* Returns the value of 10**d.  */
533
534 #define CALCULATE_EXP(x) \
535 inline static GFC_REAL_ ## x \
536 calculate_exp_ ## x  (int d)\
537 {\
538   int i;\
539   GFC_REAL_ ## x r = 1.0;\
540   for (i = 0; i< (d >= 0 ? d : -d); i++)\
541     r *= 10;\
542   r = (d >= 0) ? r : 1.0 / r;\
543   return r;\
544 }
545
546 CALCULATE_EXP(4)
547
548 CALCULATE_EXP(8)
549
550 #ifdef HAVE_GFC_REAL_10
551 CALCULATE_EXP(10)
552 #endif
553
554 #ifdef HAVE_GFC_REAL_16
555 CALCULATE_EXP(16)
556 #endif
557 #undef CALCULATE_EXP
558
559 /* Generate corresponding I/O format for FMT_G and output.
560    The rules to translate FMT_G to FMT_E or FMT_F from DEC fortran
561    LRM (table 11-2, Chapter 11, "I/O Formatting", P11-25) is:
562
563    Data Magnitude                              Equivalent Conversion
564    0< m < 0.1-0.5*10**(-d-1)                   Ew.d[Ee]
565    m = 0                                       F(w-n).(d-1), n' '
566    0.1-0.5*10**(-d-1)<= m < 1-0.5*10**(-d)     F(w-n).d, n' '
567    1-0.5*10**(-d)<= m < 10-0.5*10**(-d+1)      F(w-n).(d-1), n' '
568    10-0.5*10**(-d+1)<= m < 100-0.5*10**(-d+2)  F(w-n).(d-2), n' '
569    ................                           ..........
570    10**(d-1)-0.5*10**(-1)<= m <10**d-0.5       F(w-n).0,n(' ')
571    m >= 10**d-0.5                              Ew.d[Ee]
572
573    notes: for Gw.d ,  n' ' means 4 blanks
574           for Gw.dEe, n' ' means e+2 blanks  */
575
576 #define OUTPUT_FLOAT_FMT_G(x) \
577 static void \
578 output_float_FMT_G_ ## x (st_parameter_dt *dtp, const fnode *f, \
579                       GFC_REAL_ ## x m, char *buffer, size_t size, \
580                       int sign_bit, bool zero_flag, int ndigits, int edigits) \
581 { \
582   int e = f->u.real.e;\
583   int d = f->u.real.d;\
584   int w = f->u.real.w;\
585   fnode *newf;\
586   GFC_REAL_ ## x exp_d;\
587   int low, high, mid;\
588   int ubound, lbound;\
589   char *p;\
590   int save_scale_factor, nb = 0;\
591 \
592   save_scale_factor = dtp->u.p.scale_factor;\
593   newf = get_mem (sizeof (fnode));\
594 \
595   exp_d = calculate_exp_ ## x (d);\
596   if ((m > 0.0 && m < 0.1 - 0.05 / exp_d) || (m >= exp_d - 0.5 ) ||\
597       ((m == 0.0) && !(compile_options.allow_std & GFC_STD_F2003)))\
598     { \
599       newf->format = FMT_E;\
600       newf->u.real.w = w;\
601       newf->u.real.d = d;\
602       newf->u.real.e = e;\
603       nb = 0;\
604       goto finish;\
605     }\
606 \
607   mid = 0;\
608   low = 0;\
609   high = d + 1;\
610   lbound = 0;\
611   ubound = d + 1;\
612 \
613   while (low <= high)\
614     { \
615       GFC_REAL_ ## x temp;\
616       mid = (low + high) / 2;\
617 \
618       temp = 0.1 * calculate_exp_ ## x (mid) - 0.5\
619              * calculate_exp_ ## x (mid - d - 1);\
620 \
621       if (m < temp)\
622         { \
623           ubound = mid;\
624           if (ubound == lbound + 1)\
625             break;\
626           high = mid - 1;\
627         }\
628       else if (m > temp)\
629         { \
630           lbound = mid;\
631           if (ubound == lbound + 1)\
632             { \
633               mid ++;\
634               break;\
635             }\
636           low = mid + 1;\
637         }\
638       else\
639         break;\
640     }\
641 \
642   if (e < 0)\
643     nb = 4;\
644   else\
645     nb = e + 2;\
646 \
647   newf->format = FMT_F;\
648   newf->u.real.w = f->u.real.w - nb;\
649 \
650   if (m == 0.0)\
651     newf->u.real.d = d - 1;\
652   else\
653     newf->u.real.d = - (mid - d - 1);\
654 \
655   dtp->u.p.scale_factor = 0;\
656 \
657  finish:\
658   output_float (dtp, newf, buffer, size, sign_bit, zero_flag, ndigits, \
659                 edigits);\
660   dtp->u.p.scale_factor = save_scale_factor;\
661 \
662   free_mem(newf);\
663 \
664   if (nb > 0)\
665     { \
666       p = write_block (dtp, nb);\
667       if (p == NULL)\
668         return;\
669       memset (p, ' ', nb);\
670     }\
671 }\
672
673 OUTPUT_FLOAT_FMT_G(4)
674
675 OUTPUT_FLOAT_FMT_G(8)
676
677 #ifdef HAVE_GFC_REAL_10
678 OUTPUT_FLOAT_FMT_G(10)
679 #endif
680
681 #ifdef HAVE_GFC_REAL_16
682 OUTPUT_FLOAT_FMT_G(16)
683 #endif
684
685 #undef OUTPUT_FLOAT_FMT_G
686
687
688 /* Define a macro to build code for write_float.  */
689
690   /* Note: Before output_float is called, sprintf is used to print to buffer the
691      number in the format +D.DDDDe+ddd. For an N digit exponent, this gives us
692      (MIN_FIELD_WIDTH-5)-N digits after the decimal point, plus another one
693      before the decimal point.
694
695      #   The result will always contain a decimal point, even if no
696          digits follow it
697
698      -   The converted value is to be left adjusted on the field boundary
699
700      +   A sign (+ or -) always be placed before a number
701
702      MIN_FIELD_WIDTH  minimum field width
703
704      *   (ndigits-1) is used as the precision
705
706      e format: [-]d.ddde±dd where there is one digit before the
707        decimal-point character and the number of digits after it is
708        equal to the precision. The exponent always contains at least two
709        digits; if the value is zero, the exponent is 00.  */
710
711 #ifdef HAVE_SNPRINTF
712
713 #define DTOA \
714 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
715           "e", ndigits - 1, tmp);
716
717 #define DTOAL \
718 snprintf (buffer, size, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
719           "Le", ndigits - 1, tmp);
720
721 #else
722
723 #define DTOA \
724 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
725          "e", ndigits - 1, tmp);
726
727 #define DTOAL \
728 sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
729          "Le", ndigits - 1, tmp);
730
731 #endif
732
733 #define WRITE_FLOAT(x,y)\
734 {\
735         GFC_REAL_ ## x tmp;\
736         tmp = * (GFC_REAL_ ## x *)source;\
737         sign_bit = signbit (tmp);\
738         if (!isfinite (tmp))\
739           { \
740             write_infnan (dtp, f, isnan (tmp), sign_bit);\
741             return;\
742           }\
743         tmp = sign_bit ? -tmp : tmp;\
744         if (f->u.real.d == 0 && f->format == FMT_F)\
745           {\
746             if (tmp < 0.5)\
747               tmp = 0.0;\
748             else if (tmp < 1.0)\
749               tmp = 1.0;\
750           }\
751         zero_flag = (tmp == 0.0);\
752 \
753         DTOA ## y\
754 \
755         if (f->format != FMT_G)\
756           output_float (dtp, f, buffer, size, sign_bit, zero_flag, ndigits, \
757                         edigits);\
758         else \
759           output_float_FMT_G_ ## x (dtp, f, tmp, buffer, size, sign_bit, \
760                                     zero_flag, ndigits, edigits);\
761 }\
762
763 /* Output a real number according to its format.  */
764
765 static void
766 write_float (st_parameter_dt *dtp, const fnode *f, const char *source, int len)
767 {
768
769 #if defined(HAVE_GFC_REAL_16) && __LDBL_DIG__ > 18
770 # define MIN_FIELD_WIDTH 46
771 #else
772 # define MIN_FIELD_WIDTH 31
773 #endif
774 #define STR(x) STR1(x)
775 #define STR1(x) #x
776
777   /* This must be large enough to accurately hold any value.  */
778   char buffer[MIN_FIELD_WIDTH+1];
779   int sign_bit, ndigits, edigits;
780   bool zero_flag;
781   size_t size;
782
783   size = MIN_FIELD_WIDTH+1;
784
785   /* printf pads blanks for us on the exponent so we just need it big enough
786      to handle the largest number of exponent digits expected.  */
787   edigits=4;
788
789   if (f->format == FMT_F || f->format == FMT_EN || f->format == FMT_G 
790       || ((f->format == FMT_D || f->format == FMT_E)
791       && dtp->u.p.scale_factor != 0))
792     {
793       /* Always convert at full precision to avoid double rounding.  */
794       ndigits = MIN_FIELD_WIDTH - 4 - edigits;
795     }
796   else
797     {
798       /* The number of digits is known, so let printf do the rounding.  */
799       if (f->format == FMT_ES)
800         ndigits = f->u.real.d + 1;
801       else
802         ndigits = f->u.real.d;
803       if (ndigits > MIN_FIELD_WIDTH - 4 - edigits)
804         ndigits = MIN_FIELD_WIDTH - 4 - edigits;
805     }
806
807   switch (len)
808     {
809     case 4:
810       WRITE_FLOAT(4,)
811       break;
812
813     case 8:
814       WRITE_FLOAT(8,)
815       break;
816
817 #ifdef HAVE_GFC_REAL_10
818     case 10:
819       WRITE_FLOAT(10,L)
820       break;
821 #endif
822 #ifdef HAVE_GFC_REAL_16
823     case 16:
824       WRITE_FLOAT(16,L)
825       break;
826 #endif
827     default:
828       internal_error (NULL, "bad real kind");
829     }
830 }