home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / glibc-1.06 / sysdeps / generic / printf_fp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-04-22  |  10.3 KB  |  492 lines
  1. /* Floating-point printing for `printf'.
  2.    This is an implementation of a restricted form of the `Dragon4'
  3.    algorithm described in "How to Print Floating-Point Numbers Accurately",
  4.    by Guy L. Steele, Jr. and Jon L. White, presented at the ACM SIGPLAN '90
  5.    Conference on Programming Language Design and Implementation.
  6.  
  7. Copyright (C) 1992, 1993 Free Software Foundation, Inc.
  8. This file is part of the GNU C Library.
  9.  
  10. The GNU C Library is free software; you can redistribute it and/or
  11. modify it under the terms of the GNU Library General Public License as
  12. published by the Free Software Foundation; either version 2 of the
  13. License, or (at your option) any later version.
  14.  
  15. The GNU C Library is distributed in the hope that it will be useful,
  16. but WITHOUT ANY WARRANTY; without even the implied warranty of
  17. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18. Library General Public License for more details.
  19.  
  20. You should have received a copy of the GNU Library General Public
  21. License along with the GNU C Library; see the file COPYING.LIB.  If
  22. not, write to the Free Software Foundation, Inc., 675 Mass Ave,
  23. Cambridge, MA 02139, USA.  */
  24.  
  25. #include <ansidecl.h>
  26. #include <ctype.h>
  27. #include <stdio.h>
  28. #include <float.h>
  29. #include <math.h>
  30. #include <stdarg.h>
  31. #include <stdlib.h>
  32. #include <localeinfo.h>
  33.  
  34. #include <printf.h>
  35.  
  36. #define    outchar(x)                                  \
  37.   do                                          \
  38.     {                                          \
  39.       register CONST int outc = (x);                          \
  40.       if (putc (outc, s) == EOF)                          \
  41.     return -1;                                  \
  42.       else                                      \
  43.     ++done;                                      \
  44.     } while (0)
  45.  
  46. #if FLT_RADIX != 2
  47.  
  48. double
  49. frexp (double f, int *e)
  50. {
  51.   #error "Don't know how to extract fraction and exponent from `double'."
  52. }
  53.  
  54. #undef    ldexp
  55. #ifdef    __GNUC__
  56. inline
  57. #endif
  58. static double
  59. ldexp (double f, int e)
  60. {
  61.   while (e > 0)
  62.     {
  63.       f *= FLT_RADIX;
  64.       --e;
  65.     }
  66.   while (e < 0)
  67.     {
  68.       f /= FLT_RADIX;
  69.       ++e;
  70.     }
  71. }
  72.  
  73. #endif
  74.  
  75.  
  76. int
  77. DEFUN(__printf_fp, (s, info, args),
  78.       FILE *s AND CONST struct printf_info *info AND va_list *args)
  79. {
  80.   int done = 0;
  81.  
  82.   /* Decimal point character.  */
  83.   CONST char *CONST decimal = _numeric_info->decimal_point;
  84.  
  85.   LONG_DOUBLE fpnum;        /* Input.  */
  86.   int is_neg;
  87.  
  88.   LONG_DOUBLE f;        /* Fraction.  */
  89.  
  90.   int e;            /* Base-2 exponent of the input.  */
  91.   CONST int p = DBL_MANT_DIG;    /* Internal precision.  */
  92.   LONG_DOUBLE scale, scale10;    /* Scale factor.  */
  93.   LONG_DOUBLE loerr, hierr;    /* Potential error in the fraction.  */
  94.   int k;            /* Digits to the left of the decimal point.  */
  95.   int cutoff;            /* Where to stop generating digits.  */
  96.   LONG_DOUBLE r, r2, r10;    /* Remainder.  */
  97.   int roundup;
  98.   int low, high;
  99.   char digit;
  100.  
  101.   int j;
  102.  
  103.   char type = tolower (info->spec);
  104.   int prec = info->prec;
  105.   int width = info->width;
  106.  
  107.   /* This algorithm has the nice property of not needing a buffer.
  108.      However, to get the padding right for %g format, we need to know
  109.      the length of the number before printing it.  */
  110.  
  111. #ifndef    LDBL_DIG
  112. #define    LDBL_DIG    DBL_DIG
  113. #endif
  114. #ifndef    LDBL_MAX_10_EXP
  115. #define    LDBL_MAX_10_EXP    DBL_MAX_10_EXP
  116. #endif
  117.  
  118.   char *buf = __alloca ((prec > LDBL_DIG ? prec : LDBL_DIG) +
  119.             LDBL_MAX_10_EXP + 3); /* Dot, e, exp. sign.  */
  120.   register char *bp = buf;
  121. #define    put(c)    *bp++ = (c)
  122.  
  123.   /* Fetch the argument value.  */
  124.   if (info->is_long_double)
  125.     fpnum = va_arg (*args, LONG_DOUBLE);
  126.   else
  127.     fpnum = (LONG_DOUBLE) va_arg (*args, double);
  128.  
  129. #ifdef    HANDLE_SPECIAL
  130.   /* Allow for machine-dependent (or floating point format-dependent) code.  */
  131.   HANDLE_SPECIAL (done, s, info, fpnum);
  132. #endif
  133.  
  134. #ifndef    IS_NEGATIVE
  135. #define    IS_NEGATIVE(num)    ((num) < 0)
  136. #endif
  137.  
  138.   is_neg = IS_NEGATIVE (fpnum);
  139.   if (is_neg)
  140.     fpnum = - fpnum;
  141.  
  142.   if (prec == -1)
  143.     prec = 6;
  144.   
  145.   if (type == 'g')
  146.     {
  147.       if (prec == 0)
  148.     prec = 1;
  149.  
  150.       if (fpnum != 0)
  151.     {
  152.       if (fpnum < 1e-4)
  153.         type = 'e';
  154.       else
  155.         {
  156.           /* Is (int) floor (log10 (FPNUM)) >= PREC?  */
  157.           j = prec;
  158.           if (j > p)
  159.         j = p;
  160.           f = 10;
  161.           while (--j > 0)
  162.         {
  163.           f *= 10;
  164.           if (fpnum < f)
  165.             /* log10 (F) == floor (log10 (FPNUM)) + 1
  166.                log10 (FPNUM) == Number of iterations minus one.  */
  167.             break;
  168.         }
  169.           if (j <= 0)
  170.         /* We got all the way through the loop and F (i.e., 10**J)
  171.            never reached FPNUM, so we want to use %e format.  */
  172.         type = 'e';
  173.         }
  174.     }
  175.  
  176.       /* For 'g'/'G' format, the precision specifies "significant digits",
  177.      not digits to come after the decimal point.  */
  178.       --prec;
  179.     }
  180.  
  181.   if (fpnum == 0)
  182.     /* Special case for zero.
  183.        The general algorithm does not work for zero.  */
  184.     {
  185.       put ('0');
  186.       if (tolower (info->spec) != 'g' || info->alt)
  187.     {
  188.       if (prec > 0 || info->alt)
  189.         put (*decimal);
  190.       while (prec-- > 0)
  191.         put ('0');
  192.     }
  193.       if (type == 'e')
  194.     {
  195.       put (info->spec);
  196.       put ('+');
  197.       put ('0');
  198.       put ('0');
  199.     }
  200.     }
  201.   else
  202.     {
  203.       /* Split the number into a fraction and base-2 exponent.  */
  204.       f = frexp (fpnum, &e);
  205.  
  206.       /* Scale the fractional part by the highest possible number of
  207.      significant bits of fraction.  We want to represent the
  208.      fractional part as a (very) large integer.  */
  209.       f = ldexp (f, p);
  210.  
  211.       cutoff = -prec;
  212.  
  213.       roundup = 0;
  214.  
  215.       if (e > p)
  216.     {
  217.       /* The exponent is bigger than the number of fractional digits.  */
  218.       r = ldexp (f, e - p);
  219.       scale = 1;
  220.       /* The number is (E - P) factors of two larger than
  221.          the fraction can represent; this is the potential error.  */
  222.       loerr = ldexp (1.0, e - p);
  223.     }
  224.       else
  225.     {
  226.       /* The number of fractional digits is greater than the exponent.
  227.          Scale by the difference factors of two.  */
  228.       r = f;
  229.       scale = ldexp (1.0, p - e);
  230.       loerr = 1.0;
  231.     }
  232.       hierr = loerr;
  233.  
  234.       /* Fixup.  */
  235.  
  236.       if (f == ldexp (1.0, p - 1))
  237.     {
  238.       /* Account for unequal gaps.  */
  239.       hierr = ldexp (hierr, 1);
  240.       r = ldexp (r, 1);
  241.       scale = ldexp (scale, 1);
  242.     }
  243.  
  244.       scale10 = ceil (scale / 10.0);
  245.       k = 0;
  246.       while (r < scale10)
  247.     {
  248.       --k;
  249.       r *= 10;
  250.       loerr *= 10;
  251.       hierr *= 10;
  252.     }
  253.       do
  254.     {
  255.       r2 = 2 * r;
  256.  
  257.       while (r2 + hierr >= 2 * scale)
  258.         {
  259.           scale *= 10;
  260.           ++k;
  261.         }
  262.  
  263.       /* Perform any necessary adjustment of loerr and hierr to
  264.          take into account the formatting requirements.  */
  265.  
  266.       if (type == 'e')
  267.         cutoff += k - 1;    /* CutOffMode == "relative".  */
  268.       /* Otherwise CutOffMode == "absolute".  */
  269.  
  270.       {            /* CutOffAdjust.  */
  271.         int a = cutoff - k;
  272.         double y = scale;
  273.         while (a > 0)
  274.           {
  275.         y *= 10;
  276.         --a;
  277.           }
  278.         while (a < 0)
  279.           {
  280.         y = ceil (y / 10);
  281.         ++a;
  282.           }
  283.         /* y == ceil (scale * pow (10.0, (double) (cutoff - k))) */
  284.         if (y > loerr)
  285.           loerr = y;
  286.         if (y > hierr)
  287.           {
  288.         hierr = y;
  289.         roundup = 1;
  290.           }
  291.       }            /* End CutOffAdjust.  */
  292.     } while (r2 + hierr >= 2 * scale);
  293.  
  294.       /* End Fixup.  */
  295.  
  296.       /* First digit.  */
  297.       --k;
  298.       r10 = r * 10;
  299.       digit = '0' + (unsigned int) floor (r10 / scale);
  300.       r = fmod (r10, scale);
  301.       loerr *= 10;
  302.       hierr *= 10;
  303.       
  304.       low = 2 * r < loerr;
  305.       if (roundup)
  306.     high = 2 * r >= (2 * scale) - hierr;
  307.       else
  308.     high = 2 * r > (2 * scale) - hierr;
  309.       
  310.       if (low || high || k == cutoff)
  311.     {
  312.       if ((high && !low) || (2 * r > scale))
  313.         ++digit;
  314.     }
  315.  
  316.       if (type == 'e')
  317.     {
  318.       /* Exponential notation.  */
  319.  
  320.       int expt = k;        /* Base-10 exponent.  */
  321.       int expt_neg;
  322.  
  323.       expt_neg = k < 0;
  324.       if (expt_neg)
  325.         expt = - expt;
  326.  
  327.       /* Find the magnitude of the exponent.  */
  328.       j = 10;
  329.       while (j <= expt)
  330.         j *= 10;
  331.  
  332.       /* Write the first digit.  */
  333.       put (digit);
  334.  
  335.       if (low || high || k == cutoff)
  336.         {
  337.           if ((tolower (info->spec) != 'g' && prec > 0) || info->alt)
  338.         put (*decimal);
  339.         }
  340.       else
  341.         {
  342.           put (*decimal);
  343.  
  344.           /* Remaining digits.  */
  345.           while (1)
  346.         {
  347.           --k;
  348.           r10 = r * 10;
  349.           digit = '0' + (unsigned int) floor (r10 / scale);
  350.           r = fmod (r10, scale);
  351.           loerr *= 10;
  352.           hierr *= 10;
  353.               
  354.           low = 2 * r < loerr;
  355.           if (roundup)
  356.             high = 2 * r >= (2 * scale) - hierr;
  357.           else
  358.             high = 2 * r > (2 * scale) - hierr;
  359.               
  360.           if (low || high || k == cutoff)
  361.             {
  362.               if ((high && !low) || (2 * r > scale))
  363.             ++digit;
  364.               put (digit);
  365.               break;
  366.             }
  367.               
  368.           put (digit);
  369.         }
  370.         }
  371.  
  372.       if (tolower (info->spec) != 'g' || info->alt)
  373.         /* Pad with zeros.  */
  374.         while (--k >= cutoff)
  375.           put ('0');
  376.  
  377.       /* Write the exponent.  */
  378.       put (isupper (info->spec) ? 'E' : 'e');
  379.       put (expt_neg ? '-' : '+');
  380.       if (expt < 10)
  381.         /* Exponent always has at least two digits.  */
  382.         put ('0');
  383.       do
  384.         {
  385.           j /= 10;
  386.           put ('0' + (expt / j));
  387.           expt %= j;
  388.         }
  389.       while (j > 1);
  390.     }
  391.       else
  392.     {
  393.       /* Decimal fraction notation.  */
  394.  
  395.       if (k < 0)
  396.         {
  397.           put ('0');
  398.           if (prec > 0 || info->alt)
  399.         put (*decimal);
  400.  
  401.           /* Write leading fractional zeros.  */
  402.           j = 0;
  403.           while (--j > k)
  404.         put ('0');
  405.         }
  406.  
  407.       if (low || high || k == cutoff)
  408.         put (digit);
  409.       else
  410.         while (1)
  411.           {
  412.         put (digit);
  413.         
  414.         --k;
  415.         r10 = r * 10;
  416.         digit = '0' + (unsigned int) floor (r10 / scale);
  417.         r = fmod (r10, scale);
  418.         loerr *= 10;
  419.         hierr *= 10;
  420.         
  421.         low = 2 * r < loerr;
  422.         if (roundup)
  423.           high = 2 * r >= (2 * scale) - hierr;
  424.         else
  425.           high = 2 * r > (2 * scale) - hierr;
  426.         
  427.         if (k == -1)
  428.           put (*decimal);
  429.  
  430.         if (low || high || k == cutoff)
  431.           {
  432.             if ((high && !low) || (2 * r > scale))
  433.               ++digit;
  434.             put (digit);
  435.             break;
  436.           }
  437.           }
  438.       while (k > 0)
  439.         {
  440.           put ('0');
  441.           --k;
  442.         }
  443.       if ((type != 'g' && prec > 0) || info->alt)
  444.         {
  445.           if (k == 0)
  446.         put (*decimal);
  447.           while (k-- > -prec)
  448.         put ('0');
  449.         }
  450.     }
  451.     }
  452.  
  453. #undef    put
  454.  
  455.   /* The number is all converted in BUF.
  456.      Now write it with sign and appropriate padding.  */
  457.  
  458.   if (is_neg || info->showsign || info->space)
  459.     --width;
  460.  
  461.   width -= bp - buf;
  462.  
  463.   if (!info->left && info->pad == ' ')
  464.     /* Pad with spaces on the left.  */
  465.     while (width-- > 0)
  466.       outchar (' ');
  467.  
  468.   /* Write the sign.  */
  469.   if (is_neg)
  470.     outchar ('-');
  471.   else if (info->showsign)
  472.     outchar ('+');
  473.   else if (info->space)
  474.     outchar (' ');
  475.  
  476.   if (!info->left && info->pad == '0')
  477.     /* Pad with zeros on the left.  */
  478.     while (width-- > 0)
  479.       outchar ('0');
  480.  
  481.   if (fwrite (buf, bp - buf, 1, s) != 1)
  482.     return -1;
  483.   done += bp - buf;
  484.  
  485.   if (info->left)
  486.     /* Pad with spaces on the right.  */
  487.     while (width-- > 0)
  488.       outchar (' ');
  489.  
  490.   return done;
  491. }
  492.