[v3] FP printing, take 2

Gawain Bolton gp.bolton@computer.org
Tue Aug 31 20:10:00 GMT 2004


Hi Jerry,

Yes I have one comment...   Maybe question actually because I must 
admit, I don't know the first thing about the locale part of 
libstdc++!   The following lines strike me as odd:

*__sbuf++ = __digits[*__buf++ - '0'];

    *__tb-- = '0' + (__decp % 10);

Here the implicit assumption seems to be that the character set used 
orders the digits in successive increasing order from '0' to '9', which 
is true of ASCII.

Is this a valid assumption for the locale library?


Gawain

Jerry Quinn wrote:

>Here's a revision.  I fixed the long lines and replaced log10 in
>locale_facets.tcc with __builtin_log10 to avoid including cmath there.
>
>I also uglified everything that shows up as a symbol in floatconv.o.
>
>At this point I see two questions remaining.  What should the linker map
>section actually be called, and does floatconv.cc need to be uglified further.
>
>I think further discussion of GMP can wait until this is in.
>
>Final comments?  Approvals :-)?
>
>
>
>2004-08-29  Jerry Quinn  jlquinn@optonline.net
>
>	* config/linker-map.gnu (GLIBCXX_3.4.2): Add dec_exponent,
>        float_to_char_sprintf, dtoa.
>	* include/bits/locale_facets.tcc (__float_to_char_fixed,
>	__float_to_char_scientific, __float_to_char): New.
>	(__dec_exponent, __dtoa, __float_to_char_sprintf): Declare.
>	(num_put::_M_insert_float): Use __float_to_char instead of
>        __convert_from_v.
>	* Makefile.am (sources):  Add floatconv.cc.
>	* Makefile.in: Regenerate.
>	* src/floatconv.cc.
>
>
>Index: config/linker-map.gnu
>===================================================================
>RCS file: /cvs/gcc/gcc/libstdc++-v3/config/linker-map.gnu,v
>retrieving revision 1.68
>diff -u -r1.68 linker-map.gnu
>--- config/linker-map.gnu	30 Aug 2004 18:22:00 -0000	1.68
>+++ config/linker-map.gnu	31 Aug 2004 00:43:48 -0000
>@@ -259,6 +259,10 @@
>     _ZN9__gnu_cxx17__pool_alloc_base16_M_get_free_listE[jm];
>     _ZN9__gnu_cxx17__pool_alloc_base12_M_get_mutexEv;
> 
>+    _ZSt14__dec_exponentI[de]EiT_RiRb;
>+    _ZSt23__float_to_char_sprintfI[de]EiPcT_iiSt13_Ios_Fmtflags;
>+    _ZSt6__dtoaI[de]EPcS0_T_iiibPi;
>+
> } GLIBCXX_3.4.1;
> 
> # Symbols in the support library (libsupc++) have their own tag.
>Index: include/bits/locale_facets.tcc
>===================================================================
>RCS file: /cvs/gcc/gcc/libstdc++-v3/include/bits/locale_facets.tcc,v
>retrieving revision 1.203
>diff -u -r1.203 locale_facets.tcc
>--- include/bits/locale_facets.tcc	29 Aug 2004 09:04:33 -0000	1.203
>+++ include/bits/locale_facets.tcc	31 Aug 2004 00:43:49 -0000
>@@ -1037,16 +1037,254 @@
>       __len = __newlen;
>     }
> 
>-  // The following code uses snprintf (or sprintf(), when
>-  // _GLIBCXX_USE_C99 is not defined) to convert floating point values
>-  // for insertion into a stream.  An optimization would be to replace
>-  // them with code that works directly on a wide buffer and then use
>-  // __pad to do the padding.  It would be good to replace them anyway
>-  // to gain back the efficiency that C++ provides by knowing up front
>-  // the type of the values to insert.  Also, sprintf is dangerous
>-  // since may lead to accidental buffer overruns.  This
>-  // implementation follows the C++ standard fairly directly as
>-  // outlined in 22.2.2.2 [lib.locale.num.put]
>+  // Prints fixed point number to sbuf
>+  //
>+  // atoms contains numeric and exponent wide chars to use for the result.
>+  // dpchar is the decimal point char to use.
>+  //
>+  // The number consists of buf and decp.  buf contains the sigificant digits
>+  // of the number with no trailing 0's.  decp is the position of the decimal
>+  // point.  It is the number of digits in buf to the left of the decimal
>+  // point.
>+  //
>+  // prec specifies the maximum number of digits to print after the decimal
>+  // point.  flags is provided for showpoint.  showpoint prints a point even
>+  // if there are no more digits to print after the point.
>+  //
>+  // pad is true if exactly prec digits are to be printed after the point.  If
>+  // there are not enough, 0's are added to the end.
>+  template<typename _CharT>
>+    int
>+    __float_to_char_fixed(_CharT* __sbuf, const _CharT* __atoms,
>+			  _CharT __dpchar, const char* __buf, int __decp, 
>+			  int __prec, ios_base::fmtflags __flags, bool __pad)
>+    {
>+      _CharT* __sb = __sbuf;
>+      const _CharT* __digits = __atoms + __num_base::_S_odigits;
>+      bool __showpoint = __flags & ios_base::showpoint;
>+
>+      // Print integer portion
>+      if (__decp > 0)
>+	{
>+	  while (*__buf && __decp) {
>+	    *__sbuf++ = __digits[*__buf++ - '0'];
>+	    __decp--;
>+	  }
>+	  // Fill with 0's if we ran out of digits.  If there were enough,
>+	  // decp will be 0.
>+	  for (; __decp; __decp--)
>+	    *__sbuf++ = __digits[0];
>+	}
>+      else
>+	// Leading 0 for fractional only numbers and zero.
>+	*__sbuf++ = __digits[0];
>+
>+      // Print decimal point
>+      // Showpoint always prints point.  Otherwise when prec > 0, there must
>+      // be digits to print or forced padding.
>+      if ((__prec && (__pad || *__buf))	|| __showpoint)
>+	*__sbuf++ = __dpchar;
>+
>+
>+      // Print 0's between decimal point and first significant digit
>+      if (__decp < 0) {
>+	while (__decp < 0 && __prec > 0) {
>+	  *__sbuf++ = __digits[0];
>+	  __decp++;
>+	  __prec--;
>+	}
>+      }
>+
>+      // Print significant fractional digits
>+      // At this point decp <= 0
>+      for (; __prec > 0 && *__buf; __prec--)
>+	*__sbuf++ = __digits[*__buf++ - '0'];
>+      if (__pad) {
>+	while (__prec-- > 0)
>+	  *__sbuf++ = __digits[0];
>+      }
>+
>+      return __sbuf - __sb;
>+    }
>+
>+  // Prints scientific format number to sbuf
>+  //
>+  // Prints a number in scientific format i.e. +-d.ddddE+ee.
>+  //
>+  // atoms contains numeric and exponent wide chars to use for the result.
>+  // dpchar is the decimal point char to use.
>+  //
>+  // The number consists of buf and decp.
>+  // buf contains the sigificant digits of the number with no trailing 0's.
>+  // decp is the position of the decimal point.  It is the number of digits
>+  // to the left of the decimal point.
>+  //
>+  // prec specifies the maximum number of digits to print after the decimal
>+  // point in scientific form.  showpoint prints a point even if there are no
>+  // more digits to print after the point.
>+  //
>+  // pad is true if exactly prec digits are to be printed after the point.  If
>+  // there are not enough, 0's are added to the end.
>+  template<typename _CharT>
>+    int
>+    __float_to_char_scientific(_CharT* __sbuf, const _CharT* __atoms,
>+			       _CharT __dpchar, const char* __buf, int __decp,
>+			       int __prec, ios_base::fmtflags __flags, 
>+			       bool __pad)
>+    {
>+      const _CharT* __digits = __atoms + __num_base::_S_odigits;
>+      _CharT* __sb = __sbuf;
>+      const char* __bp = __buf;
>+      bool __showpoint = __flags & ios_base::showpoint;
>+
>+      // Insert first digit.
>+      if (!__buf[0])
>+	*__sbuf++ = __digits[0];	// Insert implicit 0.
>+      else
>+	*__sbuf++= __digits[(*__bp++) - '0'];
>+
>+      // Insert decimal point.
>+      // Showpoint always prints point.  Otherwise when prec > 0, there must
>+      // be digits to print or forced padding.
>+      if ((__prec && (*__bp || __pad)) || __showpoint)
>+	*__sbuf++= __dpchar;
>+
>+
>+      // Insert digits after decimal point
>+      for (;__prec && *__bp; __prec--)
>+	*__sbuf++ = __digits[*__bp++ - '0'];
>+      if (__pad)
>+	while (__prec--)
>+	  *__sbuf++ = __digits[0];
>+
>+      // Insert exponent symbol and sign
>+      *__sbuf++ = (__flags & ios_base::uppercase) ? __atoms[__num_base::_S_oE]
>+	: __atoms[__num_base::_S_oe];
>+
>+      __decp--;
>+      if (__decp < 0) 
>+	{
>+	  *__sbuf++ = __atoms[__num_base::_S_ominus];
>+	  __decp = -__decp;
>+	}
>+      else 
>+	*__sbuf++ = __atoms[__num_base::_S_oplus];
>+
>+
>+      // Write the exponent.  Large enough for ieee quad (15 bit exponent)
>+      char __tbuf[6];
>+      char *__tb = __tbuf + 5;
>+      *__tb-- = 0;
>+      while (__decp) 
>+	{
>+	  *__tb-- = '0' + (__decp % 10);
>+	  __decp /= 10;
>+	}
>+
>+      // We need at least 2 digits in exponent.  Pad if it's missing
>+      while (__tb - __tbuf > 2) 
>+	*__tb-- = '0';
>+      while (*++__tb)
>+	*__sbuf++ = __digits[*__tb - '0'];
>+      return __sbuf - __sb;
>+    }
>+
>+  template<typename _Tp>
>+    char* __dtoa(char*, _Tp, int, int, int, bool, int*);
>+
>+  template<typename _Tp>
>+    int __float_to_char_sprintf(char*, _Tp, int, int, ios_base::fmtflags);
>+
>+  // Convert floating point value to string.
>+  //
>+  // Takes a floating point value and converts it to a string using dtoa,
>+  // falling back to sprintf if needed.
>+  //
>+  // buf is the wide char buffer to print into, bufsize is the length of buf.
>+  //
>+  // atoms contains numeric and exponent wide chars to use for the result.
>+  // dpchar is the decimal point char to use.
>+  //
>+  // mode is 0 for %f style, 1 for %e style, and 2 for %g style printing.
>+  // v is the value to print, and exp is the exponent estimate from
>+  // dec_exponent.
>+  //
>+  // prec is the precision to use for printing, and is interpreted according
>+  // to the mode.
>+  //
>+  // flags contains the relevant printing flags: showpoint, showpos, and
>+  // uppercase.
>+  template<typename _CharT, typename _ValueT>
>+    inline int
>+    __float_to_char(_CharT* __buf, int __bufsize, const _CharT* __atoms,
>+		    _CharT __dpchar, int __mode, _ValueT __v, int __exp,
>+		    int __prec, ios_base::fmtflags __flags)
>+    {
>+      char __lbuf[__bufsize];
>+      int __decp;
>+
>+      // Insert sign.
>+      int __len = 0;
>+      if (__v < 0)
>+	{
>+	  *__buf++ = __atoms[__num_base::_S_ominus];
>+	  __len = 1;
>+	}
>+      else if (__flags & ios_base::showpos)
>+	{
>+	  __len = 1;
>+	  *__buf++ = __atoms[__num_base::_S_oplus];
>+	}
>+      
>+      // Dispatch number conversion to the correct function.
>+      if (__mode == 0)
>+	{
>+	  if (!__dtoa(__lbuf, __v, __exp, 3, __prec, false, &__decp))
>+	    __len += __float_to_char_sprintf(__lbuf, __v, __mode, __prec,
>+					     __flags);
>+	  else
>+	    __len +=__float_to_char_fixed(__buf, __atoms, __dpchar,
>+					  &__lbuf[0], __decp, __prec, __flags,
>+					  true);
>+	}
>+      else if (__mode == 1)
>+	{
>+	  if (!__dtoa(__lbuf, __v, __exp, 2, __prec + 1, false, &__decp))
>+	    __len += __float_to_char_sprintf(__lbuf, __v, __mode, __prec,
>+					     __flags);
>+	  else
>+	    __len += __float_to_char_scientific(__buf, __atoms, __dpchar,
>+						&__lbuf[0], __decp, __prec,
>+					      __flags, true);
>+	}
>+      else
>+	{
>+	  if (!__dtoa(__lbuf, __v, __exp, 2, __prec, true, &__decp))
>+	    __len += __float_to_char_sprintf(__lbuf, __v, __mode, __prec,
>+					     __flags);
>+	  else
>+	    {
>+	      if (__decp <= -4 || __decp > __prec)
>+		__len += __float_to_char_scientific(__buf, __atoms, __dpchar,
>+						    &__lbuf[0], __decp,
>+						    __prec - 1, __flags,
>+						    false);
>+	      else
>+		__len += __float_to_char_fixed(__buf, __atoms, __dpchar, 
>+					       &__lbuf[0], __decp,
>+					       __prec - __decp, __flags,
>+					       false);
>+	    }
>+	}
>+      return __len;
>+    }
>+
>+
>+  template<typename _Tp>
>+    int
>+    __dec_exponent(_Tp d, int& special, bool& sign);
>+
>+
>   template<typename _CharT, typename _OutIter>
>     template<typename _ValueT>
>       _OutIter
>@@ -1071,58 +1309,81 @@
> 	// Long enough for the max format spec.
> 	char __fbuf[16];
> 
>-#ifdef _GLIBCXX_USE_C99
>-	// First try a buffer perhaps big enough (most probably sufficient
>-	// for non-ios_base::fixed outputs)
>-	int __cs_size = __max_digits * 3;
>-	char* __cs = static_cast<char*>(__builtin_alloca(__cs_size));
>-
>-	__num_base::_S_format_float(__io, __fbuf, __mod);
>-	__len = std::__convert_from_v(__cs, __cs_size, __fbuf, __v,
>-				      _S_get_c_locale(), __prec);
>-
>-	// If the buffer was not large enough, try again with the correct size.
>-	if (__len >= __cs_size)
>+	int __special;
>+	bool __sign;
>+	_CharT* __ws;
>+	const int __exponent = __dec_exponent(__v, __special, __sign);
>+	if (__special)
> 	  {
>-	    __cs_size = __len + 1;
>-	    __cs = static_cast<char*>(__builtin_alloca(__cs_size));
>-	    __len = std::__convert_from_v(__cs, __cs_size, __fbuf, __v,
>-					  _S_get_c_locale(), __prec);
>+	    const ctype<_CharT>& __ctype = use_facet<ctype<_CharT> >(__loc);
>+	    const char* __cs;
>+	    __ws = static_cast<_CharT*>(__builtin_alloca(sizeof(_CharT) * 4));
>+	    if (__special == 1)
>+	      {
>+		if (__sign)
>+		  __cs = __io.flags() & ios_base::uppercase ? "-INF" : "-inf";
>+		else if (__io.flags() & ios_base::showpos)
>+		  __cs = __io.flags() & ios_base::uppercase ? "+INF" : "+inf";
>+		else
>+		  __cs = __io.flags() & ios_base::uppercase ? "INF" : "inf";
>+	      }
>+	    else
>+	      // special == 2
>+	      __cs = __io.flags() & ios_base::uppercase ? "NAN" : "nan";
>+	    __len = strlen(__cs);
>+	    __ctype.widen(__cs, __cs + __len, __ws);
>+      	  }
>+	else
>+	  {
>+	    // Fixed needs space to hold all digits to left of decimal
>+	    // point, plus up to __prec digits to the right.  The
>+	    // other formats need __prec digits plus sign, decimal
>+	    // point, exponent char and sign, and exponent itself.
>+	    int __decp;
>+	    int __cs_size;
>+	    int __mode;
>+	    const int __maxexp10 = numeric_limits<_ValueT>::max_exponent10;
>+
>+	    if (__io.flags() & ios_base::fixed)
>+	      {
>+		// %f style
>+		// sign + (exp + 1) digits + decpoint + prec digits
>+		__cs_size = (__exponent > 0 ? __exponent : 0) + 3 + __prec;
>+		__mode = 0;
>+	      }
>+	    else if (__io.flags() & ios_base::scientific)
>+	      {
>+		// %e style
>+		// sign + digit + decpt + precision + e + sign + explen
>+		const int __max_exp_len =
>+		  int(__builtin_log10(__maxexp10)) + 1;
>+		__cs_size = __prec + 5 + __max_exp_len;
>+		__mode = 1;
>+	      }
>+	    else
>+	      {
>+		// %g style 
>+		// sign + digit + decpt + prec + e + sign + exp
>+		// When exp >-4 || < prec, use %f style.  We could need 2*prec
>+		// digits.
>+		const int __max_exp_len =
>+		  int(__builtin_log10(__maxexp10)) + 1;
>+		if (!__prec) 
>+		  __prec = 1;
>+		__cs_size = __prec * 2 + 5 + __max_exp_len;
>+		__mode = 2;
>+	      }
>+      
>+	    __ws = static_cast<_CharT*>(__builtin_alloca(sizeof(_CharT)
>+							 * __cs_size));
>+	    __len = __float_to_char(__ws, __cs_size, &__lc->_M_atoms_out[0],
>+				    __lc->_M_decimal_point, __mode, __v,
>+				    __exponent, __prec, __io.flags());
> 	  }
>-#else
>-	// Consider the possibility of long ios_base::fixed outputs
>-	const bool __fixed = __io.flags() & ios_base::fixed;
>-	const int __max_exp = numeric_limits<_ValueT>::max_exponent10;
>-
>-	// The size of the output string is computed as follows.
>-	// ios_base::fixed outputs may need up to __max_exp + 1 chars
>-	// for the integer part + __prec chars for the fractional part
>-	// + 3 chars for sign, decimal point, '\0'. On the other hand,
>-	// for non-fixed outputs __max_digits * 2 + __prec chars are
>-	// largely sufficient.
>-	const int __cs_size = __fixed ? __max_exp + __prec + 4
>-	                              : __max_digits * 2 + __prec;
>-	char* __cs = static_cast<char*>(__builtin_alloca(__cs_size));
>-
>-	__num_base::_S_format_float(__io, __fbuf, __mod);
>-	__len = std::__convert_from_v(__cs, 0, __fbuf, __v,
>-				      _S_get_c_locale(), __prec);
>-#endif
> 
>       // [22.2.2.2.2] Stage 2, convert to char_type, using correct
>       // numpunct.decimal_point() values for '.' and adding grouping.
>-      const ctype<_CharT>& __ctype = use_facet<ctype<_CharT> >(__loc);
>-
>-      _CharT* __ws = static_cast<_CharT*>(__builtin_alloca(sizeof(_CharT)
>-							   * __len));
>-      __ctype.widen(__cs, __cs + __len, __ws);
>-
>-      // Replace decimal point.
>-      const _CharT __cdec = __ctype.widen('.');
>-      const _CharT __dec = __lc->_M_decimal_point;
>-      const _CharT* __p = char_traits<_CharT>::find(__ws, __len, __cdec);
>-      if (__p)
>-	__ws[__p - __ws] = __dec;
>+      // (The conversion happened in __float_to_char(), above.)
> 
>       // Add grouping, if necessary.
>       if (__lc->_M_use_grouping)
>@@ -1131,6 +1392,8 @@
> 	  // number of digits, but no more.
> 	  _CharT* __ws2 = static_cast<_CharT*>(__builtin_alloca(sizeof(_CharT)
> 								* __len * 2));
>+	  const _CharT* __p = 
>+	    char_traits<_CharT>::find(__ws, __len, __lc->_M_decimal_point);
> 	  _M_group_float(__lc->_M_grouping, __lc->_M_grouping_size,
> 			 __lc->_M_thousands_sep, __p, __ws2, __ws, __len);
> 	  __ws = __ws2;
>Index: src/Makefile.am
>===================================================================
>RCS file: /cvs/gcc/gcc/libstdc++-v3/src/Makefile.am,v
>retrieving revision 1.144
>diff -u -r1.144 Makefile.am
>--- src/Makefile.am	16 Apr 2004 19:04:05 -0000	1.144
>+++ src/Makefile.am	31 Aug 2004 00:43:49 -0000
>@@ -102,6 +102,7 @@
> 	ctype.cc \
> 	debug.cc \
> 	debug_list.cc \
>+	floatconv.cc \
> 	functexcept.cc \
> 	globals_locale.cc \
> 	globals_io.cc \
>Index: src/Makefile.in
>===================================================================
>RCS file: /cvs/gcc/gcc/libstdc++-v3/src/Makefile.in,v
>retrieving revision 1.208
>diff -u -r1.208 Makefile.in
>--- src/Makefile.in	23 Aug 2004 10:18:30 -0000	1.208
>+++ src/Makefile.in	31 Aug 2004 00:43:50 -0000
>@@ -65,7 +65,7 @@
> 	numeric_members.lo time_members.lo
> am__objects_2 = basic_file.lo c++locale.lo
> am__objects_3 = allocator.lo codecvt.lo complex_io.lo ctype.lo \
>-	debug.lo debug_list.lo functexcept.lo globals_locale.lo \
>+	debug.lo debug_list.lo floatconv.lo functexcept.lo globals_locale.lo \
> 	globals_io.lo ios.lo ios_failure.lo ios_init.lo ios_locale.lo \
> 	limits.lo list.lo locale.lo locale_init.lo locale_facets.lo \
> 	localename.lo stdexcept.lo strstream.lo tree.lo \
>@@ -309,6 +309,7 @@
> 	ctype.cc \
> 	debug.cc \
> 	debug_list.cc \
>+	floatconv.cc \
> 	functexcept.cc \
> 	globals_locale.cc \
> 	globals_io.cc \
>Index: src/locale-inst.cc
>===================================================================
>RCS file: /cvs/gcc/gcc/libstdc++-v3/src/locale-inst.cc,v
>retrieving revision 1.50
>diff -u -r1.50 locale-inst.cc
>--- src/locale-inst.cc	27 Feb 2004 10:12:01 -0000	1.50
>+++ src/locale-inst.cc	31 Aug 2004 00:43:50 -0000
>@@ -32,6 +32,12 @@
> // ISO C++ 14882: 22.1  Locales
> //
> 
>+// This file instantiates locale-related classes and functions for char and
>+// wchar_t.  Use C in place of the template char type.  wlocale-inst.cc
>+// changes the definition of C to reinstantiate everything with wchar_t.  For
>+// locale-related functions that don't depend on the char type, see
>+// locale-misc-inst.cc.
>+
> #include <locale>
> 
> // Instantiation configuration.
>@@ -313,5 +319,17 @@
>     __int_to_char(C*, unsigned long long, const C*, 
> 		  ios_base::fmtflags, bool);
> #endif
>+
>+
>+  template
>+    int
>+    __float_to_char_fixed(C*, const C*, C, const char*, int, int,
>+			  ios_base::fmtflags, bool);
>+
>+  template
>+    int
>+    __float_to_char_scientific(C*, const C*, C, const char*, int, int,
>+			       ios_base::fmtflags, bool);
>+
> } // namespace std
>--- /dev/null	2004-08-19 18:00:15.000000000 -0400
>+++ src/floatconv.cc	2004-08-30 21:31:17.000000000 -0400
>@@ -0,0 +1,704 @@
>+// Numeric printing support -*- C++ -*-
>+
>+// Copyright (C) 1993, 1994, 2004
>+// Free Software Foundation, Inc.
>+//
>+// This file is part of the GNU ISO C++ Library.  This library is free
>+// software; you can redistribute it and/or modify it under the
>+// terms of the GNU General Public License as published by the
>+// Free Software Foundation; either version 2, or (at your option)
>+// any later version.
>+
>+// This library is distributed in the hope that it will be useful,
>+// but WITHOUT ANY WARRANTY; without even the implied warranty of
>+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>+// GNU General Public License for more details.
>+
>+// You should have received a copy of the GNU General Public License along
>+// with this library; see the file COPYING.  If not, write to the Free
>+// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
>+// USA.
>+
>+// As a special exception, you may use this file as part of a free software
>+// library without restriction.  Specifically, if other files instantiate
>+// templates or use macros or inline functions from this file, or you compile
>+// this file and link it with other files to produce an executable, this
>+// file does not by itself cause the resulting executable to be covered by
>+// the GNU General Public License.  This exception does not however
>+// invalidate any other reasons why the executable file might be covered by
>+// the GNU General Public License.
>+
>+// Original copyright below
>+
>+/* 
>+Copyright (C) 1993, 1994 Free Software Foundation
>+
>+This file is part of the GNU IO Library.  This library is free
>+software; you can redistribute it and/or modify it under the
>+terms of the GNU General Public License as published by the
>+Free Software Foundation; either version 2, or (at your option)
>+any later version.
>+
>+This library is distributed in the hope that it will be useful,
>+but WITHOUT ANY WARRANTY; without even the implied warranty of
>+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
>+GNU General Public License for more details.
>+
>+You should have received a copy of the GNU General Public License
>+along with this library; see the file COPYING.  If not, write to the Free
>+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
>+
>+As a special exception, if you link this library with files
>+compiled with a GNU compiler to produce an executable, this does not cause
>+the resulting executable to be covered by the GNU General Public License.
>+This exception does not however invalidate any other reasons why
>+the executable file might be covered by the GNU General Public License. */
>+
>+// #include <libioP.h>
>+/****************************************************************
>+ *
>+ * The author of this software is David M. Gay.
>+ *
>+ * Copyright (c) 1991 by AT&T.
>+ *
>+ * Permission to use, copy, modify, and distribute this software for any
>+ * purpose without fee is hereby granted, provided that this entire notice
>+ * is included in all copies of any software which is or includes a copy
>+ * or modification of this software and in all copies of the supporting
>+ * documentation for such software.
>+ *
>+ * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
>+ * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
>+ * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
>+ * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
>+ *
>+ ***************************************************************/
>+
>+/* Some cleaning up by Per Bothner, bothner@cygnus.com, 1992, 1993.
>+   Re-written to not need static variables
>+   (except result, result_k, HIWORD, LOWORD). */
>+
>+/* Note that the checking of _DOUBLE_IS_32BITS is for use with the
>+   cross targets that employ the newlib ieeefp.h header.  -- brendan */
>+
>+/* Please send bug reports to
>+        David M. Gay
>+        AT&T Bell Laboratories, Room 2C-463
>+        600 Mountain Avenue
>+        Murray Hill, NJ 07974-2070
>+        U.S.A.
>+        dmg@research.att.com or research!dmg
>+ */
>+
>+/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
>+ *
>+ * This strtod returns a nearest machine number to the input decimal
>+ * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
>+ * broken by the IEEE round-even rule.  Otherwise ties are broken by
>+ * biased rounding (add half and chop).
>+ *
>+ * Inspired loosely by William D. Clinger's paper "How to Read Floating
>+ * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
>+ *
>+ * Modifications:
>+ *
>+ *      1. We only require IEEE, IBM, or VAX double-precision
>+ *              arithmetic (not IEEE double-extended).
>+ *      2. We get by with floating-point arithmetic in a case that
>+ *              Clinger missed -- when we're computing d * 10^n
>+ *              for a small integer d and the integer n is not too
>+ *              much larger than 22 (the maximum integer k for which
>+ *              we can represent 10^k exactly), we may be able to
>+ *              compute (d*10^k) * 10^(e-k) with just one roundoff.
>+ *      3. Rather than a bit-at-a-time adjustment of the binary
>+ *              result in the hard case, we use floating-point
>+ *              arithmetic to determine the adjustment to within
>+ *              one bit; only in really hard cases do we need to
>+ *              compute a second residual.
>+ *      4. Because of 3., we don't need a large table of powers of 10
>+ *              for ten-to-e (just some small tables, e.g. of 10^k
>+ *              for 0 <= k <= 22).
>+ */
>+
>+/*
>+ * #define IEEE for IEEE floating-point arithmetic
>+ * #define IBM for IBM mainframe-style floating-point arithmetic.
>+ * #define VAX for VAX-style floating-point arithmetic.
>+ * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
>+ * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
>+ */
>+
>+#include <clocale>
>+#include <cstring>
>+#include <cmath>		// log, floor, frexp
>+#include <limits>
>+#include <ios>
>+
>+namespace std
>+{
>+
>+// Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16
>+// Quick_max = floor((mant_size-1)*log(FLT_RADIX)/log(10) - 1)
>+// Int_max = floor(mant_size*log(FLT_RADIX)/log(10) - 1)
>+template<typename _T>
>+struct __fptype
>+{
>+  typedef std::numeric_limits<_T> nld;
>+
>+  static const int mant_size = nld::digits;
>+  static const int Bletch;
>+  static const int Quick_max;
>+  static const int Int_max;
>+};
>+
>+template<typename _T>
>+const int __fptype<_T>::Bletch
>+  = (1 << __builtin_ffs(__fptype::nld::max_exponent10)) / 16;
>+
>+template<typename _T>
>+const int __fptype<_T>::Quick_max
>+  = int(floor((__fptype::nld::digits-1) *
>+	      log(double(__fptype::nld::radix))/log(10) - 1));
>+
>+// For secondary fast codepath
>+template<typename _T>
>+const int __fptype<_T>::Int_max
>+  = int(floor(__fptype::nld::digits * log(__fptype::nld::radix)/log(10) - 1));
>+
>+template
>+struct __fptype<double>;
>+
>+template
>+struct __fptype<long double>;
>+
>+
>+// Ten_pmax = floor(mant_size*log(2)/log(5))
>+// IEEE & IBM double is 22
>+// VAX double is 24
>+// IEEE quad is 48
>+// This is the number of elements in tens[] to use.
>+static const int __Ten_pmax = int(__fptype<double>::mant_size * log(2)/log(5));
>+
>+// Enough entries to handle long double with 112-bit mantissa
>+static const double
>+__tens[] = {
>+                1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
>+                1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
>+                1e20, 1e21, 1e22, 1e23, 1e24, 1e25, 1e26, 1e27, 1e28, 1e29,
>+                1e30, 1e31, 1e32, 1e33, 1e34, 1e35, 1e36, 1e37, 1e38, 1e39,
>+                1e40, 1e41, 1e42, 1e43, 1e44, 1e45, 1e46, 1e47, 1e48, 1e49
>+};
>+static const double __bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256, 1e512,
>+				  1e1024, 1e2048, 1e4096 };
>+// 5 for ieee double
>+// 3 for s370 double
>+// 2 for vax double
>+static const int __n_bigtens =
>+  int(log(std::numeric_limits<double>::max_exponent10)/log(2.0)) - 3;
>+
>+static const long double __bigtensl[] =
>+  { 1e16L, 1e32L, 1e64L, 1e128L, 1e256L, 1e512L, 1e1024L, 1e2048L, 1e4096L,
>+    1e8192L, 1e16384L };
>+
>+static const int __n_bigtensl =
>+  int(log(std::numeric_limits<long double>::max_exponent10)/log(2.0)) - 3;
>+
>+
>+template<typename _T>
>+_T __bigten(int i);
>+
>+template<>
>+static inline double
>+__bigten<double>(int i) { return __bigtens[i]; }
>+
>+template<>
>+static inline long double
>+__bigten<long double>(int i) { return __bigtensl[i]; }
>+
>+template<typename _T>
>+_T __maxbigten();
>+
>+template<>
>+static inline double __maxbigten<double>()
>+{ return __bigtens[__n_bigtens-1]; }
>+
>+template<>
>+static inline long double __maxbigten<long double>()
>+{ return __bigtensl[__n_bigtensl-1]; }
>+
>+
>+
>+// dec_exponent performs the first portion of the conversion to a string.  It
>+// computes an upper bound on the decimal exponent that can be used to size a
>+// buffer to print into.  It must be passed back into dtoa.
>+//
>+// This will only be used for double and long double
>+template<typename _T>
>+int
>+__dec_exponent(_T d, int &special, bool &sign)
>+{
>+  /*     Arguments ndigits, decpt, sign are similar to those
>+	 of ecvt and fcvt; trailing zeros are suppressed from
>+	 the returned string.  If not null, *rve is set to point
>+	 to the end of the return value.  If d is +-Infinity or NaN,
>+	 then *decpt is set to 9999.
>+
>+	 special returns true if the number is nan, infinity, or 0.
>+	 type contains -1 for nan, 0 for 0, and 1 for infinity.  This allows
>+	 the calling routine to handle the special cases.
>+  */
>+
>+  int i;
>+  int k;			// Holds base 10 exponent
>+
>+  _T d2, ds;
>+
>+  // Handle NAN.
>+  if (__builtin_isnan(d))
>+    {
>+      special = 2;
>+      return 3;
>+    }
>+
>+  // Save sign value into user var
>+  if (__builtin_signbit(d)) {
>+    // set sign for everything, including 0's and NaNs
>+    sign = true;
>+    d = -d;		// Clear sign bit.
>+  }
>+  else
>+    sign = false;
>+
>+  // Handle IEEE and VAX infinity and nan.  Apparently IBM has none.
>+  if (__builtin_isinf(d))
>+    {
>+      special = 1;
>+      return 8;
>+    }
>+
>+  special = 0;
>+
>+  if (d == 0.0) return 0;
>+
>+  // Extract base-2 exponent and normalized fraction.
>+  d2 = frexp(d, &i);
>+  
>+  // Compute integer base 10 exponent into k
>+
>+  /* log(x)       ~=~ log(1.5) + (x-1.5)/1.5
>+   * log10(x)      =  log(x) / log(10)
>+   *              ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
>+   * log10(d) = i*log(2)/log(10) + log10(d2)
>+   *
>+   * This suggests computing an approximation k to log10(d) by
>+   *
>+   * k = i*0.301029995663981
>+   *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
>+   *
>+   * We want k to be too large rather than too small.
>+   * The error in the first-order Taylor series approximation
>+   * is in our favor, so we just round up the constant enough
>+   * to compensate for any error in the multiplication of
>+   * (i) by 0.301029995663981; since |i| <= 1077,
>+   * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
>+   * adding 1e-13 to the constant term more than suffices.
>+   * Hence we adjust the constant term to 0.1760912590558.
>+   * (We could get a more accurate k by invoking log10,
>+   *  but this is probably not worthwhile.)
>+   *
>+   * (From Richard Henderson)
>+   * The comment should be updated for long double.  I guess the
>+   * overestimate is still ok, since while |i| <= 16881 is much
>+   * larger, 2^-64 and 2^-113 are much smaller.
>+   */
>+  ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
>+  k = (int)ds;
>+  if (ds < 0. && ds != k)
>+    k--;    // want k = floor(ds)
>+  return k;
>+}
>+
>+template
>+int
>+__dec_exponent<double>(double d, int &special, bool &sign);
>+
>+template
>+int
>+__dec_exponent<long double>(long double d, int &special, bool &sign);
>+
>+/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
>+ *
>+ * Inspired by "How to Print Floating-Point Numbers Accurately" by
>+ * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
>+ *
>+ * Modifications:
>+ *      1. Rather than iterating, we use a simple numeric overestimate
>+ *         to determine k = floor(log10(d)).  We scale relevant
>+ *         quantities using O(log2(k)) rather than O(k) multiplications.
>+ *      2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
>+ *         try to generate digits strictly left to right.  Instead, we
>+ *         compute with fewer bits and propagate the carry if necessary
>+ *         when rounding the final digit up.  This is often faster.
>+ *      3. Under the assumption that input will be rounded nearest,
>+ *         mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
>+ *         That is, we allow equality in stopping tests when the
>+ *         round-nearest rule will give the same floating-point value
>+ *         as would satisfaction of the stopping test with strict
>+ *         inequality.
>+ *      4. We remove common factors of powers of 2 from relevant
>+ *         quantities.
>+ *      5. When converting floating-point integers less than 1e16,
>+ *         we use floating-point arithmetic rather than resorting
>+ *         to multiple-precision integers.
>+ *      6. When asked to produce fewer than 15 digits, we first try
>+ *         to get by with floating-point arithmetic; we resort to
>+ *         multiple-precision integer arithmetic only if we cannot
>+ *         guarantee that the floating-point calculation has given
>+ *         the correctly rounded result.  For k requested digits and
>+ *         "uniformly" distributed input, the probability is
>+ *         something like 10^(k-15) that we must resort to the long
>+ *         calculation.
>+ */
>+
>+// Slow method requiring big int math has been removed, returning 0.  This is
>+// because printf was about 30% faster by comparison.
>+//
>+// buf is a buffer whose size was computed using dec_exponent(), so as to be
>+// large enough to hold the number if needed.
>+//
>+// exponent is the value computed by dec_exponent
>+//
>+// The mode is as described below, but only 2 and 3 are still used.
>+//
>+// Mode 2 ndigits gives the total significant digits to print.
>+// Mode 3 ndigits gives the total digits after the decimal point to print.
>+//
>+// If trim is true, trailing zeros will be trimmed off.
>+//
>+// decpt returns the actual decimal point location.  It will often be
>+// exponent, but may be adjusted in some cases.
>+//
>+// The return value is 0 if the routine gave up, otherwise it is buf.
>+
>+template<typename _T>
>+char *
>+__dtoa(char* buf, _T d, int exponent, int mode, int ndigits, bool trim,
>+       int *decpt)
>+{
>+  /*     Arguments ndigits, decpt, sign are similar to those of ecvt and fcvt;
>+	 trailing zeros are suppressed from the returned string.  If d is
>+	 +-Infinity or NaN, then *decpt is set to 9999.
>+
>+	 mode:
>+	 0 ==> shortest string that yields d when read in
>+	       and rounded to nearest.
>+	 1 ==> like 0, but with Steele & White stopping rule;
>+	       e.g. with IEEE P754 arithmetic , mode 0 gives
>+	       1e23 whereas mode 1 gives 9.999999999999999e22.
>+	 2 ==> max(1,ndigits) significant digits.  This gives a
>+	       return value similar to that of ecvt, except
>+	       that trailing zeros are suppressed.
>+	 3 ==> through ndigits past the decimal point.  This
>+	       gives a return value similar to that from fcvt,
>+	       except that trailing zeros are suppressed, and
>+	       ndigits can be negative.
>+
>+	 4-9 should give the same return values as 2-3, i.e.,
>+	       4 <= mode <= 9 ==> same return as mode
>+	       2 + (mode & 1).  These modes are mainly for
>+	       debugging; often they run slower but sometimes
>+	       faster than modes 2-3.
>+	 4,5,8,9 ==> left-to-right digit generation.
>+	 6-9 ==> don't try fast floating-point estimate
>+	 	(if applicable).
>+
>+	 Values of mode other than 0-9 are treated as mode 0.
>+
>+	 Sufficient space is allocated to the return value
>+	 to hold the suppressed trailing zeros.
>+  */
>+
>+  int i, ilim=0, ilim1=0;
>+  int k;			// Holds base 10 exponent
>+  int k_check;
>+  int L;
>+  _T ds, eps;
>+  char *s, *s0, *stmp;
>+
>+  // Handle 0 for all cases.
>+  if (d == 0.0) {
>+    *decpt = 1;
>+    buf[0] = '0';
>+    buf[1] = 0;
>+    return buf;
>+  }
>+
>+  // Clear the sign.
>+  d = fabs(d);
>+
>+  // Extract base-2 exponent.
>+  frexp(d, &i);
>+
>+  // k is an upper bound on the integer base 10 exponent of d.  We may need to
>+  // tweak it a little.
>+  k = exponent;
>+  k_check = 1;
>+  if (k >= 0 && k < __Ten_pmax) {
>+    if (d < __tens[k])
>+      k--;
>+    k_check = 0;
>+  }
>+
>+  if (mode == 2) {
>+    if (ndigits <= 0)
>+      ndigits = 1;
>+    ilim = ilim1 = i = ndigits;
>+  }
>+  else {
>+    i = ndigits + k + 1;
>+    ilim = i;
>+    ilim1 = i - 1;
>+    if (i <= 0)
>+      i = 1;
>+  }
>+  // i is now an upper bound of the number of digits to generate.
>+
>+  s = s0 = buf;
>+
>+  if (ilim >= 0 && ilim <= __fptype<_T>::Quick_max) {
>+
>+    // Try to get by with floating-point arithmetic.
>+
>+    _T d2 = d;
>+    int i = 0;
>+    int k0 = k;
>+    int ilim0 = ilim;
>+    int ieps = 2; // conservative
>+    int j1;
>+    int j;
>+
>+    if (k > 0) {
>+      ds = __tens[k&0xf];		// 1e0 -> 1e15
>+      j = k >> 4;
>+      if (j & __fptype<_T>::Bletch) {
>+	// I've never seen this code hit in .5 million random numbers
>+	// prevent overflows
>+	// n_bigtens
>+	j &= __fptype<_T>::Bletch - 1;
>+	d /= __maxbigten<_T>();
>+	ieps++;
>+      }
>+      for(; j; j >>= 1, i++)
>+	if (j & 1) {
>+	  ieps++;
>+	  ds *= __bigten<_T>(i);
>+	}
>+      d /= ds;
>+    }
>+    else if ((j1 = -k)) {
>+      d *= __tens[j1 & 0xf];
>+      for(j = j1 >> 4; j; j >>= 1, i++)
>+	if (j & 1) {
>+	  ieps++;
>+	  d *= __bigten<_T>(i);
>+	}
>+    }
>+    if (k_check && d < 1. && ilim > 0) {
>+      if (ilim1 <= 0)
>+	goto fast_failed;
>+      ilim = ilim1;
>+      k--;
>+      d *= 10.;
>+      ieps++;
>+    }
>+    eps = ieps*d + 7.;
>+    static const _T e2m = ldexp(2, -(__fptype<_T>::mant_size-1));
>+    eps *= e2m;
>+    if (ilim == 0) {
>+      d -= 5.;
>+      if (d > eps)
>+	goto one_digit;
>+      if (d < -eps)
>+	goto no_digits;
>+      goto fast_failed;
>+    }
>+
>+    // Generate ilim digits, then fix them up.
>+    eps *= __tens[ilim-1];
>+    for(i = 1;; i++, d *= 10.) {
>+      L = (int)d;
>+      d -= L;
>+      *s++ = '0' + (int)L;
>+      if (i == ilim) {
>+	if (d > 0.5 + eps)
>+	  goto bump_up;
>+	else if (d < 0.5 - eps) {
>+	  // Strips off trailing 0's
>+	  if (trim) {
>+	    while(*--s == '0');
>+	    s++;
>+	  }
>+	  goto ret;
>+	}
>+	// else fast approach failed
>+	break;
>+      }
>+    }
>+  fast_failed:
>+    s = s0;
>+    d = d2;
>+    k = k0;
>+    ilim = ilim0;
>+  }
>+
>+  // Do we have a "small" integer?
>+
>+  if (k <= __fptype<_T>::Int_max && floor(d) == d) {
>+    // Yes.
>+    ds = __tens[k];
>+    if (ndigits < 0 && ilim <= 0) {
>+      if (ilim < 0 || d <= 5*ds)
>+	goto no_digits;
>+      goto one_digit;
>+    }
>+    for(i = 1;; i++) {
>+      L = (int)(d / ds);
>+      d -= L*ds;
>+
>+      if (std::numeric_limits<_T>::round_style
>+	  == std::round_toward_infinity ||
>+	  std::numeric_limits<_T>::round_style
>+	  == std::round_toward_neg_infinity)
>+	{
>+	// If FLT_ROUNDS == 2, L will usually be high by 1
>+	if (d < 0) {
>+	  L--;
>+	  d += ds;
>+	}
>+      }
>+
>+      *s++ = '0' + (int)L;
>+      if (i == ilim) {
>+	d += d;
>+	if (d > ds || (d == ds && L & 1)) {
>+	bump_up:
>+	  stmp = s;
>+	  while( *--s == '9' ) {
>+	    *s='0';
>+	    if (s == s0) {
>+	      k++;
>+	      *s = '0';
>+	      break;
>+	    }
>+	  }
>+	  ++*s++;
>+	  if (!trim) s = stmp;
>+	}
>+	break;
>+      }
>+      if (!(d *= 10.))
>+	break;
>+    }
>+    goto ret;
>+  }
>+
>+  return 0;			// XXX bail out to avoid mp int calcs
>+
>+  /* no digits, fcvt style */
>+ no_digits:
>+  k = -1 - ndigits;
>+  goto ret;
>+
>+ one_digit:
>+  *s++ = '1';
>+  k++;
>+  goto ret;
>+
>+ ret:
>+  *s = 0;
>+  *decpt = k + 1;
>+  return s0;
>+}
>+
>+template
>+char *
>+__dtoa<double>(char* buf, double d, int exponent, int mode, int ndigits,
>+	       bool trim, int *decpt);
>+
>+template
>+char *
>+__dtoa<long double>(char* buf, long double d, int exponent, int mode,
>+		    int ndigits, bool trim, int *decpt);
>+
>+
>+// Call to detect if we're working on double or long double.
>+template<typename _T>
>+static inline bool __is_long();
>+
>+template<>
>+static inline bool __is_long<double>() { return false; }
>+
>+template<>
>+static inline bool __is_long<long double>() { return true; }
>+
>+// The last three are in flags.
>+template<typename _T>
>+int
>+__float_to_char_sprintf(char* buf, _T d, int mode, int prec,
>+			ios_base::fmtflags flags)
>+{
>+  static const char lmodes[] = "feg";
>+  static const char umodes[] = "FEG";
>+
>+  bool point = flags & ios_base::showpoint;
>+  bool pos = flags & ios_base::showpos;
>+  bool up = flags & ios_base::uppercase;
>+
>+  // Build the format string
>+  char fmt[8] = "%";
>+  char* fp = &fmt[1];
>+  if (pos) *fp++ = '+';
>+  if (point) *fp++ = '#';
>+  *fp++ = '.';
>+  *fp++ = '*';
>+
>+  if (__is_long<_T>())
>+    *fp++ = 'L';
>+
>+  if (up)
>+    *fp++ = umodes[mode];
>+  else
>+    *fp++ = lmodes[mode];
>+
>+  *fp=0;
>+
>+  // This is safe since the buffer will be big enough.
>+  sprintf(&buf[0], fmt, prec, d);
>+
>+  // Standard says we need to print as if in the C locale.  The only
>+  // difference that may show up is that the point char can be different.
>+  // Change it to '.'.
>+  char dp = localeconv()->decimal_point[0];
>+  char* bp = buf;
>+  for (; *bp; bp++)
>+    if (*bp == dp) {
>+      *bp = '.';
>+      break;
>+    }
>+  // Return the full length of the string
>+  int len = bp - buf;
>+  len += strlen(bp);
>+  return len;
>+}
>+
>+template
>+int
>+__float_to_char_sprintf<double>(char* buf, double d, int mode, int prec,
>+				ios_base::fmtflags flags);
>+
>+template
>+int
>+__float_to_char_sprintf<long double>(char* buf, long double d, int mode,
>+				     int prec, ios_base::fmtflags flags);
>+
>+}
>
>
>  
>




More information about the Libstdc++ mailing list