/* *---------------------------------------------------------------------- * * 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.5 2005/06/26 22:10:08 dkf Exp $ * *---------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include /* * 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 /* * HP's PA_RISC architecture uses 7ff4000000000000 to represent a * quiet NaN. Everyone else uses 7ff8000000000000. (Why, HP, why?) */ #ifdef __hppa # define NAN_START 0x7ff4 # define NAN_MASK (((Tcl_WideUInt) 1) << 50) #else # define NAN_START 0x7ff8 # define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif /* * There now follows a lot of static variables that are shared across all * threads but which are not guarded by mutexes. This is OK, because they are * only ever assigned _ONCE_ during Tcl's library initialization sequence. */ /* 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( UCHAR(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( UCHAR( 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( UCHAR( 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( UCHAR(c) ) ) { continue; } else if ( c == ')' ) { *endPtr = p; break; } else if ( isdigit( UCHAR(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) NAN_START) << 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) NAN_START) << 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 }