diff options
Diffstat (limited to 'generic/tclStrToD.c')
| -rwxr-xr-x | generic/tclStrToD.c | 2562 |
1 files changed, 1919 insertions, 643 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c index 68f6c30..7f8ecd5 100755 --- a/generic/tclStrToD.c +++ b/generic/tclStrToD.c @@ -1,21 +1,20 @@ /* *---------------------------------------------------------------------- * - * tclStrToD.c -- + * tclDouble.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 collection of procedures for managing + * conversions to/from floating-point in Tcl. They include + * TclParseNumber, which parses numbers from strings; TclDoubleDigits, + * which formats numbers into strings of digits, and procedures for + * interconversion among 'double' and 'mp_int' 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.8 2005/08/24 15:15:45 kennykb Exp $ + * RCS: @(#) $Id: tclStrToD.c,v 1.9 2005/10/08 14:42:45 dgp Exp $ * *---------------------------------------------------------------------- */ @@ -30,34 +29,42 @@ #include <tommath.h> /* - * The stuff below is a bit of a hack so that this file can be used in - * environments that include no UNIX, i.e. no errno: just arrange to use the - * errno from tclExecute.c here. + * Define TIP_114_FORMATS to accept 0b and 0o for binary and octal strings. + * Define KILL_OCTAL as well as TIP_114_FORMATS to suppress interpretation + * of numbers with leading zero as octal. (Ceterum censeo: numeros octonarios + * delendos esse.) */ -#ifdef TCL_GENERIC_ONLY -# define NO_ERRNO_H -#endif +#define TIP_114_FORMATS +#undef KILL_OCTAL -#ifdef NO_ERRNO_H -extern int errno; /* Use errno from tclExecute.c. */ -# define ERANGE 34 +#ifndef TIP_114_FORMATS +#undef KILL_OCTAL #endif +/* + * This code supports (at least hypothetically), IBM, Cray, VAX and + * IEEE-754 floating point; of these, only IEEE-754 can represent NaN. + * IEEE-754 can be uniquely determined by radix and by the widths of + * significand and exponent. + */ + #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. Most gcc-isms and ix86-isms are factored out here. + * gcc on x86 needs access to rounding controls, because of a questionable + * feature where it retains intermediate results as IEEE 'long double' values + * somewhat unpredictably. 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__ __volatile__ ("fnstcw %0" : "=m" (*&cw)) +#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw)) # define FPU_IEEE_ROUNDING 0x027f # define ADJUST_FPU_CONTROL_WORD #endif @@ -75,13 +82,26 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); # define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif -/* - * There now follows a lot of static variables that are shared across all - * threads but which are not guarded by mutexes. This is OK, because they are - * only ever assigned _ONCE_ during Tcl's library initialization sequence. - */ +/* The powers of ten that can be represented exactly as wide integers */ + +static int maxpow10_wide; +static Tcl_WideUInt *pow10_wide; + +/* The number of decimal digits that fit in an mp_digit */ + +static int log10_DIGIT_MAX; + +/* 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'. */ -static const double pow_10_2_n[] = { /* Inexact higher powers of ten */ +/* Inexact higher powers of ten */ + +static CONST double pow_10_2_n [] = { 1.0, 100.0, 10000.0, @@ -92,440 +112,1414 @@ static const double pow_10_2_n[] = { /* Inexact higher powers of ten */ 1.0e+128, 1.0e+256 }; -#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. */ + + /* Logarithm of the floating point radix. */ + +static int log2FLT_RADIX; + +/* Number of bits in a double's significand */ + +static int mantBits; + +/* Table of powers of 5**(2**n), up to 5**256 */ + +static mp_int pow5[9]; + +/* The smallest representable double */ + +static double tiny; + +/* The maximum number of digits to the left of the decimal point of a + * double. */ + +static int maxDigits; + +/* The maximum number of digits to the right of the decimal point in a + * double. */ + +static int minDigits; + +/* Number of mp_digit's needed to hold the significand of a double */ + +static int mantDIGIT; /* Static functions defined in this file */ -static double RefineResult(double approx, CONST char *start, int nDigits, - long exponent); -static double ParseNaN(int signum, CONST char **end); -static double SafeLdExp(double fraction, int exponent); +static int AccumulateDecimalDigit _ANSI_ARGS_((unsigned, int, + Tcl_WideUInt*, mp_int*, int)); +static double MakeLowPrecisionDouble _ANSI_ARGS_((int signum, + Tcl_WideUInt significand, + int nSigDigs, + int exponent)); +static double MakeHighPrecisionDouble _ANSI_ARGS_((int signum, + mp_int* significand, + int nSigDigs, + int exponent)); +static double MakeNaN _ANSI_ARGS_(( int signum, Tcl_WideUInt tag )); +static double RefineApproximation _ANSI_ARGS_((double approx, + mp_int* exactSignificand, + int exponent)); +static double AbsoluteValue(double v, int* signum); +static int GetIntegerTimesPower(double v, mp_int* r, int* e); +static double BignumToBiasedFrExp _ANSI_ARGS_(( mp_int* big, int* machexp )); +static double Pow10TimesFrExp _ANSI_ARGS_(( int exponent, + double fraction, + int* machexp )); +static double SafeLdExp _ANSI_ARGS_(( double fraction, int exponent )); + /* *---------------------------------------------------------------------- * - * TclStrToD -- + * TclParseNumber -- * - * Scans a double from a string. + * Place a "numeric" internal representation on a Tcl object. * * 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 a standard Tcl result. * * 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 an internal representation appropriate to the string. + * The internal representation may be an integer, a wide integer, + * a bignum, or a double. + * + * TclMakeObjNumeric is called as a common scanner in routines + * that expect numbers in Tcl_Obj's. It scans the string representation + * of a given Tcl_Obj and stores an internal rep that represents + * a "canonical" version of its numeric value. The value of the + * canonicalization is that a routine can determine simply by + * examining the type pointer whether an object LooksLikeInt, + * what size of integer is needed to hold it, and similar questions, + * and never needs to refer back to the string representation, even + * for "impure" objects. + * + * The 'strPtr' and 'endPtrPtr' arguments allow for recognizing a number + * that is in a substring of a Tcl_Obj, for example a screen metric or + * "end-" index. If 'strPtr' is not NULL, it designates where the + * number begins within the string. (The default is the start of + * objPtr's string rep, which will be constructed if necessary.) * - *------------------------------------------------------------------------ + * If 'strPtr' is supplied, 'objPtr' may be NULL. In this case, + * no internal representation will be generated; instead, the routine + * will simply check for a syntactically correct number, returning + * TCL_OK or TCL_ERROR as appropriate, and setting *endPtrPtr if + * necessary. + * + * If 'endPtrPtr' is not NULL, it designates the first character + * after the scanned number. In this case, successfully recognizing + * any digits will yield a return code of TCL_OK. Only in the case + * where no leading string of 'strPtr' (or of objPtr's internal rep) + * represents a number will TCL_ERROR be returned. + * + * When only a partial string is being recognized, it is the caller's + * responsibility to destroy the internal representation, or at + * least change its type. Failure to do so will lead to subsequent + * problems where a string that does not represent a number will + * be recognized as one because it has a numeric internal representation. + * + * When the 'flags' word includes TCL_PARSE_DECIMAL_ONLY, only decimal + * numbers are recognized; leading 0 has no special interpretation as + * octal and leading '0x' is forbidden. + * + *---------------------------------------------------------------------- */ -double -TclStrToD(CONST char *s, /* String to scan. */ - CONST char **endPtr) /* Pointer to the end of the scanned number. */ +int +TclParseNumber( Tcl_Interp* interp, + /* Tcl interpreter for error reporting. + * May be NULL */ + Tcl_Obj* objPtr, + /* Object to receive the internal rep */ + CONST char* type, + /* Type of number being parsed ("integer", + * "wide integer", etc. */ + CONST char* string, + /* Pointer to the start of the string to + * scan, see above */ + size_t length, /* Maximum length of the string to scan, + * see above. */ + CONST char** endPtrPtr, + /* (Output) pointer to the end of the + * scanned number, see above */ + int flags) /* Flags governing the parse */ { - 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. */ - 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; - /* - * With gcc on x86, the floating point rounding mode is double-extended. - * This causes the result of double-precision calculations to be rounded - * twice: once to the precision of double-extended and then again to the - * precision of double. Double-rounding introduces gratuitous errors of - * one ulp, so we need to change rounding mode to 53-bits. - */ - -#ifdef ADJUST_FPU_CONTROL_WORD - fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING; - fpu_control_t oldRoundingMode; - _FPU_GETCW(oldRoundingMode); - _FPU_SETCW(roundTo53Bits); -# define RestoreRoundingMode() _FPU_SETCW(oldRoundingMode) -#else -# define RestoreRoundingMode() (void) 0 /* Do nothing */ + enum State { + INITIAL, SIGNUM, ZERO, ZERO_X, +#ifdef TIP_114_FORMATS + ZERO_O, ZERO_B, BINARY, +#endif + HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL, + LEADING_RADIX_POINT, FRACTION, + EXPONENT_START, EXPONENT_SIGNUM, EXPONENT, + sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY +#ifdef IEEE_FLOATING_POINT + , sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH +#endif + } state = INITIAL; + enum State acceptState = INITIAL; + + int signum = 0; /* Sign of the number being parsed */ + Tcl_WideUInt significandWide = 0; + /* Significand of the number being + * parsed (if no overflow) */ + mp_int significandBig; /* Significand of the number being + * parsed (if it overflows significandWide) */ + int significandOverflow = 0; + /* Flag==1 iff significandBig is used */ + Tcl_WideUInt octalSignificandWide = 0; + /* Significand of an octal number; needed + * because we don't know whether a number + * with a leading zero is octal or decimal + * until we've scanned forward to a '.' or + * 'e' */ + mp_int octalSignificandBig; /* Significand of octal number once + * octalSignificandWide overflows */ + int octalSignificandOverflow = 0; + /* Flag==1 if octalSignificandBig is used */ + int numSigDigs = 0; /* Number of significant digits in the + * decimal significand */ + int numTrailZeros = 0; /* Number of trailing zeroes at the + * current point in the parse. */ + int numDigitsAfterDp = 0; /* Number of digits scanned after the + * decimal point */ + int exponentSignum = 0; /* Signum of the exponent of a floating + * point number */ + long exponent = 0; /* Exponent of a floating point number */ + CONST char* p; /* Pointer to next character to scan */ + size_t len; /* Number of characters remaining after p */ + CONST char* acceptPoint; /* Pointer to position after last character + * in an acceptable number */ + size_t acceptLen; /* Number of characters following that point */ + int status = TCL_OK; /* Status to return to caller */ + char d; /* Last hexadecimal digit scanned */ + int shift = 0; /* Amount to shift when accumulating binary */ +#ifdef TIP_114_FORMATS + int explicitOctal = 0; #endif - /* - * Discard leading whitespace from input. + /* + * Initialize string to start of the object's string rep if + * the caller didn't pass anything else. */ - while (isspace(UCHAR(*p))) { - ++p; + if ( string == NULL ) { + string = Tcl_GetStringFromObj( objPtr, NULL ); } - /* - * Determine the sign of the significand. - */ + p = string; + len = length; + acceptPoint = p; + acceptLen = len; + while ( 1 ) { + char c = len ? *p : '\0'; + switch (state) { - switch (*p) { - case '-': - signum = 1; + case INITIAL: + /* + * Initial state. Acceptable characters are +, -, digits, + * period, I, N, and whitespace. + */ + if (isspace(UCHAR(c))) { + break; + } else if (c == '+') { + state = SIGNUM; + break; + } else if (c == '-') { + signum = 1; + state = SIGNUM; + break; + } /* FALLTHROUGH */ - case '+': - ++p; - } - - /* - * Discard leading zeroes from input. - */ + + case SIGNUM: + /* + * Scanned a leading + or -. Acceptable characters are + * digits, period, I, and N. + */ + if (c == '0') { + if (flags & TCL_PARSE_DECIMAL_ONLY) { + state = DECIMAL; + } else { + state = ZERO; + } + break; + } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) { + goto zerox; + } else if (flags & TCL_PARSE_OCTAL_ONLY) { + goto zeroo; + } else if (isdigit(UCHAR(c))) { + significandWide = c - '0'; + numSigDigs = 1; + state = DECIMAL; + break; + } else if (flags & TCL_PARSE_INTEGER_ONLY) { + goto endgame; + } else if (c == '.') { + state = LEADING_RADIX_POINT; + break; + } else if (c == 'I' || c == 'i') { + state = sI; + break; +#ifdef IEEE_FLOATING_POINT + } else if (c == 'N' || c == 'n') { + state = sN; + break; +#endif + } + goto endgame; - while (*p == '0') { - seenDigit = 1; - ++p; - } + case ZERO: + /* + * Scanned a leading zero (perhaps with a + or -). + * Acceptable inputs are digits, period, X, and E. + * If 8 or 9 is encountered, the number can't be + * octal. This state and the OCTAL state differ only + * in whether they recognize 'X'. + */ + acceptState = state; + acceptPoint = p; + acceptLen = len; + if (c == 'x' || c == 'X') { + state = ZERO_X; + break; + } + if (flags & TCL_PARSE_HEXADECIMAL_ONLY) { + goto zerox; + } +#ifdef TIP_114_FORMATS + if (flags & TCL_PARSE_SCAN_PREFIXES) { + goto zeroo; + } + if (c == 'b' || c == 'B') { + state = ZERO_B; + break; + } + if (c == 'o' || c == 'O') { + explicitOctal = 1; + state = ZERO_O; + break; + } +#ifdef KILL_OCTAL + goto decimal; +#endif +#endif + /* FALLTHROUGH */ - /* - * 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. - */ + case OCTAL: + /* + * Scanned an optional + or -, followed by a string of + * octal digits. Acceptable inputs are more digits, + * period, or E. If 8 or 9 is encountered, commit to + * floating point. + */ + acceptState = state; + acceptPoint = p; + acceptLen = len; +#ifdef TIP_114_FORMATS + /* FALLTHROUGH */ + case ZERO_O: +#endif + zeroo: + if (c == '0') { + ++numTrailZeros; + state = OCTAL; + break; + } else if (c >= '1' && c <= '7') { + if (objPtr != NULL) { + shift = 3 * (numTrailZeros + 1); + significandOverflow = + AccumulateDecimalDigit((unsigned)(c-'0'), + numTrailZeros, + &significandWide, + &significandBig, + significandOverflow); + + if (!octalSignificandOverflow) { + /* + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check + * for too large shifts first. + */ + if ((octalSignificandWide != 0) + && ((shift >= CHAR_BIT*sizeof(Tcl_WideUInt)) + || (octalSignificandWide + > (~(Tcl_WideUInt)0 >> shift)))) { + octalSignificandOverflow = 1; + TclBNInitBignumFromWideUInt(&octalSignificandBig, + octalSignificandWide); + } + } + if (!octalSignificandOverflow) { + octalSignificandWide + = (octalSignificandWide << shift) + (c - '0'); + } else { + mp_mul_2d(&octalSignificandBig, shift, + &octalSignificandBig); + mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'), + &octalSignificandBig); + } + } + if ( numSigDigs != 0 ) { + numSigDigs += ( numTrailZeros + 1 ); + } else { + numSigDigs = 1; + } + numTrailZeros = 0; + state = OCTAL; + break; + } + /* FALLTHROUGH */ - for (;;) { - c = *p; - if (c == '.' && !seenDp) { - seenDp = 1; - ++p; - } else if (isdigit(UCHAR(c))) { + case BAD_OCTAL: +#ifdef TIP_114_FORMATS + if (explicitOctal) { + /* No forgiveness for bad digits in explicitly octal numbers */ + goto endgame; + } +#endif + if (flags & TCL_PARSE_INTEGER_ONLY) { + /* No seeking floating point when parsing only integer */ + goto endgame; + } +#ifndef KILL_OCTAL + /* + * Scanned a number with a leading zero that contains an + * 8, 9, radix point or E. This is an invalid octal number, + * but might still be floating point. + */ if (c == '0') { - if (startOfSignificand != NULL) { - ++nTrailZero; + ++numTrailZeros; + state = BAD_OCTAL; + break; + } else if (isdigit(UCHAR(c))) { + if (objPtr != NULL) { + significandOverflow = + AccumulateDecimalDigit((unsigned)(c-'0'), + numTrailZeros, + &significandWide, + &significandBig, + significandOverflow); + } + if ( numSigDigs != 0 ) { + numSigDigs += ( numTrailZeros + 1 ); + } else { + numSigDigs = 1; } + numTrailZeros = 0; + state = BAD_OCTAL; + break; + } else if (c == '.') { + state = FRACTION; + break; + } else if (c == 'E' || c == 'e') { + state = EXPONENT_START; + break; + } +#endif + goto endgame; + + /* + * Scanned 0x. If state is HEXADECIMAL, scanned at least + * one character following the 0x. The only acceptable + * inputs are hexadecimal digits. + */ + case HEXADECIMAL: + acceptState = state; + acceptPoint = p; + acceptLen = len; + /* FALLTHROUGH */ + case ZERO_X: + zerox: + if (c == '0') { + ++numTrailZeros; + state = HEXADECIMAL; + break; + } else if (isdigit(UCHAR(c))) { + d = (c-'0'); + } else if (c >= 'A' && c <= 'F') { + d = (c-'A'+10); + } else if (c >= 'a' && c <= 'f') { + d = (c-'a'+10); } else { - if (startOfSignificand == NULL) { - startOfSignificand = p; - } else if (nTrailZero) { - if (nTrailZero + nSigDigs < DBL_DIG) { - exactSignificand *= pow10[nTrailZero]; - } else if (nSigDigs < DBL_DIG) { - exactSignificand *= pow10[DBL_DIG - nSigDigs]; + goto endgame; + } + if (objPtr != NULL) { + shift = 4 * (numTrailZeros + 1); + if (!significandOverflow) { + /* + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check + * for too large shifts first. + */ + if (significandWide != 0 + && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) + || significandWide > (~(Tcl_WideUInt)0 >> shift))) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); } - nSigDigs += nTrailZero; } - if (nSigDigs < DBL_DIG) { - exactSignificand = 10. * exactSignificand + (c - '0'); + if (!significandOverflow) { + significandWide + = (significandWide << shift) + d; + } else { + mp_mul_2d(&significandBig, shift, + &significandBig); + mp_add_d(&significandBig, (mp_digit) d, + &significandBig); } - ++nSigDigs; - nTrailZero = 0; } - if (seenDp) { - ++nDigitsAfterDp; + numTrailZeros = 0; + state = HEXADECIMAL; + break; + +#ifdef TIP_114_FORMATS + case BINARY: + acceptState = state; + acceptPoint = p; + acceptLen = len; + case ZERO_B: + if (c == '0') { + ++numTrailZeros; + state = BINARY; + break; + } else if (c != '1') { + goto endgame; } - seenDigit = 1; - ++p; - } else { + if (objPtr != NULL) { + shift = numTrailZeros + 1; + if (!significandOverflow) { + /* + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check + * for too large shifts first. + */ + if (significandWide != 0 + && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) + || significandWide > (~(Tcl_WideUInt)0 >> shift))) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); + } + } + if (!significandOverflow) { + significandWide + = (significandWide << shift) + 1; + } else { + mp_mul_2d(&significandBig, shift, + &significandBig); + mp_add_d(&significandBig, (mp_digit) 1, + &significandBig); + } + } + numTrailZeros = 0; + state = BINARY; break; - } - } +#endif - /* - * 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 - */ + case DECIMAL: + /* + * Scanned an optional + or - followed by a string of + * decimal digits. + */ +#ifdef KILL_OCTAL + decimal: +#endif + acceptState = state; + acceptPoint = p; + acceptLen = len; + if (c == '0') { + ++numTrailZeros; + state = DECIMAL; + break; + } else if (isdigit(UCHAR(c))) { + if (objPtr != NULL) { + significandOverflow = + AccumulateDecimalDigit((unsigned)(c - '0'), + numTrailZeros, + &significandWide, + &significandBig, + significandOverflow); + } + numSigDigs += ( numTrailZeros + 1 ); + numTrailZeros = 0; + state = DECIMAL; + break; + } else if (flags & TCL_PARSE_INTEGER_ONLY) { + goto endgame; + } else if (c == '.') { + state = FRACTION; + break; + } else if (c == 'E' || c == 'e') { + state = EXPONENT_START; + break; + } + goto endgame; - exponent = 0; - if (seenDigit && (*p == 'e' || *p == 'E')) { - const char* stringSave = p; - ++p; - c = *p; - if (isdigit(UCHAR(c)) || c == '+' || c == '-') { - errno = 0; - exponent = strtol(p, (char**)&p, 10); - if (errno == ERANGE) { - if (exponent > 0) { - v = HUGE_VAL; + /* + * Found a decimal point. If no digits have yet been scanned, + * E is not allowed; otherwise, it introduces the exponent. + * If at least one digit has been found, we have a possible + * complete number. + */ + case FRACTION: + acceptState = state; + acceptPoint = p; + acceptLen = len; + if (c == 'E' || c=='e') { + state = EXPONENT_START; + break; + } + /* FALLTHROUGH */ + case LEADING_RADIX_POINT: + if (c == '0') { + ++numDigitsAfterDp; + ++numTrailZeros; + state = FRACTION; + break; + } else if (isdigit(UCHAR(c))) { + ++numDigitsAfterDp; + if (objPtr != NULL) { + significandOverflow = + AccumulateDecimalDigit((unsigned)(c-'0'), + numTrailZeros, + &significandWide, + &significandBig, + significandOverflow); + } + if ( numSigDigs != 0 ) { + numSigDigs += ( numTrailZeros + 1 ); } else { - v = 0.0; + numSigDigs = 1; } - *endPtr = p; - goto returnValue; + numTrailZeros = 0; + state = FRACTION; + break; } - } - if (p == stringSave+1) { - p = stringSave; - exponent = 0; - } - } - exponent += nTrailZero - nDigitsAfterDp; + goto endgame; - /* - * If we come here with no significant digits, we might still be looking - * at Inf or NaN. Go parse them. - */ + case EXPONENT_START: + /* + * Scanned the E at the start of an exponent. Make sure + * a legal character follows before using the C library + * strtol routine, which allows whitespace. + */ + if (c == '+') { + state = EXPONENT_SIGNUM; + break; + } else if (c == '-') { + exponentSignum = 1; + state = EXPONENT_SIGNUM; + break; + } + /* FALLTHROUGH */ - if (!seenDigit) { - /* - * Test for Inf or Infinity (in any case). - */ + case EXPONENT_SIGNUM: + /* + * Found the E at the start of the exponent, followed by + * a sign character. + */ + if (isdigit(UCHAR(c))) { + exponent = c - '0'; + state = EXPONENT; + break; + } + goto endgame; - if (c == 'I' || c == 'i') { - if ((p[1] == 'N' || p[1] == 'n') - && (p[2] == 'F' || p[2] == 'f')) { - p += 3; - if ((p[0] == 'I' || p[0] == 'i') - && (p[1] == 'N' || p[1] == 'n') - && (p[2] == 'I' || p[2] == 'i') - && (p[3] == 'T' || p[3] == 't') - && (p[4] == 'Y' || p[1] == 'y')) { - p += 5; - } - errno = ERANGE; - v = HUGE_VAL; - if (endPtr != NULL) { - *endPtr = p; + case EXPONENT: + /* + * Found an exponent with at least one digit. + * Accumulate it, making sure to hard-pin it to LONG_MAX + * on overflow. + */ + acceptState = state; + acceptPoint = p; + acceptLen = len; + if (isdigit(UCHAR(c))) { + if (exponent < (LONG_MAX - 9) / 10) { + exponent = 10 * exponent + (c - '0'); + } else { + exponent = LONG_MAX; } - goto returnValue; + state = EXPONENT; + break; } + goto endgame; -#ifdef IEEE_FLOATING_POINT /* - * Only IEEE floating point supports NaN + * Parse out INFINITY by simply spelling it out. + * INF is accepted as an abbreviation; other prefices are + * not. */ - } else if ((c == 'N' || c == 'n') - && (sizeof(Tcl_WideUInt) == sizeof(double))) { - if ((p[1] == 'A' || p[1] == 'a') - && (p[2] == 'N' || p[2] == 'n')) { - p += 3; - - if (endPtr != NULL) { - *endPtr = p; - } + case sI: + if ( c == 'n' || c == 'N' ) { + state = sIN; + break; + } + goto endgame; + case sIN: + if ( c == 'f' || c == 'F' ) { + state = sINF; + break; + } + goto endgame; + case sINF: + acceptState = state; + acceptPoint = p; + acceptLen = len; + if ( c == 'i' || c == 'I' ) { + state = sINFI; + break; + } + goto endgame; + case sINFI: + if ( c == 'n' || c == 'N' ) { + state = sINFIN; + break; + } + goto endgame; + case sINFIN: + if ( c == 'i' || c == 'I' ) { + state = sINFINI; + break; + } + goto endgame; + case sINFINI: + if ( c == 't' || c == 'T' ) { + state = sINFINIT; + break; + } + goto endgame; + case sINFINIT: + if ( c == 'y' || c == 'Y' ) { + state = sINFINITY; + break; + } + goto endgame; - /* - * Restore FPU mode word. - */ + /* + * Parse NaN's. + */ +#ifdef IEEE_FLOATING_POINT + case sN: + if ( c == 'a' || c == 'A' ) { + state = sNA; + break; + } + goto endgame; + case sNA: + if ( c == 'n' || c == 'N' ) { + state = sNAN; + break; + } + case sNAN: + acceptState = state; + acceptPoint = p; + acceptLen = len; + if ( c == '(' ) { + state = sNANPAREN; + break; + } + goto endgame; - RestoreRoundingMode(); - return ParseNaN(signum, endPtr); + /* + * Parse NaN(hexdigits) + */ + case sNANHEX: + if ( c == ')' ) { + state = sNANFINISH; + break; } + /* FALLTHROUGH */ + case sNANPAREN: + if ( isspace(UCHAR(c)) ) { + break; + } + if ( numSigDigs < 13 ) { + if ( c >= '0' && c <= '9' ) { + d = c - '0'; + } else if ( c >= 'a' && c <= 'f' ) { + d = 10 + c - 'a'; + } else if ( c >= 'A' && c <= 'F' ) { + d = 10 + c - 'A'; + } + significandWide = (significandWide << 4) + d; + state = sNANHEX; + break; + } + goto endgame; + case sNANFINISH: #endif + case sINFINITY: + acceptState = state; + acceptPoint = p; + acceptLen = len; + goto endgame; } + ++p; + --len; + } + + endgame: - goto error; + /* Back up to the last accepting state in the lexer */ + + if (acceptState == INITIAL) { + status = TCL_ERROR; } + p = acceptPoint; + len = acceptLen; - /* - * We've successfully scanned; update the end-of-element pointer. - */ + /* Skip past trailing whitespace */ - if (endPtr != NULL) { - *endPtr = p; + if (endPtrPtr != NULL) { + *endPtrPtr = p; } - /* - * Test for zero. - */ + while (len > 0 && isspace(UCHAR(*p))) { + ++p; + --len; + } - if (nSigDigs == 0) { - v = 0.0; - goto returnValue; + /* Determine whether a partial string is acceptable. */ + + if (endPtrPtr == NULL && len != 0 && *p != '\0') { + status = TCL_ERROR; } - /* - * 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. - */ + /* Generate and store the appropriate internal rep */ - if (nSigDigs <= DBL_DIG) { - if (exponent >= 0) { - if (exponent <= mmaxpow) { - v = exactSignificand * pow10[exponent]; - goto returnValue; - } else { - int diff = DBL_DIG - nSigDigs; - if (exponent - diff <= mmaxpow) { - volatile double factor = exactSignificand * pow10[diff]; - v = factor * pow10[exponent - diff]; - goto returnValue; + if (status == TCL_OK && objPtr != NULL) { + if ( acceptState != INITIAL ) { + TclFreeIntRep( objPtr ); + } + switch (acceptState) { + + case INITIAL: + status = TCL_ERROR; + break; + + case SIGNUM: + case BAD_OCTAL: + case ZERO_X: +#ifdef TIP_114_FORMATS + case ZERO_O: + case ZERO_B: +#endif + case LEADING_RADIX_POINT: + case EXPONENT_START: + case EXPONENT_SIGNUM: + case sI: + case sIN: + case sINFI: + case sINFIN: + case sINFINI: + case sINFINIT: + case sN: + case sNA: + case sNANPAREN: + case sNANHEX: + panic("in TclParseNumber: bad acceptState, can't happen."); + +#ifdef TIP_114_FORMATS + case BINARY: + shift = numTrailZeros; + if (!significandOverflow) { + if (significandWide !=0 + && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) + || significandWide + > (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); } } - } else if (exponent >= -mmaxpow) { - v = exactSignificand / pow10[-exponent]; - goto returnValue; + if (shift) { + if ( !significandOverflow ) { + significandWide <<= shift; + } else { + mp_mul_2d( &significandBig, shift, &significandBig ); + } + } + goto returnInteger; +#endif + case HEXADECIMAL: + /* Returning a hex integer. Final scaling step */ + shift = 4 * numTrailZeros; + if (!significandOverflow) { + if (significandWide !=0 + && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) + || significandWide + > (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); + } + } + if (shift) { + if ( !significandOverflow ) { + significandWide <<= shift; + } else { + mp_mul_2d( &significandBig, shift, &significandBig ); + } + } + goto returnInteger; + + case OCTAL: + /* Returning an octal integer. Final scaling step */ + shift = 3 * numTrailZeros; + if (!octalSignificandOverflow) { + if (octalSignificandWide != 0 + && (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) + || octalSignificandWide + > (((~(Tcl_WideUInt)0) >> 1) + signum) >> shift )) { + octalSignificandOverflow = 1; + TclBNInitBignumFromWideUInt(&octalSignificandBig, + octalSignificandWide); + } + } + if ( shift ) { + if ( !octalSignificandOverflow ) { + octalSignificandWide <<= shift; + } else { + mp_mul_2d( &octalSignificandBig, shift, + &octalSignificandBig ); + } + } + if (!octalSignificandOverflow) { + if (octalSignificandWide > + (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { +#ifndef NO_WIDE_TYPE + if (octalSignificandWide + <= (((~(Tcl_WideUInt)0) >> 1) + signum)) { + objPtr->typePtr = &tclWideIntType; + if (signum) { + objPtr->internalRep.wideValue = + - (Tcl_WideInt) octalSignificandWide; + } else { + objPtr->internalRep.wideValue = + (Tcl_WideInt) octalSignificandWide; + } + break; + } +#endif + TclBNInitBignumFromWideUInt(&octalSignificandBig, + octalSignificandWide); + octalSignificandOverflow = 1; + } else { + objPtr->typePtr = &tclIntType; + if (signum) { + objPtr->internalRep.longValue = + - (long) octalSignificandWide; + } else { + objPtr->internalRep.longValue = + (long) octalSignificandWide; + } + } + } + if (octalSignificandOverflow) { + if (signum) { + mp_neg(&octalSignificandBig, &octalSignificandBig); + } + TclSetBignumIntRep(objPtr, &octalSignificandBig); + } + break; + + case ZERO: + case DECIMAL: + significandOverflow = + AccumulateDecimalDigit( 0, numTrailZeros-1, + &significandWide, &significandBig, + significandOverflow ); + if (!significandOverflow + && (significandWide + > (((~(Tcl_WideUInt)0) >> 1) + signum))) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); + } + returnInteger: + if (!significandOverflow) { + if (significandWide > + (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { +#ifndef NO_WIDE_TYPE + if (significandWide + <= (((~(Tcl_WideUInt)0) >> 1) + signum)) { + objPtr->typePtr = &tclWideIntType; + if (signum) { + objPtr->internalRep.wideValue = + - (Tcl_WideInt) significandWide; + } else { + objPtr->internalRep.wideValue = + (Tcl_WideInt) significandWide; + } + break; + } +#endif + TclBNInitBignumFromWideUInt(&significandBig, + significandWide); + significandOverflow = 1; + } else { + objPtr->typePtr = &tclIntType; + if (signum) { + objPtr->internalRep.longValue = + - (long) significandWide; + } else { + objPtr->internalRep.longValue = + (long) significandWide; + } + } + } + if (significandOverflow) { + if (signum) { + mp_neg(&significandBig, &significandBig); + } + TclSetBignumIntRep(objPtr, &significandBig); + } + break; + + case FRACTION: + case EXPONENT: + + /* + * Here, we're parsing a floating-point number. + * 'significandWide' or 'significandBig' contains the + * exact significand, according to whether + * 'significandOverflow' is set. The desired floating + * point value is significand * 10**k, where + * k = numTrailZeros+exponent-numDigitsAfterDp. + */ + + objPtr->typePtr = &tclDoubleType; + if ( exponentSignum ) { + exponent = - exponent; + } + if ( !significandOverflow ) { + objPtr->internalRep.doubleValue = + MakeLowPrecisionDouble( signum, + significandWide, + numSigDigs, + ( numTrailZeros + + exponent + - numDigitsAfterDp ) ); + } else { + objPtr->internalRep.doubleValue = + MakeHighPrecisionDouble( signum, + &significandBig, + numSigDigs, + ( numTrailZeros + + exponent + - numDigitsAfterDp ) ); + } + break; + + case sINF: + case sINFINITY: + if ( signum ) { + objPtr->internalRep.doubleValue = -HUGE_VAL; + } else { + objPtr->internalRep.doubleValue = HUGE_VAL; + } + objPtr->typePtr = &tclDoubleType; + break; + + case sNAN: + case sNANFINISH: + objPtr->internalRep.doubleValue + = MakeNaN( signum, significandWide ); + objPtr->typePtr = &tclDoubleType; + break; + } } - /* - * 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. - */ + /* Format an error message when an invalid number is encountered. */ + + if ( status != TCL_OK ) { + if ( interp != NULL ) { + Tcl_Obj *msg = Tcl_NewStringObj( "expected ", -1 ); + Tcl_AppendToObj( msg, type, -1 ); + Tcl_AppendToObj( msg, " but got \"", -1 ); + TclAppendLimitedToObj( msg, string, length, 50, "" ); + Tcl_AppendToObj( msg, "\"", -1 ); + if ( state == BAD_OCTAL ) { + Tcl_AppendToObj( msg, " (looks like invalid octal number)", + -1 ); + } + Tcl_SetObjResult( interp, msg ); + } + } - if (nSigDigs + exponent - 1 > maxDigits) { - v = HUGE_VAL; - errno = ERANGE; - goto returnValue; + /* Free memory */ + + if (octalSignificandOverflow) { + mp_clear(&octalSignificandBig); } - if (nSigDigs + exponent - 1 < minDigits) { - errno = ERANGE; - v = 0.; - goto returnValue; + if (significandOverflow) { + mp_clear(&significandBig); } + return status; +} + +/* + *---------------------------------------------------------------------- + * + * AccumulateDecimalDigit -- + * + * Consume a decimal digit in a number being scanned. + * + * Results: + * Returns 1 if the number has overflowed to a bignum, 0 if it + * still fits in a wide integer. + * + * Side effects: + * Updates either the wide or bignum representation. + * + *---------------------------------------------------------------------- + */ - /* - * Nothing exceeds the boundaries of the tables, at least. Compute an - * approximate value for the number, with no possibility of overflow - * because we manage the exponent separately. - */ - - if (nSigDigs > DBL_DIG) { - expt2 = exponent + nSigDigs - DBL_DIG; - } else { - expt2 = exponent; - } - v = frexp(exactSignificand, &machexp); - if (expt2 > 0) { - v = frexp(v * pow10[expt2 & 0xf], &j); - machexp += j; - for (i=4 ; i<9 ; ++i) { - if (expt2 & (1 << i)) { - v = frexp(v * pow_10_2_n[i], &j); - machexp += j; +static int +AccumulateDecimalDigit( unsigned digit, + /* Digit being scanned */ + int numZeros, + /* Count of zero digits preceding the + * digit being scanned */ + Tcl_WideUInt* wideRepPtr, + /* Representation of the partial number + * as a wide integer */ + mp_int* bignumRepPtr, + /* Representation of the partial number + * as a bignum */ + int bignumFlag ) + /* Flag == 1 if the number overflowed + * previous to this digit. */ +{ + int i, n; + + /* Check if the number still fits in a wide */ + + if (!bignumFlag) { + if (*wideRepPtr != 0) { + if ((numZeros >= maxpow10_wide) + || (*wideRepPtr > (((~(Tcl_WideUInt)0) - digit) + / pow10_wide[numZeros+1]))) { + /* Oops, it's overflowed, have to allocate a bignum */ + TclBNInitBignumFromWideUInt (bignumRepPtr, *wideRepPtr); + bignumFlag = 1; } } + } + + /* Multiply the number by 10**numZeros+1 and add in the new digit. */ + + if (!bignumFlag) { + + /* Wide multiplication */ + + *wideRepPtr = *wideRepPtr * pow10_wide[numZeros+1] + digit; + } else if (numZeros < log10_DIGIT_MAX ) { + + /* Up to about 8 zeros - single digit multiplication */ + + mp_mul_d (bignumRepPtr, (mp_digit) pow10_wide[numZeros+1], + bignumRepPtr); + mp_add_d (bignumRepPtr, (mp_digit) digit, bignumRepPtr); + } else { - v = frexp(v / pow10[(-expt2) & 0xf], &j); - machexp += j; - for (i=4 ; i<9 ; ++i) { - if ((-expt2) & (1 << i)) { - v = frexp(v / pow_10_2_n[i], &j); - machexp += j; + + /* + * More than single digit multiplication. Multiply by the appropriate + * small powers of 5, and then shift. Large strings of zeroes are + * eaten 256 at a time; this is less efficient than it could be, + * but seems implausible. We presume that DIGIT_BIT is at least 27. + * The first multiplication, by up to 10**7, is done with a + * one-DIGIT multiply (this presumes that DIGIT_BIT >= 24). + */ + + n = numZeros + 1; + mp_mul_d (bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr); + for (i = 3; i <= 7; ++i) { + if (n & (1 << i)) { + mp_mul (bignumRepPtr, pow5+i, bignumRepPtr); } } + while (n >= 256) { + mp_mul (bignumRepPtr, pow5+8, bignumRepPtr); + n -= 256; + } + mp_mul_2d (bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr); } + return bignumFlag; +} + +/* + *---------------------------------------------------------------------- + * + * MakeLowPrecisionDouble -- + * + * Makes the double precision number, signum*significand*10**exponent. + * + * Results: + * Returns the constructed number. + * + * Common cases, where there are few enough digits that the number can + * be represented with at most roundoff, are handled specially here. + * If the number requires more than one rounded operation to compute, + * the code promotes the significand to a bignum and calls + * MakeHighPrecisionDouble to do it instead. + * + *---------------------------------------------------------------------- + */ + +static double +MakeLowPrecisionDouble( int signum, + /* 1 if the number is negative, 0 otherwise */ + Tcl_WideUInt significand, + /* Significand of the number */ + int numSigDigs, + /* Number of digits in the significand */ + int exponent ) + /* Power of ten */ +{ + double retval; /* Value of the number */ + mp_int significandBig; /* Significand expressed as a bignum */ + /* - * 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. + * With gcc on x86, the floating point rounding mode is double-extended. + * This causes the result of double-precision calculations to be rounded + * twice: once to the precision of double-extended and then again to the + * precision of double. Double-rounding introduces gratuitous errors of + * 1 ulp, so we need to change rounding mode to 53-bits. */ - if (machexp > DBL_MAX_EXP * log2FLT_RADIX) { - v = HUGE_VAL; - errno = ERANGE; - goto returnValue; - } +#if defined(__GNUC__) && defined(__i386) + fpu_control_t roundTo53Bits = 0x027f; + fpu_control_t oldRoundingMode; + _FPU_GETCW( oldRoundingMode ); + _FPU_SETCW( roundTo53Bits ); +#endif - v = SafeLdExp(v, machexp); - if (v < tiny) { - v = tiny; - } + /* Test for the easy cases */ - /* - * We have a first approximation in v. Now we need to refine it. - */ + if ( numSigDigs <= DBL_DIG ) { + if ( exponent >= 0 ) { + if ( exponent <= mmaxpow ) { - v = RefineResult(v, startOfSignificand, nSigDigs, exponent); + /* + * The significand is an exact integer, and so is + * 10**exponent. The product will be correct to within + * 1/2 ulp without special handling. + */ - /* - * In a very few cases, a second iteration is needed. e.g., 457e-102 - */ + retval = (double)(Tcl_WideInt)significand * pow10[ exponent ]; + goto returnValue; - v = RefineResult(v, startOfSignificand, nSigDigs, exponent); + } else { + int diff = DBL_DIG - numSigDigs; + if ( exponent-diff <= mmaxpow ) { + + /* + * 10**exponent is not an exact integer, but + * 10**(exponent-diff) is exact, and so is + * significand*10**diff, so we can still compute + * the value with only one roundoff. + */ + volatile double factor + = (double)(Tcl_WideInt)significand * pow10[diff]; + retval = factor * pow10[exponent-diff]; + goto returnValue; + } + } + } else { + if ( exponent >= -mmaxpow ) { + + /* + * 10**-exponent is an exact integer, and so is the + * significand. Compute the result by one division, + * again with only one rounding. + */ + + retval = (double)(Tcl_WideInt)significand / pow10[-exponent]; + goto returnValue; + } + } + } /* - * Handle underflow. + * All the easy cases have failed. Promote ths significand + * to bignum and call MakeHighPrecisionDouble to do it the hard way. */ + TclBNInitBignumFromWideUInt (&significandBig, significand); + retval = MakeHighPrecisionDouble( 0, &significandBig, numSigDigs, + exponent ); + + /* Come here to return the computed value */ + returnValue: - if (nSigDigs != 0 && v == 0.0) { - errno = ERANGE; + + if ( signum ) { + retval = -retval; } + /* On gcc on x86, restore the floating point mode word. */ + +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW( oldRoundingMode ); +#endif + + return retval; +} + +/* + *---------------------------------------------------------------------- + * + * MakeHighPrecisionDouble -- + * + * Makes the double precision number, signum*significand*10**exponent. + * + * Results: + * Returns the constructed number. + * + * MakeHighPrecisionDouble is used when arbitrary-precision arithmetic + * is needed to ensure correct rounding. It begins by calculating a + * low-precision approximation to the desired number, and then refines + * the answer in high precision. + * + *---------------------------------------------------------------------- + */ + +static double +MakeHighPrecisionDouble( int signum, + /* 1=negative, 0=nonnegative */ + mp_int* significand, + /* Exact significand of the number */ + int numSigDigs, + /* Number of significant digits */ + int exponent ) + /* Power of 10 by which to multiply */ +{ + + double retval; + int machexp; /* Machine exponent of a power of 10 */ + /* - * Return a number with correct sign. + * With gcc on x86, the floating point rounding mode is double-extended. + * This causes the result of double-precision calculations to be rounded + * twice: once to the precision of double-extended and then again to the + * precision of double. Double-rounding introduces gratuitous errors of + * 1 ulp, so we need to change rounding mode to 53-bits. */ - if (signum) { - v = -v; +#if defined(__GNUC__) && defined(__i386) + fpu_control_t roundTo53Bits = 0x027f; + fpu_control_t oldRoundingMode; + _FPU_GETCW( oldRoundingMode ); + _FPU_SETCW( roundTo53Bits ); +#endif + + /* Quick checks for over/underflow */ + + if ( numSigDigs + exponent - 1 > maxDigits ) { + retval = HUGE_VAL; + goto returnValue; + } + if ( numSigDigs + exponent - 1 < minDigits ) { + retval = 0; + goto returnValue; } - /* - * Restore FPU mode word and return. + /* + * Develop a first approximation to the significand. It is tempting + * simply to force bignum to double, but that will overflow on input + * numbers like 1.[string repeat 0 1000]1; while this is a not terribly + * likely scenario, we still have to deal with it. Use fraction and + * exponent instead. Once we have the significand, multiply by + * 10**exponent. Test for overflow. Convert back to a double, and + * test for underflow. */ - RestoreRoundingMode(); - return v; + retval = BignumToBiasedFrExp( significand, &machexp ); + retval = Pow10TimesFrExp( exponent, retval, &machexp ); + if ( machexp > DBL_MAX_EXP * log2FLT_RADIX ) { + retval = HUGE_VAL; + goto returnValue; + } + retval = SafeLdExp( retval, machexp ); + if ( retval < tiny ) { + retval = tiny; + } - /* - * Come here on an invalid input. + /* + * Refine the result twice. (The second refinement should be + * necessary only if the best approximation is a power of 2 + * minus 1/2 ulp). */ - error: - if (endPtr != NULL) { - *endPtr = s; + retval = RefineApproximation( retval, significand, exponent ); + retval = RefineApproximation( retval, significand, exponent ); + + /* Come here to return the computed value */ + + returnValue: + if ( signum ) { + retval = -retval; } - /* - * Restore FPU mode word and return. - */ + /* On gcc on x86, restore the floating point mode word. */ - RestoreRoundingMode(); - return 0.0; +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW( oldRoundingMode ); +#endif + return retval; } /* *---------------------------------------------------------------------- * - * RefineResult -- + * MakeNaN -- + * + * Makes a "Not a Number" given a set of bits to put in the + * tag bits * - * 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.) + * Note that a signalling NaN is never returned. + * + *---------------------------------------------------------------------- + */ + +#ifdef IEEE_FLOATING_POINT +static double +MakeNaN( int signum, /* Sign bit (1=negative, 0=nonnegative */ + Tcl_WideUInt tags ) /* Tag bits to put in the NaN */ +{ + union { + Tcl_WideUInt iv; + double dv; + } theNaN; + + theNaN.iv = tags; + theNaN.iv &= ( ((Tcl_WideUInt) 1) << 51 ) - 1; + if ( signum ) { + theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48; + } else { + theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48; + } + + return theNaN.dv; +} +#endif + +/* + *---------------------------------------------------------------------- + * + * RefineApproximation -- + * + * 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. @@ -534,119 +1528,109 @@ TclStrToD(CONST char *s, /* String to scan. */ */ 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. */ +RefineApproximation( double approxResult, + /* Approximate result of conversion */ + mp_int* exactSignificand, + /* Integer significand */ + int exponent ) + /* Power of 10 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; /* - * 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. + * Compute the large integer significand*2**(binExponent+M2+1) + * 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); - i = nSigDigs; - for (p=sigStart ;; ++p) { - char c = *p; - if (isdigit(UCHAR(c))) { - mp_mul_d(&twoMd, (unsigned) 10, &twoMd); - mp_add_d(&twoMd, (unsigned) (c - '0'), &twoMd); - --i; - if (i == 0) { - break; - } - } - } - for (i=0 ; i<=8 ; ++i) { - if ((M5+exponent) & (1 << i)) { - mp_mul(&twoMd, pow5+i, &twoMd); + mp_init_copy( &twoMd, exactSignificand ); + 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 @@ -656,135 +1640,55 @@ RefineResult(double approxResult, /* Approximate result of conversion. */ 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) { - mp_clear(&twoMd); - mp_clear(&twoMv); + 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; } /* *---------------------------------------------------------------------- * - * ParseNaN -- - * - * Parses a "not a number" from an input string, and returns the double - * precision NaN corresponding to it. - * - * Side effects: - * Advances endPtr to follow any (hex) in the input string. - * - * If the NaN is followed by a left paren, a string of spaes and - * hexadecimal digits, and a right paren, endPtr is advanced to follow - * it. - * - * The string of hexadecimal digits is OR'ed into the resulting NaN, and - * the signum is set as well. Note that a signalling NaN is never - * returned. - * - *---------------------------------------------------------------------- - */ - -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; - char c; - union { - Tcl_WideUInt iv; - double dv; - } theNaN; - - /* - * Scan off a hex number in parentheses. Embedded blanks are ok. - */ - - theNaN.iv = 0; - if (*p == '(') { - ++p; - for (;;) { - c = *p++; - if (isspace(UCHAR(c))) { - continue; - } else if (c == ')') { - *endPtr = p; - break; - } else if (isdigit(UCHAR(c))) { - c -= '0'; - } else if (c >= 'A' && c <= 'F') { - c -= 'A' + 10; - } else if (c >= 'a' && c <= 'f') { - c -= 'a' + 10; - } else { - theNaN.iv = (((Tcl_WideUInt) NAN_START) << 48) - | (((Tcl_WideUInt) signum) << 63); - return theNaN.dv; - } - theNaN.iv = (theNaN.iv << 4) | c; - } - } - - /* - * Mask the hex number down to the least significant 51 bits. - */ - - theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1; - if (signum) { - theNaN.iv |= ((Tcl_WideUInt) 0xfff8) << 48; - } else { - theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48; - } - - *endPtr = p; - return theNaN.dv; -} - -/* - *---------------------------------------------------------------------- - * * TclDoubleDigits -- * * Converts a double to a string of digits. @@ -803,20 +1707,21 @@ ParseNaN(int signum, /* Flag == 1 if minus sign has been seen in */ 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 * string, /* Buffer in which to store the result, + * must have at least 18 chars */ + double v, /* Number to convert. Must be + * finite, and not NaN */ + int *signum ) /* Output: 1 if the number is negative. + * Should handle -0 correctly on the + * IEEE architecture. */ { - double f; /* Significand of v. */ int e; /* Power of FLT_RADIX that satisfies * v = f * FLT_RADIX**e */ int lowOK, highOK; mp_int r; /* Scaled significand. */ mp_int s; /* Divisor such that v = r / s */ + int smallestSig; /* Flag == 1 iff v's significand is + * the smallest that can be represented. */ mp_int mplus; /* Scaled epsilon: (r + 2* mplus) == v(+) * where v(+) is the floating point successor * of v. */ @@ -830,103 +1735,31 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must 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.) - */ + /* Split the number into absolute value and signum. */ -#ifndef IEEE_FLOATING_POINT - if (v >= 0.0) { - *signum = 0; - } else { - *signum = 1; - v = -v; - } -#else - union { - Tcl_WideUInt iv; - double dv; - } bitwhack; - bitwhack.dv = v; - if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { - *signum = 1; - bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); - v = bitwhack.dv; - } else { - *signum = 0; - } -#endif + v = AbsoluteValue(v, signum); /* * Handle zero specially. */ - if (v == 0.0) { - *strPtr++ = '0'; - *strPtr++ = '\0'; + if ( v == 0.0 ) { + *string++ = '0'; + *string++ = '\0'; return 1; } - /* - * Develop f and e such that v = f * FLT_RADIX**e, with - * 1.0/FLT_RADIX <= f < 1. + /* + * Find a large integer r, and integer e, such that + * v = r * FLT_RADIX**e + * and r is as small as possible. Also determine whether the + * significand is the smallest possible. */ - f = frexp(v, &e); - n = e % log2FLT_RADIX; - if (n > 0) { - n -= log2FLT_RADIX; - e += 1; - } - f *= ldexp(1.0, n); - e = (e - n) / log2FLT_RADIX; - if (f == 1.0) { - f = 1.0 / FLT_RADIX; - e += 1; - } - - /* - * If the original number was denormalized, adjust e and f to be denormal - * as well. - */ - - if (e < DBL_MIN_EXP) { - n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX; - f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX); - e = DBL_MIN_EXP; - n = (n + DIGIT_BIT - 1) / DIGIT_BIT; - } else { - n = mantDIGIT; - } - - /* - * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision - * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by - * adjusting e. - */ - - a = f; - n = mantDIGIT; - mp_init_size(&r, n); - r.used = n; - r.sign = MP_ZPOS; - i = (mantBits % DIGIT_BIT); - if (i == 0) { - i = DIGIT_BIT; - } - while (n > 0) { - a *= ldexp(1.0, i); - i = DIGIT_BIT; - r.dp[--n] = (mp_digit) a; - a -= (mp_digit) a; - } - e -= DBL_MANT_DIG; + smallestSig = GetIntegerTimesPower(v, &r, &e); lowOK = highOK = (mp_iseven(&r)); @@ -943,12 +1776,12 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must if (e >= 0) { int bits = e * log2FLT_RADIX; - if (f != 1.0/FLT_RADIX) { + if (!smallestSig) { /* * Normal case, m+ and m- are both FLT_RADIX**e */ - rfac2 += bits + 1; + rfac2 = bits + 1; sfac2 = 1; mplusfac2 = bits; mminusfac2 = bits; @@ -959,7 +1792,7 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must * smaller exponent when going to e's predecessor. */ - rfac2 += bits + log2FLT_RADIX + 1; + rfac2 = bits + log2FLT_RADIX + 1; sfac2 = 1 + log2FLT_RADIX; mplusfac2 = bits + log2FLT_RADIX; mminusfac2 = bits; @@ -969,13 +1802,13 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must * v has digits after the binary point */ - if (e <= DBL_MIN_EXP-DBL_MANT_DIG || f != 1.0/FLT_RADIX) { + if (e <= DBL_MIN_EXP-DBL_MANT_DIG || !smallestSig) { /* * Either f isn't the smallest significand or e is the smallest * exponent. mplus and mminus will both be 1. */ - rfac2 += 1; + rfac2 = 1; sfac2 = 1 - e * log2FLT_RADIX; mplusfac2 = 0; mminusfac2 = 0; @@ -986,7 +1819,7 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must * fact that v's predecessor has a smaller exponent. */ - rfac2 += 1 + log2FLT_RADIX; + rfac2 = 1 + log2FLT_RADIX; sfac2 = 1 + log2FLT_RADIX * (1 - e); mplusfac2 = FLT_RADIX; mminusfac2 = 0; @@ -1081,9 +1914,9 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must } else { tc2= (tc2 > 0); } - if (!tc1) { - if (!tc2) { - *strPtr++ = '0' + i; + if ( ! tc1 ) { + if ( !tc2 ) { + *string++ = '0' + i; } else { c = (char) (i + '1'); break; @@ -1103,8 +1936,8 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must break; } }; - *strPtr++ = c; - *strPtr++ = '\0'; + *string++ = c; + *string++ = '\0'; /* * Free memory, and return. @@ -1117,6 +1950,148 @@ TclDoubleDigits(char * strPtr, /* Buffer in which to store the result, must /* *---------------------------------------------------------------------- * + * AbsoluteValue -- + * + * Splits a 'double' into its absolute value and sign. + * + * Results: + * Returns the absolute value. + * + * Side effects: + * Stores the signum in '*signum'. + * + *---------------------------------------------------------------------- + */ + +static double +AbsoluteValue (double v, /* Number to split */ + int* signum) /* (Output) Sign of the number 1=-, 0=+ */ +{ + /* + * Take the absolute value of the number, and report the number's sign. + * Take special steps to preserve signed zeroes in IEEE floating point. + * (We can't use fpclassify, because that's a C9x feature and we still + * have to build on C89 compilers.) + */ + +#ifndef IEEE_FLOATING_POINT + if (v >= 0.0) { + *signum = 0; + } else { + *signum = 1; + v = -v; + } +#else + union { + Tcl_WideUInt iv; + double dv; + } bitwhack; + bitwhack.dv = v; + if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { + *signum = 1; + bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); + v = bitwhack.dv; + } else { + *signum = 0; + } +#endif + return v; +} + +/* + *---------------------------------------------------------------------- + * + * GetIntegerTimesPower -- + * + * Converts a floating point number to an exact integer times a + * power of the floating point radix. + * + * Results: + * Returns 1 if it converted the smallest significand, 0 otherwise. + * + * Side effects: + * Initializes the integer value (does not just assign it), + * and stores the exponent. + * + *---------------------------------------------------------------------- + */ + +static int +GetIntegerTimesPower(double v, /* Value to convert */ + mp_int* rPtr, + /* (Output) Integer value */ + int* ePtr) /* (Output) Power of FLT_RADIX by which + * r must be multiplied to yield v*/ +{ + + double a; + double f; + int e; + int i; + int n; + + /* + * Develop f and e such that v = f * FLT_RADIX**e, with + * 1.0/FLT_RADIX <= f < 1. + */ + + f = frexp(v, &e); +#if FLT_RADIX > 2 + n = e % log2FLT_RADIX; + if (n > 0) { + n -= log2FLT_RADIX; + e += 1; + f *= ldexp(1.0, n); + } + e = (e - n) / log2FLT_RADIX; +#endif + if (f == 1.0) { + f = 1.0 / FLT_RADIX; + e += 1; + } + + /* + * If the original number was denormalized, adjust e and f to be denormal + * as well. + */ + + if (e < DBL_MIN_EXP) { + n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX; + f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX); + e = DBL_MIN_EXP; + n = (n + DIGIT_BIT - 1) / DIGIT_BIT; + } else { + n = mantDIGIT; + } + + /* + * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision + * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by + * adjusting e. + */ + + a = f; + n = mantDIGIT; + mp_init_size(rPtr, n); + rPtr->used = n; + rPtr->sign = MP_ZPOS; + i = (mantBits % DIGIT_BIT); + if (i == 0) { + i = DIGIT_BIT; + } + while (n > 0) { + a *= ldexp(1.0, i); + i = DIGIT_BIT; + rPtr->dp[--n] = (mp_digit) a; + a -= (mp_digit) a; + } + *ePtr = e - DBL_MANT_DIG; + return (f == 1.0 / FLT_RADIX); +} + +/* + *---------------------------------------------------------------------- + * * TclInitDoubleConversion -- * * Initializes constants that are needed for conversions to and from @@ -1138,16 +2113,43 @@ TclInitDoubleConversion(void) { int i; int x; + Tcl_WideUInt u; double d; - if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) { - Tcl_Panic("This code doesn't work on a decimal machine!"); + /* + * Initialize table of powers of 10 expressed as wide integers. + */ + + maxpow10_wide = + (int) floor(sizeof (Tcl_WideUInt) * CHAR_BIT * log (2.) / log (10.)); + pow10_wide = (Tcl_WideUInt*) Tcl_Alloc ((maxpow10_wide + 1) + * sizeof (Tcl_WideUInt)); + u = 1; + for (i = 0; i < maxpow10_wide; ++i) { + pow10_wide[i] = u; + u *= 10; + } + pow10_wide[i] = u; + + /* + * Determine how many bits of precision a double has, and how many + * decimal digits that represents. + */ + + 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) { + + /* + * Initialize a table of powers of ten that can be exactly represented + * in a double. + */ + + x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log( 5.0 )); + if ( x < MAXPOW ) { mmaxpow = x; } else { mmaxpow = MAXPOW; @@ -1156,19 +2158,32 @@ TclInitDoubleConversion(void) pow10[i] = d; d *= 10.0; } - for (i=0 ; i<9 ; ++i) { - mp_init(pow5 + i); + + /* Initialize a table of large powers of five. */ + + 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) + 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; + + /* + * Determine the number of decimal digits to the left and right of the + * decimal point in the largest and smallest double, the smallest double + * that differs from zero, and the number of mp_digits needed to represent + * the significand of a double. + */ + + 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; + log10_DIGIT_MAX = (int) floor (DIGIT_BIT * log(2.) / log (10.)); } /* @@ -1191,9 +2206,62 @@ void TclFinalizeDoubleConversion() { int i; - for (i=0 ; i<9 ; ++i) { - mp_clear(pow5 + i); + Tcl_Free ((char*)pow10_wide); + for ( i = 0; i < 9; ++i ) { + mp_clear( pow5 + i ); + } +} + +/* + *---------------------------------------------------------------------- + * + * TclInitBignumFromDouble -- + * + * Extracts the integer part of a double and converts it to + * an arbitrary precision integer. + * + * Results: + * None. + * + * Side effects: + * Initializes the bignum supplied, and stores the converted number + * in it. + * + *---------------------------------------------------------------------- + */ + +int +TclInitBignumFromDouble(Tcl_Interp *interp, /* For error message */ + double d, /* Number to convert */ + mp_int* b) /* Place to store the result */ +{ + double fract; + int expt; + + /* Infinite values can't convert to bignum */ + if (TclIsInfinite(d)) { + if (interp != NULL) { + char *s = "integer value too large to represent"; + Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1)); + Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL); + } + return TCL_ERROR; + } + fract = frexp(d,&expt); + if (expt <= 0) { + mp_init(b); + mp_zero(b); + } else { + Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits); + int shift = expt - mantBits; + TclBNInitBignumFromWideInt(b, w); + if (shift < 0) { + mp_div_2d(b, -shift, b, NULL); + } else if (shift > 0) { + mp_mul_2d(b, shift, b); + } } + return TCL_OK; } /* @@ -1228,7 +2296,11 @@ TclBignumToDouble(mp_int *a) /* Integer to convert. */ bits = mp_count_bits(a); if (bits > DBL_MAX_EXP*log2FLT_RADIX) { errno = ERANGE; - return HUGE_VAL; + if (a->sign == MP_ZPOS) { + return HUGE_VAL; + } else { + return -HUGE_VAL; + } } shift = mantBits + 1 - bits; mp_init(&b); @@ -1268,6 +2340,210 @@ TclBignumToDouble(mp_int *a) /* Integer to convert. */ return -r; } } + +double +TclCeil(mp_int *a) /* Integer to convert. */ +{ + double r = 0.0; + mp_int b; + + mp_init(&b); + if (mp_cmp_d(a, 0) == MP_LT) { + mp_neg(a, &b); + r = -TclFloor(&b); + } else { + int bits = mp_count_bits(a); + + if (bits > DBL_MAX_EXP*log2FLT_RADIX) { + r = HUGE_VAL; + } else { + int i, exact = 1, shift = mantBits - bits; + + if (shift > 0) { + mp_mul_2d(a, shift, &b); + } else if (shift < 0) { + mp_int d; + mp_init(&d); + mp_div_2d(a, -shift, &b, &d); + exact = mp_iszero(&d); + mp_clear(&d); + } else { + mp_copy(a, &b); + } + if (!exact) { + mp_add_d(&b, 1, &b); + } + for (i=b.used-1 ; i>=0 ; --i) { + r = ldexp(r, DIGIT_BIT) + b.dp[i]; + } + r = ldexp(r, bits - mantBits); + } + } + mp_clear(&b); + return r; +} + +double +TclFloor(mp_int *a) /* Integer to convert. */ +{ + double r = 0.0; + mp_int b; + + mp_init(&b); + if (mp_cmp_d(a, 0) == MP_LT) { + mp_neg(a, &b); + r = -TclCeil(&b); + } else { + int bits = mp_count_bits(a); + + if (bits > DBL_MAX_EXP*log2FLT_RADIX) { + r = DBL_MAX; + } else { + int i, shift = mantBits - bits; + + 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); + } + for (i=b.used-1 ; i>=0 ; --i) { + r = ldexp(r, DIGIT_BIT) + b.dp[i]; + } + r = ldexp(r, bits - mantBits); + } + } + mp_clear(&b); + return r; +} + +/* + *---------------------------------------------------------------------- + * + * BignumToBiasedFrExp -- + * + * Convert an arbitrary-precision integer to a native floating + * point number in the range [0.5,1) times a power of two. + * NOTE: Intentionally converts to a number that's a few + * ulp too small, so that RefineApproximation will not overflow + * near the high end of the machine's arithmetic range. + * + * Results: + * Returns the converted number. + * + * Side effects: + * Stores the exponent of two in 'machexp'. + * + *---------------------------------------------------------------------- + */ + +static double +BignumToBiasedFrExp( mp_int* a, + /* Integer to convert */ + int* machexp ) + /* Power of two */ +{ + mp_int b; + int bits; + int shift; + int i; + double r; + + /* Determine how many bits we need, and extract that many from + * the input. Round to nearest unit in the last place. */ + + bits = mp_count_bits( a ); + shift = mantBits - 2 - bits; + mp_init( &b ); + if ( shift > 0 ) { + mp_mul_2d( a, shift, &b ); + } else if ( shift < 0 ) { + mp_div_2d( a, -shift, &b, NULL ); + } else { + mp_copy( a, &b ); + } + + /* Accumulate the result, one mp_digit at a time */ + + r = 0.0; + for ( i = b.used-1; i >= 0; --i ) { + r = ldexp( r, DIGIT_BIT ) + b.dp[i]; + } + mp_clear( &b ); + + /* Return the result with the appropriate sign. */ + + *machexp = bits - mantBits + 2; + if ( a->sign == MP_ZPOS ) { + return r; + } else { + return -r; + } +} + +/* + *---------------------------------------------------------------------- + * + * Pow10TimesFrExp -- + * + * Multiply a power of ten by a number expressed as fraction and + * exponent. + * + * Results: + * Returns the significand of the result. + * + * Side effects: + * Overwrites the 'machexp' parameter with the exponent of the + * result. + * + * Assumes that 'exponent' is such that 10**exponent would be a double, + * even though 'fraction*10**(machexp+exponent)' might overflow. + * + *---------------------------------------------------------------------- + */ + +static double +Pow10TimesFrExp( int exponent, /* Power of 10 to multiply by */ + double fraction, + /* Significand of multiplicand */ + int* machexp ) /* On input, exponent of multiplicand. + * On output, exponent of result. */ +{ + int i, j; + int expt = *machexp; + double retval = fraction; + + if ( exponent > 0 ) { + + /* Multiply by 10**exponent */ + + retval = frexp( retval * pow10[ exponent & 0xf ], &j ); + expt += j; + for ( i = 4; i < 9; ++i ) { + if ( exponent & (1<<i) ) { + retval = frexp( retval * pow_10_2_n[ i ], &j ); + expt += j; + } + } + } else if ( exponent < 0 ) { + + /* Divide by 10**-exponent */ + + retval = frexp( retval / pow10[ (-exponent) & 0xf ], &j ); + expt += j; + for ( i = 4; i < 9; ++i ) { + if ( (-exponent) & (1<<i) ) { + retval = frexp( retval / pow_10_2_n[ i ], &j ); + expt += j; + } + } + } + + *machexp = expt; + return retval; + +} /* *---------------------------------------------------------------------- |