summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
Diffstat (limited to 'generic/tclStrToD.c')
-rwxr-xr-xgeneric/tclStrToD.c1361
1 files changed, 1361 insertions, 0 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
new file mode 100755
index 0000000..2d622a7
--- /dev/null
+++ b/generic/tclStrToD.c
@@ -0,0 +1,1361 @@
+/*
+ *----------------------------------------------------------------------
+ *
+ * tclStrToD.c --
+ *
+ * This file contains a TclStrToD procedure that handles conversion
+ * of string to double, with correct rounding even where extended
+ * precision is needed to achieve that. It also contains a
+ * TclDoubleDigits procedure that handles conversion of double
+ * to string (at least the significand), and several utility functions
+ * for interconverting 'double' and the integer types.
+ *
+ * Copyright (c) 2005 by Kevin B. Kenny. All rights reserved.
+ *
+ * See the file "license.terms" for information on usage and redistribution
+ * of this file, and for a DISCLAIMER OF ALL WARRANTIES.
+ *
+ * RCS: @(#) $Id: tclStrToD.c,v 1.2 2005/05/10 18:34:49 kennykb Exp $
+ *
+ *----------------------------------------------------------------------
+ */
+
+#include <tclInt.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <ctype.h>
+#include <tommath.h>
+
+/*
+ * The stuff below is a bit of a hack so that this file can be used in
+ * environments that include no UNIX, i.e. no errno: just arrange to use
+ * the errno from tclExecute.c here.
+ */
+
+#ifdef TCL_GENERIC_ONLY
+#define NO_ERRNO_H
+#endif
+
+#ifdef NO_ERRNO_H
+extern int errno; /* Use errno from tclExecute.c. */
+#define ERANGE 34
+#endif
+
+#if ( FLT_RADIX == 2 ) && ( DBL_MANT_DIG == 53 ) && ( DBL_MAX_EXP == 1024 )
+#define IEEE_FLOATING_POINT
+#endif
+
+/*
+ * gcc on x86 needs access to rounding controls. It is tempting to
+ * include fpu_control.h, but that file exists only on Linux; it is
+ * missing on Cygwin and MinGW.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
+#define _FPU_GETCW(cw) __asm__ ("fnstcw %0" : "=m" (*&cw))
+#define _FPU_SETCW(cw) __asm__ ("fldcw %0" : : "m" (*&cw))
+#endif
+
+/* The powers of ten that can be represented exactly as IEEE754 doubles. */
+
+#define MAXPOW 22
+static double pow10 [MAXPOW+1];
+
+static int mmaxpow; /* Largest power of ten that can be
+ * represented exactly in a 'double'. */
+
+/* Inexact higher powers of ten */
+
+static CONST double pow_10_2_n [] = {
+ 1.0,
+ 100.0,
+ 10000.0,
+ 1.0e+8,
+ 1.0e+16,
+ 1.0e+32,
+ 1.0e+64,
+ 1.0e+128,
+ 1.0e+256
+};
+
+/* Logarithm of the floating point radix. */
+
+static int log2FLT_RADIX;
+
+/* Number of bits in a double's significand */
+
+static int mantBits;
+
+/* Table of powers of 5**(2**n), up to 5**256 */
+
+static mp_int pow5[9];
+
+/* The smallest representable double */
+
+static double tiny;
+
+/* The maximum number of digits to the left of the decimal point of a
+ * double. */
+
+static int maxDigits;
+
+/* The maximum number of digits to the right of the decimal point in a
+ * double. */
+
+static int minDigits;
+
+/* Number of mp_digit's needed to hold the significand of a double */
+
+static int mantDIGIT;
+
+/* Static functions defined in this file */
+
+static double RefineResult _ANSI_ARGS_((double approx, CONST char* start,
+ int nDigits, long exponent));
+static double ParseNaN _ANSI_ARGS_(( int signum, CONST char** end ));
+static double SafeLdExp _ANSI_ARGS_(( double fraction, int exponent ));
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclStrToD --
+ *
+ * Scans a double from a string.
+ *
+ * Results:
+ * Returns the scanned number. In the case of underflow, returns
+ * an appropriately signed zero; in the case of overflow, returns
+ * an appropriately signed HUGE_VAL.
+ *
+ * Side effects:
+ * Stores a pointer to the end of the scanned number in '*endPtr',
+ * if endPtr is not NULL. If '*endPtr' is equal to 's' on return from
+ * this function, it indicates that the input string could not be
+ * recognized as a number.
+ * In the case of underflow or overflow, 'errno' is set to ERANGE.
+ *
+ *------------------------------------------------------------------------
+ */
+
+double
+TclStrToD( CONST char* s,
+ /* String to scan */
+ CONST char ** endPtr )
+ /* Pointer to the end of the scanned number */
+{
+
+ CONST char* p = s;
+ CONST char* startOfSignificand = NULL;
+ /* Start of the significand in the
+ * string */
+ int signum = 0; /* Sign of the significand */
+ double exactSignificand = 0.0;
+ /* Significand, represented exactly
+ * as a floating-point number */
+ int seenDigit = 0; /* Flag == 1 if a digit has been seen */
+ int nSigDigs = 0; /* Number of significant digits presented */
+ int nDigitsAfterDp = 0; /* Number of digits after the decimal point */
+ int nTrailZero = 0; /* Number of trailing zeros in the
+ * significand */
+ long exponent = 0; /* Exponent */
+ int seenDp = 0; /* Flag == 1 if decimal point has been seen */
+
+ char c; /* One character extracted from the input */
+
+ /*
+ * v must be 'volatile double' on gc-ix86 to force correct rounding
+ * to IEEE double and not Intel double-extended.
+ */
+
+ volatile double v; /* Scanned value */
+ int machexp; /* Exponent of the machine rep of the
+ * scanned value */
+ int expt2; /* Exponent for computing first
+ * approximation to the true value */
+ int i, j;
+
+ /*
+ * With gcc on x86, the floating point rounding mode is double-extended.
+ * This causes the result of double-precision calculations to be rounded
+ * twice: once to the precision of double-extended and then again to the
+ * precision of double. Double-rounding introduces gratuitous errors of
+ * 1 ulp, so we need to change rounding mode to 53-bits.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+ fpu_control_t roundTo53Bits = 0x027f;
+ fpu_control_t oldRoundingMode;
+ _FPU_GETCW( oldRoundingMode );
+ _FPU_SETCW( roundTo53Bits );
+#endif
+
+ /* Discard leading whitespace */
+
+ while ( isspace( *p ) ) {
+ ++p;
+ }
+
+ /* Determine the sign of the significand */
+
+ switch( *p ) {
+ case '-':
+ signum = 1;
+ /* FALLTHROUGH */
+ case '+':
+ ++p;
+ }
+
+ /* Discard leading zeroes */
+
+ while ( *p == '0' ) {
+ seenDigit = 1;
+ ++p;
+ }
+
+ /*
+ * Scan digits from the significand. Simultaneously, keep track
+ * of the number of digits after the decimal point. Maintain
+ * a pointer to the start of the significand. Keep "exactSignificand"
+ * equal to the conversion of the DBL_DIG most significant digits.
+ */
+
+ for ( ; ; ) {
+ c = *p;
+ if ( c == '.' && !seenDp ) {
+ seenDp = 1;
+ ++p;
+ } else if ( isdigit(c) ) {
+ if ( c == '0' ) {
+ if ( startOfSignificand != NULL ) {
+ ++nTrailZero;
+ }
+ } else {
+ if ( startOfSignificand == NULL ) {
+ startOfSignificand = p;
+ } else if ( nTrailZero ) {
+ if ( nTrailZero + nSigDigs < DBL_DIG ) {
+ exactSignificand *= pow10[ nTrailZero ];
+ } else if ( nSigDigs < DBL_DIG ) {
+ exactSignificand *= pow10[ DBL_DIG - nSigDigs ];
+ }
+ nSigDigs += nTrailZero;
+ }
+ if ( nSigDigs < DBL_DIG ) {
+ exactSignificand = 10. * exactSignificand + (c - '0');
+ }
+ ++nSigDigs;
+ nTrailZero = 0;
+ }
+ if ( seenDp ) {
+ ++nDigitsAfterDp;
+ }
+ seenDigit = 1;
+ ++p;
+ } else {
+ break;
+ }
+ }
+
+ /*
+ * At this point, we've scanned the significand, and p points
+ * to the character beyond it. "startOfSignificand" is the first
+ * non-zero character in the significand. "nSigDigs" is the number
+ * of significant digits of the significand, not including any
+ * trailing zeroes. "exactSignificand" is a floating point number
+ * that represents, without loss of precision, the first
+ * min(DBL_DIG,n) digits of the significand. "nDigitsAfterDp"
+ * is the number of digits after the decimal point, again excluding
+ * trailing zeroes.
+ *
+ * Now scan 'E' format
+ */
+
+ exponent = 0;
+ if ( seenDigit && ( *p == 'e' || *p == 'E' ) ) {
+ CONST char* stringSave = p;
+ ++p;
+ c = *p;
+ if ( isdigit( c ) || c == '+' || c == '-' ) {
+ errno = 0;
+ exponent = strtol( p, (char**)&p, 10 );
+ if ( errno == ERANGE ) {
+ if ( exponent > 0 ) {
+ v = HUGE_VAL;
+ } else {
+ v = 0.0;
+ }
+ *endPtr = p;
+ goto returnValue;
+ }
+ }
+ if ( p == stringSave + 1 ) {
+ p = stringSave;
+ exponent = 0;
+ }
+ }
+ exponent = exponent + nTrailZero - nDigitsAfterDp;
+
+ /*
+ * If we come here with no significant digits, we might still be
+ * looking at Inf or NaN. Go parse them.
+ */
+
+ if ( !seenDigit ) {
+
+ /* Test for Inf */
+
+ if ( c == 'I' || c == 'i' ) {
+
+ if ( ( p[1] == 'N' || p[1] == 'n' )
+ && ( p[2] == 'F' || p[2] == 'f' ) ) {
+ p += 3;
+ if ( ( p[0] == 'I' || p[0] == 'i' )
+ && ( p[1] == 'N' || p[1] == 'n' )
+ && ( p[2] == 'I' || p[2] == 'i' )
+ && ( p[3] == 'T' || p[3] == 't' )
+ && ( p[4] == 'Y' || p[1] == 'y' ) ) {
+ p += 5;
+ }
+ errno = ERANGE;
+ v = HUGE_VAL;
+ if ( endPtr != NULL ) {
+ *endPtr = p;
+ }
+ goto returnValue;
+ }
+
+
+#ifdef IEEE_FLOATING_POINT
+
+ /* IEEE floating point supports NaN */
+
+ } else if ( (c == 'N' || c == 'n' )
+ && ( sizeof(Tcl_WideUInt) == sizeof( double ) ) ) {
+
+ if ( ( p[1] == 'A' || p[1] == 'a' )
+ && ( p[2] == 'N' || p[2] == 'n' ) ) {
+ p += 3;
+
+ if ( endPtr != NULL ) {
+ *endPtr = p;
+ }
+
+ /* Restore FPU mode word */
+
+#if defined(__GNUC__) && defined(__i386)
+ _FPU_SETCW( oldRoundingMode );
+#endif
+ return ParseNaN( signum, endPtr );
+
+ }
+#endif
+
+ }
+
+ goto error;
+ }
+
+ /*
+ * We've successfully scanned; update the end-of-element pointer.
+ */
+
+ if ( endPtr != NULL ) {
+ *endPtr = p;
+ }
+
+ /* Test for zero. */
+
+ if ( nSigDigs == 0 ) {
+ v = 0.0;
+ goto returnValue;
+ }
+
+ /*
+ * The easy cases are where we have an exact significand and
+ * the exponent is small enough that we can compute the value
+ * with only one roundoff. In addition to the cases where we
+ * can multiply or divide an exact-integer significand by an
+ * exact-integer power of 10, there is also David Gay's case
+ * where we can scale the significand by a power of 10 (still
+ * keeping it exact) and then multiply by an exact power of 10.
+ * The last case enables combinations like 83e25 that would
+ * otherwise require high precision arithmetic.
+ */
+
+ if ( nSigDigs <= DBL_DIG ) {
+ if ( exponent >= 0 ) {
+ if ( exponent <= mmaxpow ) {
+ v = exactSignificand * pow10[ exponent ];
+ goto returnValue;
+ } else {
+ int diff = DBL_DIG - nSigDigs;
+ if ( exponent - diff <= mmaxpow ) {
+ volatile double factor = exactSignificand * pow10[ diff ];
+ v = factor * pow10[ exponent - diff ];
+ goto returnValue;
+ }
+ }
+ } else {
+ if ( exponent >= -mmaxpow ) {
+ v = exactSignificand / pow10[ -exponent ];
+ goto returnValue;
+ }
+ }
+ }
+
+ /*
+ * We don't have one of the easy cases, so we can't compute the
+ * scanned number exactly, and have to do it in multiple precision.
+ * Begin by testing for obvious overflows and underflows.
+ */
+
+ if ( nSigDigs + exponent - 1 > maxDigits ) {
+ v = HUGE_VAL;
+ errno = ERANGE;
+ goto returnValue;
+ }
+ if ( nSigDigs + exponent - 1 < minDigits ) {
+ errno = ERANGE;
+ v = 0.;
+ goto returnValue;
+ }
+
+ /*
+ * Nothing exceeds the boundaries of the tables, at least.
+ * Compute an approximate value for the number, with
+ * no possibility of overflow because we manage the exponent
+ * separately.
+ */
+
+ if ( nSigDigs > DBL_DIG ) {
+ expt2 = exponent + nSigDigs - DBL_DIG;
+ } else {
+ expt2 = exponent;
+ }
+ v = frexp( exactSignificand, &machexp );
+ if ( expt2 > 0 ) {
+ v = frexp( v * pow10[ expt2 & 0xf ], &j );
+ machexp += j;
+ for ( i = 4; i < 9; ++i ) {
+ if ( expt2 & ( 1 << i ) ) {
+ v = frexp( v * pow_10_2_n[ i ], &j );
+ machexp += j;
+ }
+ }
+ } else {
+ v = frexp( v / pow10[ (-expt2) & 0xf ], &j );
+ machexp += j;
+ for ( i = 4; i < 9; ++i ) {
+ if ( (-expt2) & ( 1 << i ) ) {
+ v = frexp( v / pow_10_2_n[ i ], &j );
+ machexp += j;
+ }
+ }
+ }
+
+ /*
+ * A first approximation is that the result will be v * 2 ** machexp.
+ * v is greater than or equal to 0.5 and less than 1.
+ * If machexp > DBL_MAX_EXP * log2(FLT_RADIX), there is an overflow.
+ * Constrain the result to the smallest representible number to avoid
+ * premature underflow.
+ */
+
+ if ( machexp > DBL_MAX_EXP * log2FLT_RADIX ) {
+ v = HUGE_VAL;
+ errno = ERANGE;
+ goto returnValue;
+ }
+
+ v = SafeLdExp( v, machexp );
+ if ( v < tiny ) {
+ v = tiny;
+ }
+
+ /* We have a first approximation in v. Now we need to refine it. */
+
+ v = RefineResult( v, startOfSignificand, nSigDigs, exponent );
+
+ /* In a very few cases, a second iteration is needed. e.g., 457e-102 */
+
+ v = RefineResult( v, startOfSignificand, nSigDigs, exponent );
+
+ /* Handle underflow */
+
+ returnValue:
+ if ( nSigDigs != 0 && v == 0.0 ) {
+ errno = ERANGE;
+ }
+
+ /* Return a number with correct sign */
+
+ if ( signum ) {
+ v = -v;
+ }
+
+ /* Restore FPU mode word */
+
+#if defined(__GNUC__) && defined(__i386)
+ _FPU_SETCW( oldRoundingMode );
+#endif
+
+ return v;
+
+ /* Come here on an invalid input */
+
+ error:
+ if ( endPtr != NULL ) {
+ *endPtr = s;
+ }
+
+ /* Restore FPU mode word */
+
+#if defined(__GNUC__) && defined(__i386)
+ _FPU_SETCW( oldRoundingMode );
+#endif
+ return 0.0;
+
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * RefineResult --
+ *
+ * Given a poor approximation to a floating point number, returns
+ * a better one (The better approximation is correct to within
+ * 1 ulp, and is entirely correct if the poor approximation is
+ * correct to 1 ulp.)
+ *
+ * Results:
+ * Returns the improved result.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+RefineResult( double approxResult,
+ /* Approximate result of conversion */
+ CONST char* sigStart,
+ /* Pointer to start of significand in
+ * input string. */
+ int nSigDigs, /* Number of significant digits */
+ long exponent ) /* Power of ten to multiply by significand */
+{
+
+ int M2, M5; /* Powers of 2 and of 5 needed to put
+ * the decimal and binary numbers over
+ * a common denominator. */
+ double significand; /* Sigificand of the binary number */
+ int binExponent; /* Exponent of the binary number */
+
+ int msb; /* Most significant bit position of an
+ * intermediate result */
+ int nDigits; /* Number of mp_digit's in an intermediate
+ * result */
+ mp_int twoMv; /* Approx binary value expressed as an
+ * exact integer scaled by the multiplier 2M */
+ mp_int twoMd; /* Exact decimal value expressed as an
+ * exact integer scaled by the multiplier 2M */
+ int scale; /* Scale factor for M */
+ int multiplier; /* Power of two to scale M */
+ double num, den; /* Numerator and denominator of the
+ * correction term */
+ double quot; /* Correction term */
+ double minincr; /* Lower bound on the absolute value
+ * of the correction term. */
+ int i;
+ CONST char* p;
+
+ /*
+ * The first approximation is always low. If we find that
+ * it's HUGE_VAL, we're done.
+ */
+
+ if ( approxResult == HUGE_VAL ) {
+ return approxResult;
+ }
+
+ /*
+ * Find a common denominator for the decimal and binary fractions.
+ * The common denominator will be 2**M2 + 5**M5.
+ */
+
+ significand = frexp( approxResult, &binExponent );
+ i = mantBits - binExponent;
+ if ( i < 0 ) {
+ M2 = 0;
+ } else {
+ M2 = i;
+ }
+ if ( exponent > 0 ) {
+ M5 = 0;
+ } else {
+ M5 = -exponent;
+ if ( (M5-1) > M2 ) {
+ M2 = M5-1;
+ }
+ }
+
+ /*
+ * The floating point number is significand*2**binExponent.
+ * The 2**-1 bit of the significand (the most significant)
+ * corresponds to the 2**(binExponent+M2 + 1) bit of 2*M2*v.
+ * Allocate enough digits to hold that quantity, then
+ * convert the significand to a large integer, scaled
+ * appropriately. Then multiply by the appropriate power of 5.
+ */
+
+ msb = binExponent + M2; /* 1008 */
+ nDigits = msb / DIGIT_BIT + 1;
+ mp_init_size( &twoMv, nDigits );
+ i = ( msb % DIGIT_BIT + 1 );
+ twoMv.used = nDigits;
+ significand *= SafeLdExp( 1.0, i );
+ while ( -- nDigits >= 0 ) {
+ twoMv.dp[nDigits] = (mp_digit) significand;
+ significand -= (mp_digit) significand;
+ significand = SafeLdExp( significand, DIGIT_BIT );
+ }
+ for ( i = 0; i <= 8; ++i ) {
+ if ( M5 & ( 1 << i ) ) {
+ mp_mul( &twoMv, pow5+i, &twoMv );
+ }
+ }
+
+ /*
+ * Collect the decimal significand as a high precision integer.
+ * The least significant bit corresponds to bit M2+exponent+1
+ * so it will need to be shifted left by that many bits after
+ * being multiplied by 5**(M5+exponent).
+ */
+
+ mp_init( &twoMd ); mp_zero( &twoMd );
+ i = nSigDigs;
+ for ( p = sigStart ; ; ++p ) {
+ char c = *p;
+ if ( isdigit( c ) ) {
+ mp_mul_d( &twoMd, (unsigned) 10, &twoMd );
+ mp_add_d( &twoMd, (unsigned) (c - '0'), &twoMd );
+ --i;
+ if ( i == 0 ) break;
+ }
+ }
+ for ( i = 0; i <= 8; ++i ) {
+ if ( (M5+exponent) & ( 1 << i ) ) {
+ mp_mul( &twoMd, pow5+i, &twoMd );
+ }
+ }
+ mp_mul_2d( &twoMd, M2+exponent+1, &twoMd );
+ mp_sub( &twoMd, &twoMv, &twoMd );
+
+ /*
+ * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
+ * term. Because 2M may well overflow a double, we need to scale the
+ * denominator by a factor of 2**binExponent-mantBits
+ */
+
+ scale = binExponent - mantBits - 1;
+
+ mp_set( &twoMv, 1 );
+ for ( i = 0; i <= 8; ++i ) {
+ if ( M5 & ( 1 << i ) ) {
+ mp_mul( &twoMv, pow5+i, &twoMv );
+ }
+ }
+ multiplier = M2 + scale + 1;
+ if ( multiplier > 0 ) {
+ mp_mul_2d( &twoMv, multiplier, &twoMv );
+ } else if ( multiplier < 0 ) {
+ mp_div_2d( &twoMv, -multiplier, &twoMv, NULL );
+ }
+
+ /*
+ * If the result is less than unity, the error is less than 1/2 unit
+ * in the last place, so there's no correction to make.
+ */
+
+ if ( mp_cmp_mag( &twoMd, &twoMv ) == MP_LT ) {
+ return approxResult;
+ }
+
+ /*
+ * Convert the numerator and denominator of the corrector term
+ * accurately to floating point numbers.
+ */
+
+ num = TclBignumToDouble( &twoMd );
+ den = TclBignumToDouble( &twoMv );
+
+ quot = SafeLdExp( num/den, scale );
+ minincr = SafeLdExp( 1.0, binExponent - mantBits );
+
+ if ( quot < 0. && quot > -minincr ) {
+ quot = -minincr;
+ } else if ( quot > 0. && quot < minincr ) {
+ quot = minincr;
+ }
+
+ mp_clear( &twoMd );
+ mp_clear( &twoMv );
+
+
+ return approxResult + quot;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * ParseNaN --
+ *
+ * Parses a "not a number" from an input string, and returns the
+ * double precision NaN corresponding to it.
+ *
+ * Side effects:
+ * Advances endPtr to follow any (hex) in the input string.
+ *
+ * If the NaN is followed by a left paren, a string of spaes
+ * and hexadecimal digits, and a right paren, endPtr is advanced
+ * to follow it.
+ *
+ * The string of hexadecimal digits is OR'ed into the resulting
+ * NaN, and the signum is set as well. Note that a signalling NaN
+ * is never returned.
+ *
+ *----------------------------------------------------------------------
+ */
+
+double
+ParseNaN( int signum, /* Flag == 1 if minus sign has been
+ * seen in front of NaN */
+ CONST char** endPtr )
+ /* Pointer-to-pointer to char following "NaN"
+ * in the input string */
+{
+ CONST char* p = *endPtr;
+ char c;
+ union {
+ Tcl_WideUInt iv;
+ double dv;
+ } theNaN;
+
+ /* Scan off a hex number in parentheses. Embedded blanks are ok. */
+
+ theNaN.iv = 0;
+ if ( *p == '(' ) {
+ ++p;
+ for ( ; ; ) {
+ c = *p++;
+ if ( isspace(c) ) {
+ continue;
+ } else if ( c == ')' ) {
+ *endPtr = p;
+ break;
+ } else if ( isdigit(c) ) {
+ c -= '0';
+ } else if ( c >= 'A' && c <= 'F' ) {
+ c = c - 'A' + 10;
+ } else if ( c >= 'a' && c <= 'f' ) {
+ c = c - 'a' + 10;
+ } else {
+ theNaN.iv = ( ((Tcl_WideUInt) 0x7ff8) << 48 )
+ | ( ((Tcl_WideUInt) signum) << 63 );
+ return theNaN.dv;
+ }
+ theNaN.iv = (theNaN.iv << 4) | c;
+ }
+ }
+
+ /*
+ * Mask the hex number down to the least significant 51 bits.
+ */
+
+ theNaN.iv &= ( ((Tcl_WideUInt) 1) << 51 ) - 1;
+ if ( signum ) {
+ theNaN.iv |= ((Tcl_WideUInt) 0xfff8) << 48;
+ } else {
+ theNaN.iv |= ((Tcl_WideUInt) 0x7ff8) << 48;
+ }
+
+ *endPtr = p;
+ return theNaN.dv;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclDoubleDigits --
+ *
+ * Converts a double to a string of digits.
+ *
+ * Results:
+ * Returns the position of the character in the string
+ * after which the decimal point should appear. Since
+ * the string contains only significant digits, the
+ * position may be less than zero or greater than the
+ * length of the string.
+ *
+ * Side effects:
+ * Stores the digits in the given buffer and sets 'signum'
+ * according to the sign of the number.
+ *
+ *----------------------------------------------------------------------
+ */
+
+int
+TclDoubleDigits( char * strPtr, /* Buffer in which to store the result,
+ * must have at least 18 chars */
+ double v, /* Number to convert. Must be
+ * finite, and not NaN */
+ int *signum ) /* Output: 1 if the number is negative.
+ * Should handle -0 correctly on the
+ * IEEE architecture. */
+{
+
+ double f; /* Significand of v */
+
+ int e; /* Power of FLT_RADIX that satisfies
+ * v = f * FLT_RADIX**e */
+
+
+ int low_ok;
+ int high_ok;
+
+ mp_int r; /* Scaled significand */
+ mp_int s; /* Divisor such that v = r / s */
+ mp_int mplus; /* Scaled epsilon: (r + 2* mplus) ==
+ * v(+) where v(+) is the floating point
+ * successor of v. */
+ mp_int mminus; /* Scaled epsilon: (r - 2*mminus ) ==
+ * v(-) where v(-) is the floating point
+ * predecessor of v. */
+ mp_int temp;
+
+ int rfac2 = 0; /* Powers of 2 and 5 by which large */
+ int rfac5 = 0; /* integers should be scaled */
+ int sfac2 = 0;
+ int sfac5 = 0;
+ int mplusfac2 = 0;
+ int mminusfac2 = 0;
+
+ double a;
+ char c;
+ int i, k, n;
+
+ /*
+ * Take the absolute value of the number, and report the number's
+ * sign. Take special steps to preserve signed zeroes in IEEE floating
+ * point. (We can't use fpclassify, because that's a C9x feature and
+ * we still have to build on C89 compilers.)
+ */
+
+#ifndef IEEE_FLOATING_POINT
+ if ( v >= 0.0 ) {
+ *signum = 0;
+ } else {
+ *signum = 1;
+ v = -v;
+ }
+#else
+ union {
+ Tcl_WideUInt iv;
+ double dv;
+ } bitwhack;
+ bitwhack.dv = v;
+ if ( bitwhack.iv & ( (Tcl_WideUInt) 1 << 63 ) ) {
+ *signum = 1;
+ bitwhack.iv &= ~( (Tcl_WideUInt) 1 << 63 );
+ v = bitwhack.dv;
+ } else {
+ *signum = 0;
+ }
+#endif
+
+ /* Handle zero specially */
+
+ if ( v == 0.0 ) {
+ *strPtr++ = '0';
+ *strPtr++ = '\0';
+ return 1;
+ }
+
+ /*
+ * Develop f and e such that v = f * FLT_RADIX**e, with
+ * 1.0/FLT_RADIX <= f < 1.
+ */
+
+ f = frexp( v, &e );
+ n = e % log2FLT_RADIX;
+ if ( n > 0 ) {
+ n -= log2FLT_RADIX;
+ e += 1;
+ }
+ f *= ldexp( 1.0, n );
+ e = ( e - n ) / log2FLT_RADIX;
+ if ( f == 1.0 ) {
+ f = 1.0 / FLT_RADIX;
+ e += 1;
+ }
+
+ /*
+ * If the original number was denormalized, adjust e and f to be
+ * denormal as well.
+ */
+
+ if ( e < DBL_MIN_EXP ) {
+ n = mantBits + ( e - DBL_MIN_EXP ) * log2FLT_RADIX;
+ f = ldexp( f, ( e - DBL_MIN_EXP ) * log2FLT_RADIX );
+ e = DBL_MIN_EXP;
+ n = ( n + DIGIT_BIT - 1 ) / DIGIT_BIT;
+ } else {
+ n = mantDIGIT;
+ }
+
+ /*
+ * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
+ * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e
+ * by adjusting e.
+ */
+
+ a = f;
+ n = mantDIGIT;
+ mp_init_size( &r, n );
+ r.used = n;
+ r.sign = MP_ZPOS;
+ i = ( mantBits % DIGIT_BIT );
+ if ( i == 0 ) {
+ i = DIGIT_BIT;
+ }
+ while ( n > 0 ) {
+ a *= ldexp( 1.0, i );
+ i = DIGIT_BIT;
+ r.dp[--n] = (mp_digit) a;
+ a -= (mp_digit) a;
+ }
+ e -= DBL_MANT_DIG;
+
+ low_ok = high_ok = ( mp_iseven( &r ) );
+
+ /*
+ * We are going to want to develop integers r, s, mplus, and mminus
+ * such that v = r / s, v(+)-v / 2 = mplus / s; v-v(-) / 2 = mminus / s
+ * and then scale either s or r, mplus, mminus by an appropriate
+ * power of ten.
+ *
+ * We actually do this by keeping track of the powers of 2 and 5
+ * by which f is multiplied to yield v and by which 1 is multiplied
+ * to yield s, mplus, and mminus.
+ */
+
+ if ( e >= 0 ) {
+
+ int bits = e * log2FLT_RADIX;
+
+ if ( f != 1.0 / FLT_RADIX ) {
+
+ /* Normal case, m+ and m- are both FLT_RADIX**e */
+
+ rfac2 += bits + 1;
+ sfac2 = 1;
+ mplusfac2 = bits;
+ mminusfac2 = bits;
+
+ } else {
+
+ /*
+ * If f is equal to the smallest significand, then we need another
+ * factor of FLT_RADIX in s to cope with stepping to
+ * the next smaller exponent when going to e's predecessor.
+ */
+
+ rfac2 += bits + log2FLT_RADIX - 1;
+ sfac2 = 1 + log2FLT_RADIX;
+ mplusfac2 = bits + log2FLT_RADIX;
+ mminusfac2 = bits;
+
+ }
+
+ } else {
+
+ /* v has digits after the binary point */
+
+ if ( e <= DBL_MIN_EXP - DBL_MANT_DIG
+ || f != 1.0 / FLT_RADIX ) {
+
+ /*
+ * Either f isn't the smallest significand or e is
+ * the smallest exponent. mplus and mminus will both be 1.
+ */
+
+ rfac2 += 1;
+ sfac2 = 1 - e * log2FLT_RADIX;
+ mplusfac2 = 0;
+ mminusfac2 = 0;
+
+ } else {
+
+ /*
+ * f is the smallest significand, but e is not the smallest
+ * exponent. We need to scale by FLT_RADIX again to cope
+ * with the fact that v's predecessor has a smaller exponent.
+ */
+
+ rfac2 += 1 + log2FLT_RADIX;
+ sfac2 = 1 + log2FLT_RADIX * ( 1 - e );
+ mplusfac2 = FLT_RADIX;
+ mminusfac2 = 0;
+
+ }
+ }
+
+ /*
+ * Estimate the highest power of ten that will be
+ * needed to hold the result.
+ */
+
+ k = (int) ceil( log( v ) / log( 10. ) );
+ if ( k >= 0 ) {
+ sfac2 += k;
+ sfac5 = k;
+ } else {
+ rfac2 -= k;
+ mplusfac2 -= k;
+ mminusfac2 -= k;
+ rfac5 = -k;
+ }
+
+ /*
+ * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5.
+ */
+
+ mp_init_set( &mplus, 1 );
+ for ( i = 0; i <= 8; ++i ) {
+ if ( rfac5 & ( 1 << i ) ) {
+ mp_mul( &mplus, pow5+i, &mplus );
+ }
+ }
+ mp_mul( &r, &mplus, &r );
+ mp_mul_2d( &r, rfac2, &r );
+ mp_init_copy( &mminus, &mplus );
+ mp_mul_2d( &mplus, mplusfac2, &mplus );
+ mp_mul_2d( &mminus, mminusfac2, &mminus );
+ mp_init_set( &s, 1 );
+ for ( i = 0; i <= 8; ++i ) {
+ if ( sfac5 & ( 1 << i ) ) {
+ mp_mul( &s, pow5+i, &s );
+ }
+ }
+ mp_mul_2d( &s, sfac2, &s );
+
+ /*
+ * It is possible for k to be off by one because we used an
+ * inexact logarithm.
+ */
+
+ mp_init( &temp );
+ mp_add( &r, &mplus, &temp );
+ i = mp_cmp_mag( &temp, &s );
+ if ( i > 0 || ( high_ok && i == 0 ) ) {
+ mp_mul_d( &s, 10, &s );
+ ++k;
+ } else {
+ mp_mul_d( &temp, 10, &temp );
+ i = mp_cmp_mag( &temp, &s );
+ if ( i < 0 || ( high_ok && i == 0 ) ) {
+ mp_mul_d( &r, 10, &r );
+ mp_mul_d( &mplus, 10, &mplus );
+ mp_mul_d( &mminus, 10, &mminus );
+ --k;
+ }
+ }
+
+ /*
+ * At this point, k contains the power of ten by which we're
+ * scaling the result. r/s is at least 1/10 and strictly less
+ * than ten, and v = r/s * 10**k. mplus and mminus give the
+ * rounding limits.
+ */
+
+ for ( ; ; ) {
+ int tc1, tc2;
+ mp_mul_d( &r, 10, &r );
+ mp_div( &r, &s, &temp, &r ); /* temp = 10r / s; r = 10r mod s */
+ i = temp.dp[0];
+ mp_mul_d( &mplus, 10, &mplus );
+ mp_mul_d( &mminus, 10, &mminus );
+ tc1 = mp_cmp_mag( &r, &mminus );
+ if ( low_ok ) {
+ tc1 = ( tc1 <= 0 );
+ } else {
+ tc1 = ( tc1 < 0 );
+ }
+ mp_add( &r, &mplus, &temp );
+ tc2 = mp_cmp_mag( &temp, &s );
+ if ( high_ok ) {
+ tc2 = ( tc2 >= 0 );
+ } else {
+ tc2= ( tc2 > 0 );
+ }
+ if ( ! tc1 ) {
+ if ( !tc2 ) {
+ *strPtr++ = '0' + i;
+ } else {
+ c = (char) (i + '1');
+ break;
+ }
+ } else {
+ if ( !tc2 ) {
+ c = (char) (i + '0');
+ } else {
+ mp_mul_2d( &r, 1, &r );
+ n = mp_cmp_mag( &r, &s );
+ if ( n < 0 ) {
+ c = (char) (i + '0');
+ } else {
+ c = (char) (i + '1');
+ }
+ }
+ break;
+ }
+ };
+ *strPtr++ = c;
+ *strPtr++ = '\0';
+
+ /* Free memory */
+
+ mp_clear_multi( &r, &s, &mplus, &mminus, &temp, NULL );
+ return k;
+
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclInitDoubleConversion --
+ *
+ * Initializes constants that are needed for conversions to and
+ * from 'double'
+ *
+ * Results:
+ * None.
+ *
+ * Side effects:
+ * The log base 2 of the floating point radix, the number of
+ * bits in a double mantissa, and a table of the powers of five
+ * and ten are computed and stored.
+ *
+ *----------------------------------------------------------------------
+ */
+
+void
+TclInitDoubleConversion( void )
+{
+ int i;
+ int x;
+ double d;
+ if ( frexp( (double) FLT_RADIX, &log2FLT_RADIX ) != 0.5 ) {
+ Tcl_Panic( "This code doesn't work on a decimal machine!" );
+ }
+ --log2FLT_RADIX;
+ mantBits = DBL_MANT_DIG * log2FLT_RADIX;
+ d = 1.0;
+ x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log( 5.0 ));
+ if ( x < MAXPOW ) {
+ mmaxpow = x;
+ } else {
+ mmaxpow = MAXPOW;
+ }
+ for ( i = 0; i <= mmaxpow; ++i ) {
+ pow10[i] = d;
+ d *= 10.0;
+ }
+ for ( i = 0; i < 9; ++i ) {
+ mp_init( pow5 + i );
+ }
+ mp_set( pow5, 5 );
+ for ( i = 0; i < 8; ++i ) {
+ mp_sqr( pow5+i, pow5+i+1 );
+ }
+ tiny = SafeLdExp( 1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits );
+ maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX)
+ + 0.5 * log(10.))
+ / log( 10. ));
+ minDigits = (int) floor ( ( DBL_MIN_EXP - DBL_MANT_DIG )
+ * log( (double) FLT_RADIX ) / log( 10. ) );
+ mantDIGIT = ( mantBits + DIGIT_BIT - 1 ) / DIGIT_BIT;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclFinalizeDoubleConversion --
+ *
+ * Cleans up this file on exit.
+ *
+ * Results:
+ * None
+ *
+ * Side effects:
+ * Memory allocated by TclInitDoubleConversion is freed.
+ *
+ *----------------------------------------------------------------------
+ */
+
+void
+TclFinalizeDoubleConversion()
+{
+ int i;
+ for ( i = 0; i < 9; ++i ) {
+ mp_clear( pow5 + i );
+ }
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclBignumToDouble --
+ *
+ * Convert an arbitrary-precision integer to a native floating
+ * point number.
+ *
+ * Results:
+ * Returns the converted number. Sets errno to ERANGE if the
+ * number is too large to convert.
+ *
+ *----------------------------------------------------------------------
+ */
+
+double
+TclBignumToDouble( mp_int* a )
+ /* Integer to convert */
+{
+ mp_int b;
+ int bits;
+ int shift;
+ int i;
+ double r;
+
+ /* Determine how many bits we need, and extract that many from
+ * the input. Round to nearest unit in the last place. */
+
+ bits = mp_count_bits( a );
+ if ( bits > DBL_MAX_EXP * log2FLT_RADIX ) {
+ errno = ERANGE;
+ return HUGE_VAL;
+ }
+ shift = mantBits + 1 - bits;
+ mp_init( &b );
+ if ( shift > 0 ) {
+ mp_mul_2d( a, shift, &b );
+ } else if ( shift < 0 ) {
+ mp_div_2d( a, -shift, &b, NULL );
+ } else {
+ mp_copy( a, &b );
+ }
+ mp_add_d( &b, 1, &b );
+ mp_div_2d( &b, 1, &b, NULL );
+
+ /* Accumulate the result, one mp_digit at a time */
+
+ r = 0.0;
+ for ( i = b.used-1; i >= 0; --i ) {
+ r = ldexp( r, DIGIT_BIT ) + b.dp[i];
+ }
+ mp_clear( &b );
+
+ /* Scale the result to the correct number of bits. */
+
+ r = ldexp( r, bits - mantBits );
+
+ /* Return the result with the appropriate sign. */
+
+ if ( a->sign == MP_ZPOS ) {
+ return r;
+ } else {
+ return -r;
+ }
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * SafeLdExp --
+ *
+ * Do an 'ldexp' operation, but handle denormals gracefully.
+ *
+ * Results:
+ * Returns the appropriately scaled value.
+ *
+ * On some platforms, 'ldexp' fails when presented with a number
+ * too small to represent as a normalized double. This routine
+ * does 'ldexp' in two steps for those numbers, to return correctly
+ * denormalized values.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+SafeLdExp( double fract, int expt )
+{
+ int minexpt = DBL_MIN_EXP * log2FLT_RADIX;
+ volatile double a, b, retval;
+ if ( expt < minexpt ) {
+ a = ldexp( fract, expt - mantBits - minexpt );
+ b = ldexp( 1.0, mantBits + minexpt );
+ retval = a * b;
+ } else {
+ retval = ldexp( fract, expt );
+ }
+ return retval;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclFormatNaN --
+ *
+ * Makes the string representation of a "Not a Number"
+ *
+ * Results:
+ * None.
+ *
+ * Side effects:
+ * Stores the string representation in the supplied buffer,
+ * which must be at least TCL_DOUBLE_SPACE characters.
+ *
+ *----------------------------------------------------------------------
+ */
+
+void
+TclFormatNaN( double value, /* The Not-a-Number to format */
+ char* buffer ) /* String representation */
+{
+#ifndef IEEE_FLOATING_POINT
+ strcpy( buffer, "NaN" );
+ return;
+#else
+
+ union {
+ double dv;
+ Tcl_WideUInt iv;
+ } bitwhack;
+
+ bitwhack.dv = value;
+ if ( bitwhack.iv & ((Tcl_WideUInt) 1 << 63 ) ) {
+ bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63 );
+ *buffer++ = '-';
+ }
+ *buffer++ = 'N'; *buffer++ = 'a'; *buffer++ = 'N';
+ bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
+ if ( bitwhack.iv != 0 ) {
+ sprintf( buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv );
+ } else {
+ *buffer = '\0';
+ }
+
+#endif
+}