diff options
-rwxr-xr-x | generic/tclStrToD.c | 1447 |
1 files changed, 713 insertions, 734 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c index e4d7e18..7c2d6a8 100755 --- a/generic/tclStrToD.c +++ b/generic/tclStrToD.c @@ -1,20 +1,20 @@ /* *---------------------------------------------------------------------- * - * tclDouble.c -- + * tclStrToD.c -- * - * 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. + * 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. + * 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.11 2005/10/13 15:15:28 dgp Exp $ + * RCS: @(#) $Id: tclStrToD.c,v 1.12 2005/10/13 15:23:22 dkf Exp $ * *---------------------------------------------------------------------- */ @@ -30,8 +30,8 @@ /* * 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 + * 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.) */ @@ -43,10 +43,9 @@ #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. + * 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) @@ -56,9 +55,9 @@ /* * 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. + * 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) @@ -71,7 +70,7 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); /* * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN. - * Everyone else uses 7ff8000000000000. (Why, HP, why?) + * Everyone else uses 7ff8000000000000. (Why, HP, why?) */ #ifdef __hppa @@ -82,26 +81,33 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); # define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif -/* The powers of ten that can be represented exactly as wide integers */ +/* + * Constants used by this file (most of which are only ever calculated at + * runtime). + */ -static int maxpow10_wide; +static int maxpow10_wide; /* The powers of ten that can be represented + * exactly as wide integers. */ 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]; - +#define MAXPOW 22 +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'. */ - -/* Inexact higher powers of ten */ - -static CONST double pow_10_2_n [] = { +static int log10_DIGIT_MAX; /* The number of decimal digits that fit in an + * mp_digit. */ +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 CONST double pow_10_2_n[] = { /* Inexact higher powers of ten. */ 1.0, 100.0, 10000.0, @@ -113,60 +119,26 @@ static CONST double pow_10_2_n [] = { 1.0e+256 }; - /* Logarithm of the floating point radix. */ - -static int log2FLT_RADIX; - -/* Number of bits in a double's significand */ - -static int mantBits; - -/* Table of powers of 5**(2**n), up to 5**256 */ - -static mp_int pow5[9]; - -/* The smallest representable double */ - -static double tiny; - -/* The maximum number of digits to the left of the decimal point of a - * double. */ - -static int maxDigits; - -/* The maximum number of digits to the right of the decimal point in a - * double. */ - -static int minDigits; - -/* Number of mp_digit's needed to hold the significand of a double */ - -static int mantDIGIT; - -/* Static functions defined in this file */ - -static 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 )); +/* + * Static functions defined in this file. + */ +static double AbsoluteValue(double v, int *signum); +static int AccumulateDecimalDigit(unsigned, int, + Tcl_WideUInt *, mp_int *, int); +static double BignumToBiasedFrExp(mp_int *big, int* machexp); +static int GetIntegerTimesPower(double v, mp_int *r, int *e); +static double MakeHighPrecisionDouble(int signum, + mp_int *significand, int nSigDigs, int exponent); +static double MakeLowPrecisionDouble(int signum, + Tcl_WideUInt significand, int nSigDigs, + int exponent); +static double MakeNaN(int signum, Tcl_WideUInt tag); +static double Pow10TimesFrExp(int exponent, double fraction, + int *machexp); +static double RefineApproximation(double approx, + mp_int *exactSignificand, int exponent); +static double SafeLdExp(double fraction, int exponent); /* *---------------------------------------------------------------------- @@ -179,71 +151,64 @@ static double SafeLdExp _ANSI_ARGS_(( double fraction, int exponent )); * Returns a standard Tcl result. * * Side effects: - * 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. + * 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. * *---------------------------------------------------------------------- */ 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", +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 *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 */ { - enum State { INITIAL, SIGNUM, ZERO, ZERO_X, #ifdef TIP_114_FORMATS @@ -261,66 +226,69 @@ TclParseNumber( Tcl_Interp* interp, 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 */ + /* 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' */ + * 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 */ + 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 */ + 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 */ + 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 = 0; /* Last hexadecimal digit scanned; initialized + char d = 0; /* Last hexadecimal digit scanned; initialized * to avoid a compiler warning. */ int shift = 0; /* Amount to shift when accumulating binary */ #ifdef TIP_114_FORMATS int explicitOctal = 0; #endif +#define ALL_BITS (~(Tcl_WideUInt)0) +#define MOST_BITS (ALL_BITS >> 1) + /* - * Initialize string to start of the object's string rep if - * the caller didn't pass anything else. + * Initialize string to start of the object's string rep if the caller + * didn't pass anything else. */ - if ( string == NULL ) { - string = Tcl_GetStringFromObj( objPtr, NULL ); + if (string == NULL) { + string = TclGetString(objPtr); } p = string; len = length; acceptPoint = p; acceptLen = len; - while ( 1 ) { + while (1) { char c = len ? *p : '\0'; switch (state) { case INITIAL: /* - * Initial state. Acceptable characters are +, -, digits, - * period, I, N, and whitespace. + * Initial state. Acceptable characters are +, -, digits, period, + * I, N, and whitespace. */ + if (isspace(UCHAR(c))) { break; } else if (c == '+') { @@ -335,9 +303,10 @@ TclParseNumber( Tcl_Interp* interp, case SIGNUM: /* - * Scanned a leading + or -. Acceptable characters are - * digits, period, I, and N. + * Scanned a leading + or -. Acceptable characters are digits, + * period, I, and N. */ + if (c == '0') { if (flags & TCL_PARSE_DECIMAL_ONLY) { state = DECIMAL; @@ -372,12 +341,12 @@ TclParseNumber( Tcl_Interp* interp, 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'. + * 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; @@ -409,11 +378,11 @@ TclParseNumber( Tcl_Interp* interp, 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. + * 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; @@ -429,40 +398,39 @@ TclParseNumber( Tcl_Interp* interp, } else if (c >= '1' && c <= '7') { if (objPtr != NULL) { shift = 3 * (numTrailZeros + 1); - significandOverflow = - AccumulateDecimalDigit((unsigned)(c-'0'), - numTrailZeros, - &significandWide, - &significandBig, - significandOverflow); - + 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. + * 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); + octalSignificandWide); } } if (!octalSignificandOverflow) { - octalSignificandWide - = (octalSignificandWide << shift) + (c - '0'); + octalSignificandWide = + (octalSignificandWide << shift) + (c - '0'); } else { mp_mul_2d(&octalSignificandBig, shift, - &octalSignificandBig); + &octalSignificandBig); mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'), - &octalSignificandBig); + &octalSignificandBig); } } - if ( numSigDigs != 0 ) { - numSigDigs += ( numTrailZeros + 1 ); + if (numSigDigs != 0) { + numSigDigs += numTrailZeros+1; } else { numSigDigs = 1; } @@ -475,35 +443,42 @@ TclParseNumber( Tcl_Interp* interp, case BAD_OCTAL: #ifdef TIP_114_FORMATS if (explicitOctal) { - /* No forgiveness for bad digits in explicitly octal numbers */ + /* + * 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 */ + /* + * 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. + * 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') { ++numTrailZeros; state = BAD_OCTAL; break; } else if (isdigit(UCHAR(c))) { if (objPtr != NULL) { - significandOverflow = - AccumulateDecimalDigit((unsigned)(c-'0'), - numTrailZeros, - &significandWide, - &significandBig, - significandOverflow); + significandOverflow = AccumulateDecimalDigit( + (unsigned)(c-'0'), numTrailZeros, + &significandWide, &significandBig, + significandOverflow); } - if ( numSigDigs != 0 ) { - numSigDigs += ( numTrailZeros + 1 ); + if (numSigDigs != 0) { + numSigDigs += (numTrailZeros + 1); } else { numSigDigs = 1; } @@ -521,15 +496,17 @@ TclParseNumber( Tcl_Interp* interp, goto endgame; /* - * Scanned 0x. If state is HEXADECIMAL, scanned at least - * one character following the 0x. The only acceptable - * inputs are hexadecimal digits. + * 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') { @@ -550,25 +527,23 @@ TclParseNumber( Tcl_Interp* interp, 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. + * 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))) { + + if (significandWide != 0 && + (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || + significandWide > (~(Tcl_WideUInt)0 >> shift))) { significandOverflow = 1; TclBNInitBignumFromWideUInt(&significandBig, - significandWide); + significandWide); } } if (!significandOverflow) { - significandWide - = (significandWide << shift) + d; + significandWide = (significandWide << shift) + d; } else { - mp_mul_2d(&significandBig, shift, - &significandBig); - mp_add_d(&significandBig, (mp_digit) d, - &significandBig); + mp_mul_2d(&significandBig, shift, &significandBig); + mp_add_d(&significandBig, (mp_digit) d, &significandBig); } } numTrailZeros = 0; @@ -593,25 +568,23 @@ TclParseNumber( Tcl_Interp* interp, 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. + * 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))) { + + if (significandWide != 0 && + (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || + significandWide > (~(Tcl_WideUInt)0 >> shift))) { significandOverflow = 1; TclBNInitBignumFromWideUInt(&significandBig, - significandWide); + significandWide); } } if (!significandOverflow) { - significandWide - = (significandWide << shift) + 1; + significandWide = (significandWide << shift) + 1; } else { - mp_mul_2d(&significandBig, shift, - &significandBig); - mp_add_d(&significandBig, (mp_digit) 1, - &significandBig); + mp_mul_2d(&significandBig, shift, &significandBig); + mp_add_d(&significandBig, (mp_digit) 1, &significandBig); } } numTrailZeros = 0; @@ -621,9 +594,10 @@ TclParseNumber( Tcl_Interp* interp, case DECIMAL: /* - * Scanned an optional + or - followed by a string of - * decimal digits. + * Scanned an optional + or - followed by a string of decimal + * digits. */ + #ifdef KILL_OCTAL decimal: #endif @@ -636,14 +610,12 @@ TclParseNumber( Tcl_Interp* interp, break; } else if (isdigit(UCHAR(c))) { if (objPtr != NULL) { - significandOverflow = - AccumulateDecimalDigit((unsigned)(c - '0'), - numTrailZeros, - &significandWide, - &significandBig, - significandOverflow); + significandOverflow = AccumulateDecimalDigit( + (unsigned)(c - '0'), numTrailZeros, + &significandWide, &significandBig, + significandOverflow); } - numSigDigs += ( numTrailZeros + 1 ); + numSigDigs += numTrailZeros+1; numTrailZeros = 0; state = DECIMAL; break; @@ -659,11 +631,11 @@ TclParseNumber( Tcl_Interp* interp, goto endgame; /* - * 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. + * 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; @@ -673,6 +645,7 @@ TclParseNumber( Tcl_Interp* interp, break; } /* FALLTHROUGH */ + case LEADING_RADIX_POINT: if (c == '0') { ++numDigitsAfterDp; @@ -682,15 +655,13 @@ TclParseNumber( Tcl_Interp* interp, } else if (isdigit(UCHAR(c))) { ++numDigitsAfterDp; if (objPtr != NULL) { - significandOverflow = - AccumulateDecimalDigit((unsigned)(c-'0'), - numTrailZeros, - &significandWide, - &significandBig, - significandOverflow); + significandOverflow = AccumulateDecimalDigit( + (unsigned)(c-'0'), numTrailZeros, + &significandWide, &significandBig, + significandOverflow); } - if ( numSigDigs != 0 ) { - numSigDigs += ( numTrailZeros + 1 ); + if (numSigDigs != 0) { + numSigDigs += numTrailZeros+1; } else { numSigDigs = 1; } @@ -702,10 +673,11 @@ TclParseNumber( Tcl_Interp* interp, 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. + * 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; @@ -718,9 +690,10 @@ TclParseNumber( Tcl_Interp* interp, case EXPONENT_SIGNUM: /* - * Found the E at the start of the exponent, followed by - * a sign character. + * Found the E at the start of the exponent, followed by a sign + * character. */ + if (isdigit(UCHAR(c))) { exponent = c - '0'; state = EXPONENT; @@ -730,10 +703,10 @@ TclParseNumber( Tcl_Interp* interp, case EXPONENT: /* - * Found an exponent with at least one digit. - * Accumulate it, making sure to hard-pin it to LONG_MAX - * on overflow. + * 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; @@ -749,18 +722,18 @@ TclParseNumber( Tcl_Interp* interp, goto endgame; /* - * Parse out INFINITY by simply spelling it out. - * INF is accepted as an abbreviation; other prefices are - * not. + * Parse out INFINITY by simply spelling it out. INF is accepted + * as an abbreviation; other prefices are not. */ + case sI: - if ( c == 'n' || c == 'N' ) { + if (c == 'n' || c == 'N') { state = sIN; break; } goto endgame; case sIN: - if ( c == 'f' || c == 'F' ) { + if (c == 'f' || c == 'F') { state = sINF; break; } @@ -769,31 +742,31 @@ TclParseNumber( Tcl_Interp* interp, acceptState = state; acceptPoint = p; acceptLen = len; - if ( c == 'i' || c == 'I' ) { + if (c == 'i' || c == 'I') { state = sINFI; break; } goto endgame; case sINFI: - if ( c == 'n' || c == 'N' ) { + if (c == 'n' || c == 'N') { state = sINFIN; break; } goto endgame; case sINFIN: - if ( c == 'i' || c == 'I' ) { + if (c == 'i' || c == 'I') { state = sINFINI; break; } goto endgame; case sINFINI: - if ( c == 't' || c == 'T' ) { + if (c == 't' || c == 'T') { state = sINFINIT; break; } goto endgame; case sINFINIT: - if ( c == 'y' || c == 'Y' ) { + if (c == 'y' || c == 'Y') { state = sINFINITY; break; } @@ -804,13 +777,13 @@ TclParseNumber( Tcl_Interp* interp, */ #ifdef IEEE_FLOATING_POINT case sN: - if ( c == 'a' || c == 'A' ) { + if (c == 'a' || c == 'A') { state = sNA; break; } goto endgame; case sNA: - if ( c == 'n' || c == 'N' ) { + if (c == 'n' || c == 'N') { state = sNAN; break; } @@ -819,7 +792,7 @@ TclParseNumber( Tcl_Interp* interp, acceptState = state; acceptPoint = p; acceptLen = len; - if ( c == '(' ) { + if (c == '(') { state = sNANPAREN; break; } @@ -829,21 +802,21 @@ TclParseNumber( Tcl_Interp* interp, * Parse NaN(hexdigits) */ case sNANHEX: - if ( c == ')' ) { + if (c == ')') { state = sNANFINISH; break; } /* FALLTHROUGH */ case sNANPAREN: - if ( isspace(UCHAR(c)) ) { + if (isspace(UCHAR(c))) { break; } - if ( numSigDigs < 13 ) { - if ( c >= '0' && c <= '9' ) { + if (numSigDigs < 13) { + if (c >= '0' && c <= '9') { d = c - '0'; - } else if ( c >= 'a' && c <= 'f' ) { + } else if (c >= 'a' && c <= 'f') { d = 10 + c - 'a'; - } else if ( c >= 'A' && c <= 'F' ) { + } else if (c >= 'A' && c <= 'F') { d = 10 + c - 'A'; } significandWide = (significandWide << 4) + d; @@ -853,6 +826,7 @@ TclParseNumber( Tcl_Interp* interp, goto endgame; case sNANFINISH: #endif + case sINFINITY: acceptState = state; acceptPoint = p; @@ -865,7 +839,9 @@ TclParseNumber( Tcl_Interp* interp, endgame: - /* Back up to the last accepting state in the lexer */ + /* + * Back up to the last accepting state in the lexer. + */ if (acceptState == INITIAL) { status = TCL_ERROR; @@ -873,7 +849,9 @@ TclParseNumber( Tcl_Interp* interp, p = acceptPoint; len = acceptLen; - /* Skip past trailing whitespace */ + /* + * Skip past trailing whitespace. + */ if (endPtrPtr != NULL) { *endPtrPtr = p; @@ -884,17 +862,21 @@ TclParseNumber( Tcl_Interp* interp, --len; } - /* Determine whether a partial string is acceptable. */ + /* + * Determine whether a partial string is acceptable. + */ if (endPtrPtr == NULL && len != 0 && *p != '\0') { status = TCL_ERROR; } - /* Generate and store the appropriate internal rep */ + /* + * Generate and store the appropriate internal rep. + */ if (status == TCL_OK && objPtr != NULL) { - if ( acceptState != INITIAL ) { - TclFreeIntRep( objPtr ); + if (acceptState != INITIAL) { + TclFreeIntRep(objPtr); } switch (acceptState) { @@ -927,74 +909,69 @@ TclParseNumber( Tcl_Interp* interp, #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); - } + if (!significandOverflow && significandWide != 0 && + (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || + significandWide > (MOST_BITS + signum) >> shift)) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, significandWide); } if (shift) { - if ( !significandOverflow ) { + if (!significandOverflow) { significandWide <<= shift; } else { - mp_mul_2d( &significandBig, shift, &significandBig ); + mp_mul_2d(&significandBig, shift, &significandBig); } } goto returnInteger; #endif + case HEXADECIMAL: - /* Returning a hex integer. Final scaling step */ + /* + * 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 (!significandOverflow && significandWide !=0 && + (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || + significandWide > (MOST_BITS + signum) >> shift)) { + significandOverflow = 1; + TclBNInitBignumFromWideUInt(&significandBig, significandWide); } if (shift) { - if ( !significandOverflow ) { + if (!significandOverflow) { significandWide <<= shift; } else { - mp_mul_2d( &significandBig, shift, &significandBig ); + mp_mul_2d(&significandBig, shift, &significandBig); } } goto returnInteger; case OCTAL: - /* Returning an octal integer. Final scaling step */ + /* + * 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 (!octalSignificandOverflow && octalSignificandWide != 0 && + (shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || + octalSignificandWide > (MOST_BITS + signum) >> shift)) { + octalSignificandOverflow = 1; + TclBNInitBignumFromWideUInt(&octalSignificandBig, + octalSignificandWide); } - if ( shift ) { - if ( !octalSignificandOverflow ) { + if (shift) { + if (!octalSignificandOverflow) { octalSignificandWide <<= shift; } else { - mp_mul_2d( &octalSignificandBig, shift, - &octalSignificandBig ); + 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)) { + if (octalSignificandWide <= (MOST_BITS + signum)) { objPtr->typePtr = &tclWideIntType; if (signum) { objPtr->internalRep.wideValue = @@ -1007,7 +984,7 @@ TclParseNumber( Tcl_Interp* interp, } #endif TclBNInitBignumFromWideUInt(&octalSignificandBig, - octalSignificandWide); + octalSignificandWide); octalSignificandOverflow = 1; } else { objPtr->typePtr = &tclIntType; @@ -1030,24 +1007,18 @@ TclParseNumber( Tcl_Interp* interp, case ZERO: case DECIMAL: - significandOverflow = - AccumulateDecimalDigit( 0, numTrailZeros-1, - &significandWide, &significandBig, - significandOverflow ); - if (!significandOverflow - && (significandWide - > (((~(Tcl_WideUInt)0) >> 1) + signum))) { + significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1, + &significandWide, &significandBig, significandOverflow); + if (!significandOverflow && (significandWide > MOST_BITS+signum)) { significandOverflow = 1; - TclBNInitBignumFromWideUInt(&significandBig, - significandWide); + TclBNInitBignumFromWideUInt(&significandBig, significandWide); } - returnInteger: + returnInteger: if (!significandOverflow) { if (significandWide > (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { #ifndef NO_WIDE_TYPE - if (significandWide - <= (((~(Tcl_WideUInt)0) >> 1) + signum)) { + if (significandWide <= MOST_BITS+signum) { objPtr->typePtr = &tclWideIntType; if (signum) { objPtr->internalRep.wideValue = @@ -1060,7 +1031,7 @@ TclParseNumber( Tcl_Interp* interp, } #endif TclBNInitBignumFromWideUInt(&significandBig, - significandWide); + significandWide); significandOverflow = 1; } else { objPtr->typePtr = &tclIntType; @@ -1085,40 +1056,31 @@ TclParseNumber( Tcl_Interp* interp, 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 + * 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 ) { + if (exponentSignum) { exponent = - exponent; } - if ( !significandOverflow ) { - objPtr->internalRep.doubleValue = - MakeLowPrecisionDouble( signum, - significandWide, - numSigDigs, - ( numTrailZeros - + exponent - - numDigitsAfterDp ) ); + if (!significandOverflow) { + objPtr->internalRep.doubleValue = MakeLowPrecisionDouble( + signum, significandWide, numSigDigs, + (numTrailZeros + exponent - numDigitsAfterDp)); } else { - objPtr->internalRep.doubleValue = - MakeHighPrecisionDouble( signum, - &significandBig, - numSigDigs, - ( numTrailZeros - + exponent - - numDigitsAfterDp ) ); + objPtr->internalRep.doubleValue = MakeHighPrecisionDouble( + signum, &significandBig, numSigDigs, + (numTrailZeros + exponent - numDigitsAfterDp)); } break; case sINF: case sINFINITY: - if ( signum ) { + if (signum) { objPtr->internalRep.doubleValue = -HUGE_VAL; } else { objPtr->internalRep.doubleValue = HUGE_VAL; @@ -1128,32 +1090,34 @@ TclParseNumber( Tcl_Interp* interp, case sNAN: case sNANFINISH: - objPtr->internalRep.doubleValue - = MakeNaN( signum, significandWide ); + objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide); objPtr->typePtr = &tclDoubleType; break; - } } - /* Format an error message when an invalid number is encountered. */ + /* + * 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 ); + 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 ); + Tcl_SetObjResult(interp, msg); } } - /* Free memory */ + /* + * Free memory. + */ if (octalSignificandOverflow) { mp_clear(&octalSignificandBig); @@ -1172,8 +1136,8 @@ TclParseNumber( Tcl_Interp* interp, * 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. + * 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. @@ -1182,75 +1146,73 @@ TclParseNumber( Tcl_Interp* interp, */ 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. */ +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 */ + /* + * 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; - } - } + if (!bignumFlag && *wideRepPtr!=0 && ((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. */ + /* + * Multiply the number by 10**numZeros+1 and add in the new digit. + */ if (!bignumFlag) { - - /* Wide multiplication */ + /* + * 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 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 { - /* * 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). + * 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) { + 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); + mp_mul(bignumRepPtr, pow5+i, bignumRepPtr); } } while (n >= 256) { - mp_mul (bignumRepPtr, pow5+8, bignumRepPtr); + mp_mul(bignumRepPtr, pow5+8, bignumRepPtr); n -= 256; } - mp_mul_2d (bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr); + mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr); } return bignumFlag; @@ -1266,24 +1228,21 @@ AccumulateDecimalDigit( unsigned digit, * 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. + * 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 */ +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 */ @@ -1292,55 +1251,54 @@ MakeLowPrecisionDouble( int signum, * 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. + * precision of double. Double-rounding introduces gratuitous errors of 1 + * ulp, so we need to change rounding mode to 53-bits. */ #if defined(__GNUC__) && defined(__i386) fpu_control_t roundTo53Bits = 0x027f; fpu_control_t oldRoundingMode; - _FPU_GETCW( oldRoundingMode ); - _FPU_SETCW( roundTo53Bits ); + _FPU_GETCW(oldRoundingMode); + _FPU_SETCW(roundTo53Bits); #endif - /* Test for the easy cases */ - - if ( numSigDigs <= DBL_DIG ) { - if ( exponent >= 0 ) { - if ( exponent <= mmaxpow ) { + /* + * Test for the easy cases. + */ + if (numSigDigs <= DBL_DIG) { + if (exponent >= 0) { + if (exponent <= mmaxpow) { /* - * The significand is an exact integer, and so is - * 10**exponent. The product will be correct to within - * 1/2 ulp without special handling. + * The significand is an exact integer, and so is + * 10**exponent. The product will be correct to within 1/2 ulp + * without special handling. */ retval = (double)(Tcl_WideInt)significand * pow10[ exponent ]; goto returnValue; - } else { int diff = DBL_DIG - numSigDigs; - if ( exponent-diff <= mmaxpow ) { - + 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. + * significand*10**diff, so we can still compute the value + * with only one roundoff. */ - volatile double factor - = (double)(Tcl_WideInt)significand * pow10[diff]; + + volatile double factor = + (double)(Tcl_WideInt)significand * pow10[diff]; retval = factor * pow10[exponent-diff]; goto returnValue; } } } else { - if ( exponent >= -mmaxpow ) { - + 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. + * significand. Compute the result by one division, again with + * only one rounding. */ retval = (double)(Tcl_WideInt)significand / pow10[-exponent]; @@ -1350,26 +1308,29 @@ MakeLowPrecisionDouble( int signum, } /* - * All the easy cases have failed. Promote ths significand - * to bignum and call MakeHighPrecisionDouble to do it the hard way. + * 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 ); + TclBNInitBignumFromWideUInt(&significandBig, significand); + retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs, + exponent); - /* Come here to return the computed value */ + /* + * Come here to return the computed value. + */ returnValue: - - if ( signum ) { + if (signum) { retval = -retval; } - /* On gcc on x86, restore the floating point mode word. */ + /* + * On gcc on x86, restore the floating point mode word. + */ #if defined(__GNUC__) && defined(__i386) - _FPU_SETCW( oldRoundingMode ); + _FPU_SETCW(oldRoundingMode); #endif return retval; @@ -1385,25 +1346,21 @@ MakeLowPrecisionDouble( int signum, * 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. + * 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 */ +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 */ @@ -1411,69 +1368,73 @@ MakeHighPrecisionDouble( int signum, * 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. + * precision of double. Double-rounding introduces gratuitous errors of 1 + * ulp, so we need to change rounding mode to 53-bits. */ #if defined(__GNUC__) && defined(__i386) fpu_control_t roundTo53Bits = 0x027f; fpu_control_t oldRoundingMode; - _FPU_GETCW( oldRoundingMode ); - _FPU_SETCW( roundTo53Bits ); + _FPU_GETCW(oldRoundingMode); + _FPU_SETCW(roundTo53Bits); #endif - /* Quick checks for over/underflow */ + /* + * Quick checks for over/underflow. + */ - if ( numSigDigs + exponent - 1 > maxDigits ) { + if (numSigDigs+exponent-1 > maxDigits) { retval = HUGE_VAL; goto returnValue; } - if ( numSigDigs + exponent - 1 < minDigits ) { + if (numSigDigs+exponent-1 < minDigits) { retval = 0; goto returnValue; } /* - * 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. + * 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. */ - retval = BignumToBiasedFrExp( significand, &machexp ); - retval = Pow10TimesFrExp( exponent, retval, &machexp ); - if ( machexp > DBL_MAX_EXP * log2FLT_RADIX ) { + 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 = SafeLdExp(retval, machexp); + if (retval < tiny) { retval = tiny; } /* - * Refine the result twice. (The second refinement should be - * necessary only if the best approximation is a power of 2 - * minus 1/2 ulp). + * Refine the result twice. (The second refinement should be necessary + * only if the best approximation is a power of 2 minus 1/2 ulp). */ - retval = RefineApproximation( retval, significand, exponent ); - retval = RefineApproximation( retval, significand, exponent ); + retval = RefineApproximation(retval, significand, exponent); + retval = RefineApproximation(retval, significand, exponent); - /* Come here to return the computed value */ + /* + * Come here to return the computed value. + */ returnValue: - if ( signum ) { + if (signum) { retval = -retval; } - /* On gcc on x86, restore the floating point mode word. */ + /* + * On gcc on x86, restore the floating point mode word. + */ #if defined(__GNUC__) && defined(__i386) - _FPU_SETCW( oldRoundingMode ); + _FPU_SETCW(oldRoundingMode); #endif return retval; } @@ -1483,18 +1444,18 @@ MakeHighPrecisionDouble( int signum, * * MakeNaN -- * - * Makes a "Not a Number" given a set of bits to put in the - * tag bits + * Makes a "Not a Number" given a set of bits to put in the tag bits * - * Note that a signalling NaN is never returned. + * 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 */ +MakeNaN( + int signum, /* Sign bit (1=negative, 0=nonnegative */ + Tcl_WideUInt tags) /* Tag bits to put in the NaN */ { union { Tcl_WideUInt iv; @@ -1502,8 +1463,8 @@ MakeNaN( int signum, /* Sign bit (1=negative, 0=nonnegative */ } theNaN; theNaN.iv = tags; - theNaN.iv &= ( ((Tcl_WideUInt) 1) << 51 ) - 1; - if ( signum ) { + 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; @@ -1518,10 +1479,9 @@ MakeNaN( int signum, /* Sign bit (1=negative, 0=nonnegative */ * * 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.) + * 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. @@ -1530,109 +1490,104 @@ MakeNaN( int signum, /* Sign bit (1=negative, 0=nonnegative */ */ static double -RefineApproximation( double approxResult, - /* Approximate result of conversion */ - mp_int* exactSignificand, - /* Integer significand */ - int exponent ) - /* Power of 10 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. */ + int M2, M5; /* Powers of 2 and of 5 needed to put the + * decimal and binary numbers over a common + * denominator. */ double significand; /* Sigificand of the binary number */ int binExponent; /* Exponent of the binary number */ - int msb; /* Most significant bit position of an * intermediate result */ int nDigits; /* Number of mp_digit's in an intermediate * result */ - mp_int twoMv; /* Approx binary value expressed as an - * exact integer scaled by the multiplier 2M */ - mp_int twoMd; /* Exact decimal value expressed as an - * exact integer scaled by the multiplier 2M */ + 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 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. */ + double minincr; /* Lower bound on the absolute value of the + * correction term. */ int i; /* - * 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. - * 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 + * 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_copy( &twoMd, exactSignificand ); - 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 @@ -1642,49 +1597,48 @@ RefineApproximation( 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; } @@ -1697,7 +1651,7 @@ RefineApproximation( double approxResult, * * Results: * Returns the position of the character in the string after which the - * decimal point should appear. Since the string contains only + * 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. * @@ -1706,24 +1660,26 @@ RefineApproximation( double approxResult, * the sign of the number. * *---------------------------------------------------------------------- + */ int -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. */ +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. */ { 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. */ + 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. */ @@ -1740,7 +1696,9 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, char c; int i, k, n; - /* Split the number into absolute value and signum. */ + /* + * Split the number into absolute value and signum. + */ v = AbsoluteValue(v, signum); @@ -1748,7 +1706,7 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, * Handle zero specially. */ - if ( v == 0.0 ) { + if (v == 0.0) { *string++ = '0'; *string++ = '\0'; return 1; @@ -1757,8 +1715,8 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, /* * 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. + * and r is as small as possible. Also determine whether the significand + * is the smallest possible. */ smallestSig = GetIntegerTimesPower(v, &r, &e); @@ -1807,7 +1765,7 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, 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. + * exponent. mplus and mminus will both be 1. */ rfac2 = 1; @@ -1892,7 +1850,7 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, /* * 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. + * 10**k. mplus and mminus give the rounding limits. */ for (;;) { @@ -1916,8 +1874,8 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, } else { tc2= (tc2 > 0); } - if ( ! tc1 ) { - if ( !tc2 ) { + if (!tc1) { + if (!tc2) { *string++ = '0' + i; } else { c = (char) (i + '1'); @@ -1966,8 +1924,9 @@ TclDoubleDigits( char * string, /* Buffer in which to store the result, */ static double -AbsoluteValue (double v, /* Number to split */ - int* signum) /* (Output) Sign of the number 1=-, 0=+ */ +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. @@ -2005,32 +1964,28 @@ AbsoluteValue (double v, /* Number to split */ * * GetIntegerTimesPower -- * - * Converts a floating point number to an exact integer times a - * power of the floating point radix. + * 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. + * 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*/ +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; + double a, f; + int e, i, n; /* * Develop f and e such that v = f * FLT_RADIX**e, with @@ -2122,10 +2077,10 @@ TclInitDoubleConversion(void) * 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)); + 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; @@ -2138,8 +2093,8 @@ TclInitDoubleConversion(void) * decimal digits that represents. */ - 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; @@ -2150,8 +2105,8 @@ TclInitDoubleConversion(void) * in a double. */ - 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; @@ -2161,14 +2116,16 @@ TclInitDoubleConversion(void) d *= 10.0; } - /* Initialize a table of large powers of five. */ + /* + * Initialize a table of large powers of five. + */ - 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); } /* @@ -2178,14 +2135,13 @@ TclInitDoubleConversion(void) * the significand of a double. */ - tiny = SafeLdExp( 1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits ); + 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.)); + + 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.)); } /* @@ -2208,9 +2164,10 @@ void TclFinalizeDoubleConversion() { int i; - Tcl_Free ((char*)pow10_wide); - 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); } } @@ -2219,28 +2176,32 @@ TclFinalizeDoubleConversion() * * TclInitBignumFromDouble -- * - * Extracts the integer part of a double and converts it to - * an arbitrary precision integer. + * 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. + * 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 */ +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 */ + /* + * Infinite values can't convert to bignum. + */ + if (TclIsInfinite(d)) { if (interp != NULL) { char *s = "integer value too large to represent"; @@ -2249,6 +2210,7 @@ TclInitBignumFromDouble(Tcl_Interp *interp, /* For error message */ } return TCL_ERROR; } + fract = frexp(d,&expt); if (expt <= 0) { mp_init(b); @@ -2256,6 +2218,7 @@ TclInitBignumFromDouble(Tcl_Interp *interp, /* For error message */ } 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); @@ -2275,14 +2238,15 @@ TclInitBignumFromDouble(Tcl_Interp *interp, /* For error message */ * number. * * Results: - * Returns the converted number. Sets errno to ERANGE if the number is + * 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; @@ -2344,7 +2308,8 @@ TclBignumToDouble(mp_int *a) /* Integer to convert. */ } double -TclCeil(mp_int *a) /* Integer to convert. */ +TclCeil( + mp_int *a) /* Integer to convert. */ { double r = 0.0; mp_int b; @@ -2386,7 +2351,8 @@ TclCeil(mp_int *a) /* Integer to convert. */ } double -TclFloor(mp_int *a) /* Integer to convert. */ +TclFloor( + mp_int *a) /* Integer to convert. */ { double r = 0.0; mp_int b; @@ -2425,11 +2391,11 @@ TclFloor(mp_int *a) /* Integer to convert. */ * * 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. + * 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. @@ -2441,10 +2407,9 @@ TclFloor(mp_int *a) /* Integer to convert. */ */ static double -BignumToBiasedFrExp( mp_int* a, - /* Integer to convert */ - int* machexp ) - /* Power of two */ +BignumToBiasedFrExp( + mp_int *a, /* Integer to convert */ + int *machexp) /* Power of two */ { mp_int b; int bits; @@ -2452,36 +2417,38 @@ BignumToBiasedFrExp( 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 ); + 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 ); + 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); } - /* 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); - /* Return the result with the appropriate sign. */ + /* + * Return the result with the appropriate sign. + */ *machexp = bits - mantBits + 2; - if ( a->sign == MP_ZPOS ) { - return r; - } else { - return -r; - } + return ((a->sign == MP_ZPOS) ? r : -r); } /* @@ -2496,47 +2463,48 @@ BignumToBiasedFrExp( mp_int* a, * Returns the significand of the result. * * Side effects: - * Overwrites the 'machexp' parameter with the exponent of the - * result. + * 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. + * 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. */ +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 */ + if (exponent > 0) { + /* + * Multiply by 10**exponent + */ - retval = frexp( retval * pow10[ exponent & 0xf ], &j ); + 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 ); + 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 */ + } else if (exponent < 0) { + /* + * Divide by 10**-exponent + */ - retval = frexp( retval / pow10[ (-exponent) & 0xf ], &j ); + 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 ); + for (i=4; i<9; ++i) { + if ((-exponent) & (1<<i)) { + retval = frexp(retval / pow_10_2_n[i], &j); expt += j; } } @@ -2544,7 +2512,6 @@ Pow10TimesFrExp( int exponent, /* Power of 10 to multiply by */ *machexp = expt; return retval; - } /* @@ -2566,10 +2533,13 @@ Pow10TimesFrExp( int exponent, /* Power of 10 to multiply by */ */ 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); @@ -2598,8 +2568,9 @@ SafeLdExp(double fract, int expt) */ 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"); @@ -2626,3 +2597,11 @@ TclFormatNaN(double value, /* The Not-a-Number to format. */ } #endif /* IEEE_FLOATING_POINT */ } + +/* + * Local Variables: + * mode: c + * c-basic-offset: 4 + * fill-column: 78 + * End: + */ |