summaryrefslogtreecommitdiffstats
path: root/generic
diff options
context:
space:
mode:
Diffstat (limited to 'generic')
-rwxr-xr-xgeneric/tclStrToD.c1212
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 */
}
generated by cgit v0.12 at 2024-11-18 14:33:38 (GMT)