diff options
Diffstat (limited to 'generic/tclStrToD.c')
| -rwxr-xr-x | generic/tclStrToD.c | 1212 |
1 files changed, 590 insertions, 622 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c index 9f31841..ab2701d 100755 --- a/generic/tclStrToD.c +++ b/generic/tclStrToD.c @@ -3,19 +3,19 @@ * * 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. + * 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 $ + * RCS: @(#) $Id: tclStrToD.c,v 1.6 2005/07/05 11:37:02 dkf Exp $ * *---------------------------------------------------------------------- */ @@ -31,46 +31,48 @@ /* * 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. + * 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 +# define NO_ERRNO_H #endif #ifdef NO_ERRNO_H extern int errno; /* Use errno from tclExecute.c. */ -#define ERANGE 34 +# define ERANGE 34 #endif -#if ( FLT_RADIX == 2 ) && ( DBL_MANT_DIG == 53 ) && ( DBL_MAX_EXP == 1024 ) -#define IEEE_FLOATING_POINT +#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. + * 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. Most gcc-isms and ix86-isms are factored out here. */ #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)) +# define _FPU_GETCW(cw) __asm__ ("fnstcw %0" : "=m" (*&cw)) +# define _FPU_SETCW(cw) __asm__ ("fldcw %0" : : "m" (*&cw)) +# define FPU_IEEE_ROUNDING 0x027f +# define ADJUST_FPU_CONTROL_WORD #endif /* - * HP's PA_RISC architecture uses 7ff4000000000000 to represent a - * quiet NaN. Everyone else uses 7ff8000000000000. (Why, HP, why?) + * 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) +# define NAN_START 0x7ff4 +# define NAN_MASK (((Tcl_WideUInt) 1) << 50) #else -# define NAN_START 0x7ff8 -# define NAN_MASK (((Tcl_WideUInt) 1) << 51) +# define NAN_START 0x7ff8 +# define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif /* @@ -79,17 +81,7 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); * 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 [] = { +static const double pow_10_2_n[] = { /* Inexact higher powers of ten */ 1.0, 100.0, 10000.0, @@ -100,43 +92,29 @@ static CONST double pow_10_2_n [] = { 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; +#define MAXPOW 22 /* Num of exactly representable powers of 10 */ +static double pow10[MAXPOW+1]; /* The powers of ten that can be represented + * exactly as IEEE754 doubles. */ +static int mmaxpow; /* Largest power of ten that can be + * represented exactly in a 'double'. */ +static int log2FLT_RADIX; /* Logarithm of the floating point radix. */ +static int mantBits; /* Number of bits in a double's significand. */ +static mp_int pow5[9]; /* Table of powers of 5**(2**n), up to + * 5**256. */ +static double tiny; /* The smallest representable double. */ +static int maxDigits; /* The maximum number of digits to the left of + * the decimal point of a double. */ +static int minDigits; /* The maximum number of digits to the right + * of the decimal point in a double. */ +static int mantDIGIT; /* Number of mp_digit's needed to hold the + * significand of a double. */ /* 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 )); +static double RefineResult(double approx, CONST char *start, int nDigits, + long exponent); +static double ParseNaN(int signum, CONST char **end); +static double SafeLdExp(double fraction, int exponent); /* *---------------------------------------------------------------------- @@ -146,55 +124,46 @@ static double SafeLdExp _ANSI_ARGS_(( double fraction, int exponent )); * 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. + * 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. + * 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 */ +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 */ + 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 */ + /* 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. */ + volatile double v; /* Scanned value; must be 'volatile double' on + * gc-ix86 to force correct rounding to IEEE + * double and not Intel double-extended. */ + 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; /* @@ -202,25 +171,32 @@ TclStrToD( CONST char* s, * 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. + * one ulp, so we need to change rounding mode to 53-bits. */ -#if defined(__GNUC__) && defined(__i386) - fpu_control_t roundTo53Bits = 0x027f; +#ifdef ADJUST_FPU_CONTROL_WORD + fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING; fpu_control_t oldRoundingMode; - _FPU_GETCW( oldRoundingMode ); - _FPU_SETCW( roundTo53Bits ); + _FPU_GETCW(oldRoundingMode); + _FPU_SETCW(roundTo53Bits); +# define RestoreRoundingMode() _FPU_SETCW(oldRoundingMode) +#else +# define RestoreRoundingMode() (void) 0 /* Do nothing */ #endif - /* Discard leading whitespace */ + /* + * Discard leading whitespace from input. + */ - while ( isspace( *p ) ) { + while (isspace(UCHAR(*p))) { ++p; } - /* Determine the sign of the significand */ + /* + * Determine the sign of the significand. + */ - switch( *p ) { + switch (*p) { case '-': signum = 1; /* FALLTHROUGH */ @@ -228,48 +204,50 @@ TclStrToD( CONST char* s, ++p; } - /* Discard leading zeroes */ + /* + * Discard leading zeroes from input. + */ - while ( *p == '0' ) { + 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. + /* + * 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 ( ; ; ) { + for (;;) { c = *p; - if ( c == '.' && !seenDp ) { + if (c == '.' && !seenDp) { seenDp = 1; ++p; - } else if ( isdigit( UCHAR(c) ) ) { - if ( c == '0' ) { - if ( startOfSignificand != NULL ) { + } else if (isdigit(UCHAR(c))) { + if (c == '0') { + if (startOfSignificand != NULL) { ++nTrailZero; } } else { - if ( startOfSignificand == NULL ) { + 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 ]; + } 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 ) { + if (nSigDigs < DBL_DIG) { exactSignificand = 10. * exactSignificand + (c - '0'); } ++nSigDigs; nTrailZero = 0; } - if ( seenDp ) { + if (seenDp) { ++nDigitsAfterDp; } seenDigit = 1; @@ -280,29 +258,28 @@ TclStrToD( CONST char* s, } /* - * 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. + * 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; + if (seenDigit && (*p == 'e' || *p == 'E')) { + const char* stringSave = p; ++p; c = *p; - if ( isdigit( UCHAR( c ) ) || c == '+' || c == '-' ) { + if (isdigit(UCHAR(c)) || c == '+' || c == '-') { errno = 0; - exponent = strtol( p, (char**)&p, 10 ); - if ( errno == ERANGE ) { - if ( exponent > 0 ) { + exponent = strtol(p, (char**)&p, 10); + if (errno == ERANGE) { + if (exponent > 0) { v = HUGE_VAL; } else { v = 0.0; @@ -311,68 +288,64 @@ TclStrToD( CONST char* s, goto returnValue; } } - if ( p == stringSave + 1 ) { + if (p == stringSave+1) { p = stringSave; exponent = 0; } } - exponent = exponent + nTrailZero - nDigitsAfterDp; + exponent += nTrailZero - nDigitsAfterDp; /* - * If we come here with no significant digits, we might still be - * looking at Inf or NaN. Go parse them. + * 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 (!seenDigit) { + /* + * Test for Inf or Infinity (in any case). + */ - if ( ( p[1] == 'N' || p[1] == 'n' ) - && ( p[2] == 'F' || p[2] == 'f' ) ) { + 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' ) ) { + 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 ) { + 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' ) ) { + /* + * Only 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 ) { + + if (endPtr != NULL) { *endPtr = p; } - - /* Restore FPU mode word */ -#if defined(__GNUC__) && defined(__i386) - _FPU_SETCW( oldRoundingMode ); -#endif - return ParseNaN( signum, endPtr ); + /* + * Restore FPU mode word. + */ + RestoreRoundingMode(); + return ParseNaN(signum, endPtr); } #endif - } goto error; @@ -382,162 +355,167 @@ TclStrToD( CONST char* s, * We've successfully scanned; update the end-of-element pointer. */ - if ( endPtr != NULL ) { + if (endPtr != NULL) { *endPtr = p; } - /* Test for zero. */ + /* + * Test for zero. + */ - if ( nSigDigs == 0 ) { + 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. + * 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 ]; + 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 ]; + 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; - } + } 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. + /* + * 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 ) { + if (nSigDigs + exponent - 1 > maxDigits) { v = HUGE_VAL; errno = ERANGE; goto returnValue; } - if ( nSigDigs + exponent - 1 < minDigits ) { + 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. + * 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 ) { + 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 ); + 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 ); + 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 ); + 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 ); + 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. + * 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 ) { + if (machexp > DBL_MAX_EXP * log2FLT_RADIX) { v = HUGE_VAL; errno = ERANGE; goto returnValue; } - v = SafeLdExp( v, machexp ); - if ( v < tiny ) { + v = SafeLdExp(v, machexp); + if (v < tiny) { v = tiny; } - /* We have a first approximation in v. Now we need to refine it. */ + /* + * We have a first approximation in v. Now we need to refine it. + */ + + v = RefineResult(v, startOfSignificand, nSigDigs, exponent); - v = RefineResult( v, startOfSignificand, nSigDigs, exponent ); + /* + * In a very few cases, a second iteration is needed. e.g., 457e-102 + */ - /* In a very few cases, a second iteration is needed. e.g., 457e-102 */ - - v = RefineResult( v, startOfSignificand, nSigDigs, exponent ); + v = RefineResult(v, startOfSignificand, nSigDigs, exponent); - /* Handle underflow */ + /* + * Handle underflow. + */ returnValue: - if ( nSigDigs != 0 && v == 0.0 ) { + if (nSigDigs != 0 && v == 0.0) { errno = ERANGE; } - /* Return a number with correct sign */ + /* + * Return a number with correct sign. + */ - if ( signum ) { + if (signum) { v = -v; } - /* Restore FPU mode word */ - -#if defined(__GNUC__) && defined(__i386) - _FPU_SETCW( oldRoundingMode ); -#endif + /* + * Restore FPU mode word and return. + */ - return v; - - /* Come here on an invalid input */ + RestoreRoundingMode(); + return v; + + /* + * Come here on an invalid input. + */ error: - if ( endPtr != NULL ) { + if (endPtr != NULL) { *endPtr = s; } - /* Restore FPU mode word */ - -#if defined(__GNUC__) && defined(__i386) - _FPU_SETCW( oldRoundingMode ); -#endif - return 0.0; + /* + * Restore FPU mode word and return. + */ + RestoreRoundingMode(); + return 0.0; } /* @@ -545,10 +523,9 @@ TclStrToD( CONST char* s, * * 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.) + * 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. @@ -557,120 +534,119 @@ TclStrToD( CONST char* s, */ 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 */ +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 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 */ + * 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. */ + * 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; + const char* p; /* - * The first approximation is always low. If we find that - * it's HUGE_VAL, we're done. + * The first approximation is always low. If we find that it's HUGE_VAL, + * we're done. */ - if ( approxResult == HUGE_VAL ) { + 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. + * Find a common denominator for the decimal and binary fractions. The + * common denominator will be 2**M2 + 5**M5. */ - significand = frexp( approxResult, &binExponent ); + significand = frexp(approxResult, &binExponent); i = mantBits - binExponent; - if ( i < 0 ) { + if (i < 0) { M2 = 0; } else { M2 = i; } - if ( exponent > 0 ) { + if (exponent > 0) { M5 = 0; } else { M5 = -exponent; - if ( (M5-1) > M2 ) { + 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 + /* + * 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 */ + msb = binExponent + M2; /* 1008 */ nDigits = msb / DIGIT_BIT + 1; - mp_init_size( &twoMv, nDigits ); - i = ( 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 ) { + significand *= SafeLdExp(1.0, i); + while (--nDigits >= 0) { twoMv.dp[nDigits] = (mp_digit) significand; significand -= (mp_digit) significand; - significand = SafeLdExp( significand, DIGIT_BIT ); + significand = SafeLdExp(significand, DIGIT_BIT); } - for ( i = 0; i <= 8; ++i ) { - if ( M5 & ( 1 << i ) ) { - mp_mul( &twoMv, pow5+i, &twoMv ); + 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). + + /* + * 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 ); + mp_init(&twoMd); + mp_zero(&twoMd); i = nSigDigs; - for ( p = sigStart ; ; ++p ) { + 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 ); + if (isdigit(UCHAR(c))) { + mp_mul_d(&twoMd, (unsigned) 10, &twoMd); + mp_add_d(&twoMd, (unsigned) (c - '0'), &twoMd); --i; - if ( i == 0 ) break; + if (i == 0) { + break; + } } } - for ( i = 0; i <= 8; ++i ) { - if ( (M5+exponent) & ( 1 << i ) ) { - mp_mul( &twoMd, pow5+i, &twoMd ); + 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 ); + 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 @@ -680,49 +656,48 @@ RefineResult( double approxResult, scale = binExponent - mantBits - 1; - mp_set( &twoMv, 1 ); - for ( i = 0; i <= 8; ++i ) { - if ( M5 & ( 1 << i ) ) { - mp_mul( &twoMv, pow5+i, &twoMv ); + 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 (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 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 ) { + if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) { return approxResult; } - /* - * Convert the numerator and denominator of the corrector term - * accurately to floating point numbers. + /* + * Convert the numerator and denominator of the corrector term accurately + * to floating point numbers. */ - num = TclBignumToDouble( &twoMd ); - den = TclBignumToDouble( &twoMv ); + num = TclBignumToDouble(&twoMd); + den = TclBignumToDouble(&twoMv); - quot = SafeLdExp( num/den, scale ); - minincr = SafeLdExp( 1.0, binExponent - mantBits ); + quot = SafeLdExp(num/den, scale); + minincr = SafeLdExp(1.0, binExponent - mantBits); - if ( quot < 0. && quot > -minincr ) { + if (quot<0. && quot>-minincr) { quot = -minincr; - } else if ( quot > 0. && quot < minincr ) { + } else if (quot>0. && quot<minincr) { quot = minincr; } - mp_clear( &twoMd ); - mp_clear( &twoMv ); + mp_clear(&twoMd); + mp_clear(&twoMv); - return approxResult + quot; } @@ -731,70 +706,71 @@ RefineResult( double approxResult, * * ParseNaN -- * - * Parses a "not a number" from an input string, and returns the - * double precision NaN corresponding to it. + * 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. + * 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. + * 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 */ +static 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; + const char* p = *endPtr; char c; union { Tcl_WideUInt iv; double dv; } theNaN; - /* Scan off a hex number in parentheses. Embedded blanks are ok. */ + /* + * Scan off a hex number in parentheses. Embedded blanks are ok. + */ theNaN.iv = 0; - if ( *p == '(' ) { + if (*p == '(') { ++p; - for ( ; ; ) { + for (;;) { c = *p++; - if ( isspace( UCHAR(c) ) ) { + if (isspace(UCHAR(c))) { continue; - } else if ( c == ')' ) { + } else if (c == ')') { *endPtr = p; break; - } else if ( isdigit( UCHAR(c) ) ) { + } 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 if (c >= 'A' && c <= 'F') { + c -= 'A' + 10; + } else if (c >= 'a' && c <= 'f') { + c -= 'a' + 10; } else { - theNaN.iv = ( ((Tcl_WideUInt) NAN_START) << 48 ) - | ( ((Tcl_WideUInt) signum) << 63 ); + 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) 1) << 51) - 1; + if (signum) { theNaN.iv |= ((Tcl_WideUInt) 0xfff8) << 48; } else { theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48; @@ -812,68 +788,59 @@ ParseNaN( int signum, /* Flag == 1 if minus sign has been * 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. + * 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. + * 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. */ +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 */ - + 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 */ + int lowOK, highOK; + 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 + 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 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.) + /* + * 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 ) { + if (v >= 0.0) { *signum = 0; } else { *signum = 1; @@ -885,159 +852,152 @@ TclDoubleDigits( char * strPtr, /* Buffer in which to store the result, double dv; } bitwhack; bitwhack.dv = v; - if ( bitwhack.iv & ( (Tcl_WideUInt) 1 << 63 ) ) { + if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { *signum = 1; - bitwhack.iv &= ~( (Tcl_WideUInt) 1 << 63 ); + bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); v = bitwhack.dv; } else { *signum = 0; } #endif - /* Handle zero specially */ + /* + * Handle zero specially. + */ - if ( v == 0.0 ) { + 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 ); + f = frexp(v, &e); n = e % log2FLT_RADIX; - if ( n > 0 ) { + if (n > 0) { n -= log2FLT_RADIX; e += 1; } - f *= ldexp( 1.0, n ); - e = ( e - n ) / log2FLT_RADIX; - if ( f == 1.0 ) { + 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 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 ); + 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; + 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. + * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by + * adjusting e. */ - + a = f; n = mantDIGIT; - mp_init_size( &r, n ); + mp_init_size(&r, n); r.used = n; r.sign = MP_ZPOS; - i = ( mantBits % DIGIT_BIT ); - if ( i == 0 ) { + i = (mantBits % DIGIT_BIT); + if (i == 0) { i = DIGIT_BIT; } - while ( n > 0 ) { - a *= ldexp( 1.0, i ); + 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 ) ); + lowOK = highOK = (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 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. + * 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 ) { - + if (e >= 0) { int bits = e * log2FLT_RADIX; - if ( f != 1.0 / FLT_RADIX ) { - - /* Normal case, m+ and m- are both FLT_RADIX**e */ + 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. + * 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. + /* + * 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. + * 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 ); + 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. + /* + * Estimate the highest power of ten that will be needed to hold the + * result. */ - k = (int) ceil( log( v ) / log( 10. ) ); - if ( k >= 0 ) { + k = (int) ceil(log(v) / log(10.)); + if (k >= 0) { sfac2 += k; sfac5 = k; } else { @@ -1051,88 +1011,88 @@ TclDoubleDigits( char * strPtr, /* Buffer in which to store the result, * 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_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(&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 ); + mp_mul_2d(&s, sfac2, &s); /* - * It is possible for k to be off by one because we used an - * inexact logarithm. + * 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 ); + mp_init(&temp); + mp_add(&r, &mplus, &temp); + i = mp_cmp_mag(&temp, &s); + if (i>0 || (highOK && 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 ); + mp_mul_d(&temp, 10, &temp); + i = mp_cmp_mag(&temp, &s); + if (i<0 || (highOK && 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. + * 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 ( ; ; ) { + for (;;) { int tc1, tc2; - mp_mul_d( &r, 10, &r ); - mp_div( &r, &s, &temp, &r ); /* temp = 10r / s; r = 10r mod s */ + + 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 ); + mp_mul_d(&mplus, 10, &mplus); + mp_mul_d(&mminus, 10, &mminus); + tc1 = mp_cmp_mag(&r, &mminus); + if (lowOK) { + tc1 = (tc1 <= 0); } else { - tc1 = ( tc1 < 0 ); + tc1 = (tc1 < 0); } - mp_add( &r, &mplus, &temp ); - tc2 = mp_cmp_mag( &temp, &s ); - if ( high_ok ) { - tc2 = ( tc2 >= 0 ); + mp_add(&r, &mplus, &temp); + tc2 = mp_cmp_mag(&temp, &s); + if (highOK) { + tc2 = (tc2 >= 0); } else { - tc2= ( tc2 > 0 ); + tc2= (tc2 > 0); } - if ( ! tc1 ) { - if ( !tc2 ) { + if (!tc1) { + if (!tc2) { *strPtr++ = '0' + i; } else { c = (char) (i + '1'); break; } } else { - if ( !tc2 ) { + if (!tc2) { c = (char) (i + '0'); } else { - mp_mul_2d( &r, 1, &r ); - n = mp_cmp_mag( &r, &s ); - if ( n < 0 ) { + mp_mul_2d(&r, 1, &r); + n = mp_cmp_mag(&r, &s); + if (n < 0) { c = (char) (i + '0'); } else { c = (char) (i + '1'); @@ -1144,11 +1104,12 @@ TclDoubleDigits( char * strPtr, /* Buffer in which to store the result, *strPtr++ = c; *strPtr++ = '\0'; - /* Free memory */ + /* + * Free memory, and return. + */ - mp_clear_multi( &r, &s, &mplus, &mminus, &temp, NULL ); + mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL); return k; - } /* @@ -1156,56 +1117,56 @@ TclDoubleDigits( char * strPtr, /* Buffer in which to store the result, * * TclInitDoubleConversion -- * - * Initializes constants that are needed for conversions to and - * from 'double' + * 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. + * 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 ) +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!" ); + + 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 ) { + 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 ) { + for (i=0 ; i<=mmaxpow ; ++i) { pow10[i] = d; d *= 10.0; } - for ( i = 0; i < 9; ++i ) { - mp_init( pow5 + i ); + 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 ); + 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; + tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits); + maxDigits = (int) + ((DBL_MAX_EXP * log((double) FLT_RADIX) + log(10.)/2) / log(10.)); + minDigits = (int) + floor((DBL_MIN_EXP-DBL_MANT_DIG)*log((double)FLT_RADIX)/log(10.)); + mantDIGIT = (mantBits + DIGIT_BIT - 1) / DIGIT_BIT; } /* @@ -1228,8 +1189,8 @@ void TclFinalizeDoubleConversion() { int i; - for ( i = 0; i < 9; ++i ) { - mp_clear( pow5 + i ); + for (i=0 ; i<9 ; ++i) { + mp_clear(pow5 + i); } } @@ -1238,19 +1199,18 @@ TclFinalizeDoubleConversion() * * TclBignumToDouble -- * - * Convert an arbitrary-precision integer to a native floating - * point number. + * 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. + * Returns the converted number. Sets errno to ERANGE if the number is + * too large to convert. * *---------------------------------------------------------------------- */ double -TclBignumToDouble( mp_int* a ) - /* Integer to convert */ +TclBignumToDouble(mp_int *a) /* Integer to convert. */ { mp_int b; int bits; @@ -1258,46 +1218,54 @@ TclBignumToDouble( mp_int* a ) 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. */ + /* + * 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 ) { + 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 ); + 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_copy(a, &b); } - mp_add_d( &b, 1, &b ); - mp_div_2d( &b, 1, &b, NULL ); + mp_add_d(&b, 1, &b); + mp_div_2d(&b, 1, &b, NULL); - /* Accumulate the result, one mp_digit at a time */ + /* + * 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]; + for (i=b.used-1 ; i>=0 ; --i) { + r = ldexp(r, DIGIT_BIT) + b.dp[i]; } - mp_clear( &b ); + mp_clear(&b); - /* Scale the result to the correct number of bits. */ + /* + * Scale the result to the correct number of bits. + */ - r = ldexp( r, bits - mantBits ); + r = ldexp(r, bits - mantBits); - /* Return the result with the appropriate sign. */ + /* + * Return the result with the appropriate sign. + */ - if ( a->sign == MP_ZPOS ) { + if (a->sign == MP_ZPOS) { return r; } else { return -r; } -} +} /* *---------------------------------------------------------------------- @@ -1309,25 +1277,25 @@ TclBignumToDouble( mp_int* a ) * 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. + * 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 ) +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 ); + if (expt < minexpt) { + a = ldexp(fract, expt - mantBits - minexpt); + b = ldexp(1.0, mantBits + minexpt); retval = a * b; } else { - retval = ldexp( fract, expt ); + retval = ldexp(fract, expt); } return retval; } @@ -1343,38 +1311,38 @@ SafeLdExp( double fract, int expt ) * None. * * Side effects: - * Stores the string representation in the supplied buffer, - * which must be at least TCL_DOUBLE_SPACE characters. + * 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 */ +TclFormatNaN(double value, /* The Not-a-Number to format. */ + char *buffer) /* String representation. */ { #ifndef IEEE_FLOATING_POINT - strcpy( buffer, "NaN" ); + 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 ); + if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { + bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63); *buffer++ = '-'; } - *buffer++ = 'N'; *buffer++ = 'a'; *buffer++ = 'N'; + *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 ); + if (bitwhack.iv != 0) { + sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv); } else { *buffer = '\0'; } - -#endif +#endif /* IEEE_FLOATING_POINT */ } |