tclStrToD.c Go to the documentation of this file. /* *---------------------------------------------------------------------- * * 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. * * 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.32 2007/12/13 15:23:20 dgp Exp $ * *---------------------------------------------------------------------- */ #include <tclInt.h> #include <stdio.h> #include <stdlib.h> #include <float.h> #include <limits.h> #include <math.h> #include <ctype.h> #include <tommath.h> /* * Define KILL_OCTAL to suppress interpretation of numbers with leading zero * as octal. (Ceterum censeo: numeros octonarios delendos esse.) */ #undef KILL_OCTAL /* * 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, 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__ __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 /* * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN. * Everyone else uses 7ff8000000000000. (Why, HP, why?) */ #ifdef __hppa # define NAN_START 0x7ff4 # define NAN_MASK (((Tcl_WideUInt) 1) << 50) #else # define NAN_START 0x7ff8 # define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif /* * Constants used by this file (most of which are only ever calculated at * runtime). */ static int maxpow10_wide; /* The powers of ten that can be represented * exactly as wide integers. */ static Tcl_WideUInt *pow10_wide; #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'. */ 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, 1.0e+8, 1.0e+16, 1.0e+32, 1.0e+64, 1.0e+128, 1.0e+256 }; static int n770_fp; /* Flag is 1 on Nokia N770 floating point. * Nokia's floating point has the words * reversed: if big-endian is 7654 3210, * and little-endian is 0123 4567, * then Nokia's FP is 4567 0123; * little-endian within the 32-bit words * but big-endian between them. */ /* * 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 Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w); 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);  /* *---------------------------------------------------------------------- * * TclParseNumber -- * * Scans bytes, interpreted as characters in Tcl's internal encoding, and * parses the longest prefix that is the string representation of a * number in a format recognized by Tcl. * * The arguments bytes, numBytes, and objPtr are the inputs which * determine the string to be parsed. If bytes is non-NULL, it points to * the first byte to be scanned. If bytes is NULL, then objPtr must be * non-NULL, and the string representation of objPtr will be scanned * (generated first, if necessary). The numBytes argument determines the * number of bytes to be scanned. If numBytes is negative, the first NUL * byte encountered will terminate the scan. If numBytes is non-negative, * then no more than numBytes bytes will be scanned. * * The argument flags is an input that controls the numeric formats * recognized by the parser. The flag bits are: * * - TCL_PARSE_INTEGER_ONLY: accept only integer values; reject * strings that denote floating point values (or accept only the * leading portion of them that are integer values). * - TCL_PARSE_SCAN_PREFIXES: ignore the prefixes 0b and 0o that are * not part of the [scan] command's vocabulary. Use only in * combination with TCL_PARSE_INTEGER_ONLY. * - TCL_PARSE_OCTAL_ONLY: parse only in the octal format, whether * or not a prefix is present that would lead to octal parsing. * Use only in combination with TCL_PARSE_INTEGER_ONLY. * - TCL_PARSE_HEXADECIMAL_ONLY: parse only in the hexadecimal format, * whether or not a prefix is present that would lead to * hexadecimal parsing. Use only in combination with * TCL_PARSE_INTEGER_ONLY. * - TCL_PARSE_DECIMAL_ONLY: parse only in the decimal format, no * matter whether a 0 prefix would normally force a different * base. * - TCL_PARSE_NO_WHITESPACE: reject any leading/trailing whitespace * * The arguments interp and expected are inputs that control error * message generation. If interp is NULL, no error message will be * generated. If interp is non-NULL, then expected must also be non-NULL. * When TCL_ERROR is returned, an error message will be left in the * result of interp, and the expected argument will appear in the error * message as the thing TclParseNumber expected, but failed to find in * the string. * * The arguments objPtr and endPtrPtr as well as the return code are the * outputs. * * When the parser cannot find any prefix of the string that matches a * format it is looking for, TCL_ERROR is returned and an error message * may be generated and returned as described above. The contents of * objPtr will not be changed. If endPtrPtr is non-NULL, a pointer to the * character in the string that terminated the scan will be written to * *endPtrPtr. * * When the parser determines that the entire string matches a format it * is looking for, TCL_OK is returned, and if objPtr is non-NULL, then * the internal rep and Tcl_ObjType of objPtr are set to the "canonical" * numeric value that matches the scanned string. If endPtrPtr is not * NULL, a pointer to the end of the string will be written to *endPtrPtr * (that is, either bytes+numBytes or a pointer to a terminating NUL * byte). * * When the parser determines that a partial string matches a format it * is looking for, the value of endPtrPtr determines what happens: * * - If endPtrPtr is NULL, then TCL_ERROR is returned, with error message * generation as above. * * - If endPtrPtr is non-NULL, then TCL_OK is returned and objPtr * internals are set as above. Also, a pointer to the first * character following the parsed numeric string is written to * *endPtrPtr. * * In some cases where the string being scanned is the string rep of * objPtr, this routine can leave objPtr in an inconsistent state where * its string rep and its internal rep do not agree. In these cases the * internal rep will be in agreement with only some substring of the * string rep. This might happen if the caller passes in a non-NULL bytes * value that points somewhere into the string rep. It might happen if * the caller passes in a numBytes value that limits the scan to only a * prefix of the string rep. Or it might happen if a non-NULL value of * endPtrPtr permits a TCL_OK return from only a partial string match. It * is the responsibility of the caller to detect and correct such * inconsistencies when they can and do arise. * * Results: * Returns a standard Tcl result. * * Side effects: * The string representaton of objPtr may be generated. * * The internal representation and Tcl_ObjType of objPtr may be changed. * This may involve allocation and/or freeing of memory. * *---------------------------------------------------------------------- */ int TclParseNumber( Tcl_Interp *interp, /* Used for error reporting. May be NULL. */ Tcl_Obj *objPtr, /* Object to receive the internal rep. */ const char *expected, /* Description of the type of number the * caller expects to be able to parse * ("integer", "boolean value", etc.). */ const char *bytes, /* Pointer to the start of the string to * scan. */ int numBytes, /* Maximum number of bytes to scan, see * above. */ const char **endPtrPtr, /* Place to store pointer to the character * that terminated the scan. */ int flags) /* Flags governing the parse. */ { enum State { INITIAL, SIGNUM, ZERO, ZERO_X, ZERO_O, ZERO_B, BINARY, 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 = 0; /* Last hexadecimal digit scanned; initialized * to avoid a compiler warning. */ int shift = 0; /* Amount to shift when accumulating binary */ int explicitOctal = 0; #define ALL_BITS (~(Tcl_WideUInt)0) #define MOST_BITS (ALL_BITS >> 1) /* * Initialize bytes to start of the object's string rep if the caller * didn't pass anything else. */ if (bytes == NULL) { bytes = TclGetString(objPtr); } p = bytes; len = numBytes; acceptPoint = p; acceptLen = len; while (1) { char c = len ? *p : '\0'; switch (state) { case INITIAL: /* * Initial state. Acceptable characters are +, -, digits, period, * I, N, and whitespace. */ if (isspace(UCHAR(c))) { if (flags & TCL_PARSE_NO_WHITESPACE) { goto endgame; } break; } else if (c == '+') { state = SIGNUM; break; } else if (c == '-') { signum = 1; state = SIGNUM; break; } /* FALLTHROUGH */ 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; 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; } 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 /* FALLTHROUGH */ 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; /* FALLTHROUGH */ case ZERO_O: 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) && (((size_t)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 */ case BAD_OCTAL: if (explicitOctal) { /* * No forgiveness for bad digits in explicitly octal numbers. */ goto endgame; } 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') { ++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 { 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 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || significandWide > (~(Tcl_WideUInt)0 >> shift))) { significandOverflow = 1; TclBNInitBignumFromWideUInt(&significandBig, significandWide); } } if (!significandOverflow) { significandWide = (significandWide << shift) + d; } else { mp_mul_2d(&significandBig, shift, &significandBig); mp_add_d(&significandBig, (mp_digit) d, &significandBig); } } numTrailZeros = 0; state = HEXADECIMAL; break; case BINARY: acceptState = state; acceptPoint = p; acceptLen = len; case ZERO_B: if (c == '0') { ++numTrailZeros; state = BINARY; break; } else if (c != '1') { goto endgame; } 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 && ((size_t)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; 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; /* * 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 { numSigDigs = 1; } numTrailZeros = 0; state = FRACTION; break; } goto endgame; 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 */ 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; 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; } state = EXPONENT; break; } goto endgame; /* * 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') { 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; /* * 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; } goto endgame; case sNAN: acceptState = state; acceptPoint = p; acceptLen = len; if (c == '(') { state = sNANPAREN; break; } goto endgame; /* * 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: if (acceptState == INITIAL) { /* * No numeric string at all found. */ status = TCL_ERROR; if (endPtrPtr != NULL) { *endPtrPtr = p; } } else { /* * Back up to the last accepting state in the lexer. */ p = acceptPoint; len = acceptLen; if (!(flags & TCL_PARSE_NO_WHITESPACE)) { /* * Accept trailing whitespace. */ while (len != 0 && isspace(UCHAR(*p))) { ++p; --len; } } if (endPtrPtr == NULL) { if ((len != 0) && ((numBytes > 0) || (*p != '\0'))) { status = TCL_ERROR; } } else { *endPtrPtr = p; } } /* * Generate and store the appropriate internal rep. */ if (status == TCL_OK && objPtr != NULL) { TclFreeIntRep(objPtr); switch (acceptState) { case SIGNUM: case BAD_OCTAL: case ZERO_X: case ZERO_O: case ZERO_B: 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: Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'", acceptState, bytes); case BINARY: shift = numTrailZeros; if (!significandOverflow && significandWide != 0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || significandWide > (MOST_BITS + signum) >> shift)) { significandOverflow = 1; TclBNInitBignumFromWideUInt(&significandBig, significandWide); } if (shift) { if (!significandOverflow) { significandWide <<= shift; } else { mp_mul_2d(&significandBig, shift, &significandBig); } } goto returnInteger; case HEXADECIMAL: /* * Returning a hex integer. Final scaling step. */ shift = 4 * numTrailZeros; if (!significandOverflow && significandWide !=0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || significandWide > (MOST_BITS + 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 && octalSignificandWide != 0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || octalSignificandWide > (MOST_BITS + 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 <= (MOST_BITS + 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 > MOST_BITS+signum)) { significandOverflow = 1; TclBNInitBignumFromWideUInt(&significandBig, significandWide); } returnInteger: if (!significandOverflow) { if (significandWide > (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { #ifndef NO_WIDE_TYPE if (significandWide <= MOST_BITS+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; case INITIAL: /* This case only to silence compiler warning */ Tcl_Panic("TclParseNumber: state INITIAL can't happen here"); } } /* * Format an error message when an invalid number is encountered. */ if (status != TCL_OK) { if (interp != NULL) { Tcl_Obj *msg; TclNewLiteralStringObj(msg, "expected "); Tcl_AppendToObj(msg, expected, -1); Tcl_AppendToObj(msg, " but got \"", -1); Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, ""); Tcl_AppendToObj(msg, "\"", -1); if (state == BAD_OCTAL) { Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1); } Tcl_SetObjResult(interp, msg); } } /* * Free memory. */ if (octalSignificandOverflow) { mp_clear(&octalSignificandBig); } 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. * *---------------------------------------------------------------------- */ 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; Tcl_WideUInt w; /* * Try wide multiplication first */ if (!bignumFlag) { w = *wideRepPtr; if (w == 0) { /* * There's no need to multiply if the multiplicand is zero. */ *wideRepPtr = digit; return 0; } else if (numZeros >= maxpow10_wide || w > ((~(Tcl_WideUInt)0)-digit)/pow10_wide[numZeros+1]) { /* * Wide multiplication will overflow. Expand the * number to a bignum and fall through into the bignum case. */ TclBNInitBignumFromWideUInt (bignumRepPtr, w); } else { /* * Wide multiplication. */ *wideRepPtr = w * pow10_wide[numZeros+1] + digit; return 0; } } /* * Bignum multiplication. */ 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). */ 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); mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr); } return 1; }  /* *---------------------------------------------------------------------- * * 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 */ /* * With gcc on x86, the floating point rounding mode is double-extended. * This causes the result of double-precision calculations to be rounded * twice: once to the precision of double-extended and then again to the * precision of double. Double-rounding introduces gratuitous errors of 1 * ulp, so we need to change rounding mode to 53-bits. */ #if defined(__GNUC__) && defined(__i386) fpu_control_t roundTo53Bits = 0x027f; fpu_control_t oldRoundingMode; _FPU_GETCW(oldRoundingMode); _FPU_SETCW(roundTo53Bits); #endif /* * 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. */ retval = (double)(Tcl_WideInt)significand * pow10[ exponent ]; goto returnValue; } 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; } } } /* * 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); mp_clear(&significandBig); /* * Come here to return the computed value. */ returnValue: 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 */ /* * With gcc on x86, the floating point rounding mode is double-extended. * This causes the result of double-precision calculations to be rounded * twice: once to the precision of double-extended and then again to the * precision of double. Double-rounding introduces gratuitous errors of 1 * ulp, so we need to change rounding mode to 53-bits. */ #if defined(__GNUC__) && defined(__i386) fpu_control_t roundTo53Bits = 0x027f; fpu_control_t oldRoundingMode; _FPU_GETCW(oldRoundingMode); _FPU_SETCW(roundTo53Bits); #endif /* * Quick checks for over/underflow. */ if (numSigDigs+exponent-1 > maxDigits) { retval = HUGE_VAL; goto returnValue; } 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. */ 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; } /* * 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); /* * Come here to return the computed value. */ returnValue: 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; }  /* *---------------------------------------------------------------------- * * MakeNaN -- * * Makes a "Not a Number" given a set of bits to put in the tag bits * * 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; } if (n770_fp) { theNaN.iv = Nokia770Twiddle(theNaN.iv); } 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. * *---------------------------------------------------------------------- */ static double 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 msb; /* Most significant bit position of an * intermediate result */ int nDigits; /* Number of mp_digit's in an intermediate * result */ mp_int twoMv; /* Approx binary value expressed as an exact * integer scaled by the multiplier 2M */ mp_int twoMd; /* Exact decimal value expressed as an exact * integer scaled by the multiplier 2M */ int scale; /* Scale factor for M */ int multiplier; /* Power of two to scale M */ double num, den; /* Numerator and denominator of the correction * term */ double quot; /* Correction term */ double minincr; /* Lower bound on the absolute value of the * correction term. */ int i; /* * The first approximation is always low. If we find that it's HUGE_VAL, * we're done. */ if (approxResult == HUGE_VAL) { return approxResult; } /* * Find a common denominator for the decimal and binary fractions. The * common denominator will be 2**M2 + 5**M5. */ significand = frexp(approxResult, &binExponent); i = mantBits - binExponent; if (i < 0) { M2 = 0; } else { M2 = i; } if (exponent > 0) { M5 = 0; } else { M5 = -exponent; if ((M5-1) > M2) { M2 = M5-1; } } /* * The floating point number is significand*2**binExponent. 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 */ nDigits = msb / DIGIT_BIT + 1; mp_init_size(&twoMv, nDigits); i = (msb % DIGIT_BIT + 1); twoMv.used = nDigits; significand *= SafeLdExp(1.0, i); while (--nDigits >= 0) { twoMv.dp[nDigits] = (mp_digit) significand; significand -= (mp_digit) significand; significand = SafeLdExp(significand, DIGIT_BIT); } for (i = 0; i <= 8; ++i) { if (M5 & (1 << i)) { mp_mul(&twoMv, pow5+i, &twoMv); } } /* * Collect the decimal significand as a high precision integer. The least * significant bit corresponds to bit M2+exponent+1 so it will need to be * shifted left by that many bits after being multiplied by * 5**(M5+exponent). */ mp_init_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); /* * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction * term. Because 2M may well overflow a double, we need to scale the * denominator by a factor of 2**binExponent-mantBits */ scale = binExponent - mantBits - 1; mp_set(&twoMv, 1); for (i=0; i<=8; ++i) { if (M5 & (1 << i)) { mp_mul(&twoMv, pow5+i, &twoMv); } } multiplier = M2 + scale + 1; if (multiplier > 0) { mp_mul_2d(&twoMv, multiplier, &twoMv); } else if (multiplier < 0) { mp_div_2d(&twoMv, -multiplier, &twoMv, NULL); } /* * If the result is less than unity, the error is less than 1/2 unit in * the last place, so there's no correction to make. */ if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) { mp_clear(&twoMd); mp_clear(&twoMv); return approxResult; } /* * Convert the numerator and denominator of the corrector term accurately * to floating point numbers. */ num = TclBignumToDouble(&twoMd); den = TclBignumToDouble(&twoMv); quot = SafeLdExp(num/den, scale); minincr = SafeLdExp(1.0, binExponent-mantBits); if (quot<0. && quot>-minincr) { quot = -minincr; } else if (quot>0. && quot<minincr) { quot = minincr; } mp_clear(&twoMd); mp_clear(&twoMv); return approxResult + quot; }  /* *---------------------------------------------------------------------- * * TclDoubleDigits -- * * Converts a double to a string of digits. * * Results: * Returns the position of the character in the string after which the * decimal point should appear. Since the string contains only * significant digits, the position may be less than zero or greater than * the length of the string. * * Side effects: * Stores the digits in the given buffer and sets 'signum' according to * the sign of the number. * *---------------------------------------------------------------------- */ int TclDoubleDigits( char *buffer, /* 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. */ mp_int mplus; /* Scaled epsilon: (r + 2* mplus) == v(+) * where v(+) is the floating point successor * of v. */ mp_int mminus; /* Scaled epsilon: (r - 2*mminus) == v(-) * where v(-) is the floating point * predecessor of v. */ mp_int temp; int rfac2 = 0; /* Powers of 2 and 5 by which large */ int rfac5 = 0; /* integers should be scaled. */ int sfac2 = 0; int sfac5 = 0; int mplusfac2 = 0; int mminusfac2 = 0; char c; int i, k, n; /* * Split the number into absolute value and signum. */ v = AbsoluteValue(v, signum); /* * Handle zero specially. */ if (v == 0.0) { *buffer++ = '0'; *buffer++ = '\0'; return 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. */ smallestSig = GetIntegerTimesPower(v, &r, &e); lowOK = highOK = (mp_iseven(&r)); /* * We are going to want to develop integers r, s, mplus, and mminus such * that v = r / s, v(+)-v / 2 = mplus / s; v-v(-) / 2 = mminus / s and * then scale either s or r, mplus, mminus by an appropriate power of ten. * * We actually do this by keeping track of the powers of 2 and 5 by which * f is multiplied to yield v and by which 1 is multiplied to yield s, * mplus, and mminus. */ if (e >= 0) { int bits = e * log2FLT_RADIX; if (!smallestSig) { /* * Normal case, m+ and m- are both FLT_RADIX**e */ rfac2 = bits + 1; sfac2 = 1; mplusfac2 = bits; mminusfac2 = bits; } else { /* * If f is equal to the smallest significand, then we need another * factor of FLT_RADIX in s to cope with stepping to the next * smaller exponent when going to e's predecessor. */ rfac2 = bits + log2FLT_RADIX + 1; sfac2 = 1 + log2FLT_RADIX; mplusfac2 = bits + log2FLT_RADIX; mminusfac2 = bits; } } else { /* * v has digits after the binary point */ if (e <= DBL_MIN_EXP-DBL_MANT_DIG || !smallestSig) { /* * Either f isn't the smallest significand or e is the smallest * exponent. mplus and mminus will both be 1. */ rfac2 = 1; sfac2 = 1 - e * log2FLT_RADIX; mplusfac2 = 0; mminusfac2 = 0; } else { /* * f is the smallest significand, but e is not the smallest * exponent. We need to scale by FLT_RADIX again to cope with the * fact that v's predecessor has a smaller exponent. */ rfac2 = 1 + log2FLT_RADIX; sfac2 = 1 + log2FLT_RADIX * (1 - e); mplusfac2 = FLT_RADIX; mminusfac2 = 0; } } /* * Estimate the highest power of ten that will be needed to hold the * result. */ k = (int) ceil(log(v) / log(10.)); if (k >= 0) { sfac2 += k; sfac5 = k; } else { rfac2 -= k; mplusfac2 -= k; mminusfac2 -= k; rfac5 = -k; } /* * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5. */ mp_init_set(&mplus, 1); for (i=0 ; i<=8 ; ++i) { if (rfac5 & (1 << i)) { mp_mul(&mplus, pow5+i, &mplus); } } mp_mul(&r, &mplus, &r); mp_mul_2d(&r, rfac2, &r); mp_init_copy(&mminus, &mplus); mp_mul_2d(&mplus, mplusfac2, &mplus); mp_mul_2d(&mminus, mminusfac2, &mminus); mp_init_set(&s, 1); for (i=0 ; i<=8 ; ++i) { if (sfac5 & (1 << i)) { mp_mul(&s, pow5+i, &s); } } mp_mul_2d(&s, sfac2, &s); /* * It is possible for k to be off by one because we used an inexact * logarithm. */ mp_init(&temp); mp_add(&r, &mplus, &temp); i = mp_cmp_mag(&temp, &s); if (i>0 || (highOK && i==0)) { mp_mul_d(&s, 10, &s); ++k; } else { mp_mul_d(&temp, 10, &temp); i = mp_cmp_mag(&temp, &s); if (i<0 || (highOK && i==0)) { mp_mul_d(&r, 10, &r); mp_mul_d(&mplus, 10, &mplus); mp_mul_d(&mminus, 10, &mminus); --k; } } /* * At this point, k contains the power of ten by which we're scaling the * result. r/s is at least 1/10 and strictly less than ten, and v = r/s * * 10**k. mplus and mminus give the rounding limits. */ for (;;) { int tc1, tc2; mp_mul_d(&r, 10, &r); mp_div(&r, &s, &temp, &r); /* temp = 10r / s; r = 10r mod s */ i = temp.dp[0]; mp_mul_d(&mplus, 10, &mplus); mp_mul_d(&mminus, 10, &mminus); tc1 = mp_cmp_mag(&r, &mminus); if (lowOK) { tc1 = (tc1 <= 0); } else { tc1 = (tc1 < 0); } mp_add(&r, &mplus, &temp); tc2 = mp_cmp_mag(&temp, &s); if (highOK) { tc2 = (tc2 >= 0); } else { tc2= (tc2 > 0); } if (!tc1) { if (!tc2) { *buffer++ = '0' + i; } else { c = (char) (i + '1'); break; } } else { if (!tc2) { c = (char) (i + '0'); } else { mp_mul_2d(&r, 1, &r); n = mp_cmp_mag(&r, &s); if (n < 0) { c = (char) (i + '0'); } else { c = (char) (i + '1'); } } break; } }; *buffer++ = c; *buffer++ = '\0'; /* * Free memory, and return. */ mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL); return k; }  /* *---------------------------------------------------------------------- * * 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 (n770_fp) { bitwhack.iv = Nokia770Twiddle(bitwhack.iv); } if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { *signum = 1; bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); if (n770_fp) { bitwhack.iv = Nokia770Twiddle(bitwhack.iv); } 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, f; int e, i, 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 * 'double' * * Results: * None. * * Side effects: * The log base 2 of the floating point radix, the number of bits in a * double mantissa, and a table of the powers of five and ten are * computed and stored. * *---------------------------------------------------------------------- */ void TclInitDoubleConversion(void) { int i; int x; Tcl_WideUInt u; double d; #ifdef IEEE_FLOATING_POINT union { double dv; Tcl_WideUInt iv; } bitwhack; #endif /* * 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 *) ckalloc((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; /* * 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; } for (i=0 ; i<=mmaxpow ; ++i) { pow10[i] = d; d *= 10.0; } /* * 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); } /* * 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.)); /* * Nokia 770's software-emulated floating point is "middle endian": the * bytes within a 32-bit word are little-endian (like the native * integers), but the two words of a 'double' are presented most * significant word first. */ #ifdef IEEE_FLOATING_POINT bitwhack.dv = 1.000000238418579; /* 3ff0 0000 4000 0000 */ if ((bitwhack.iv >> 32) == 0x3ff00000) { n770_fp = 0; } else if ((bitwhack.iv & 0xffffffff) == 0x3ff00000) { n770_fp = 1; } else { Tcl_Panic("unknown floating point word order on this machine"); } #endif }  /* *---------------------------------------------------------------------- * * TclFinalizeDoubleConversion -- * * Cleans up this file on exit. * * Results: * None * * Side effects: * Memory allocated by TclInitDoubleConversion is freed. * *---------------------------------------------------------------------- */ void TclFinalizeDoubleConversion(void) { int i; Tcl_Free((char *) pow10_wide); for (i=0; i<9; ++i) { mp_clear(pow5 + i); } }  /* *---------------------------------------------------------------------- * * Tcl_InitBignumFromDouble -- * * 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 Tcl_InitBignumFromDouble( 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) { const 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; }  /* *---------------------------------------------------------------------- * * TclBignumToDouble -- * * Convert an arbitrary-precision integer to a native floating point * number. * * Results: * Returns the converted number. Sets errno to ERANGE if the number is * too large to convert. * *---------------------------------------------------------------------- */ double TclBignumToDouble( mp_int *a) /* Integer to convert. */ { mp_int b; int bits, shift, i; double r; /* * Determine how many bits we need, and extract that many from the input. * Round to nearest unit in the last place. */ bits = mp_count_bits(a); if (bits > DBL_MAX_EXP*log2FLT_RADIX) { errno = ERANGE; if (a->sign == MP_ZPOS) { return HUGE_VAL; } else { return -HUGE_VAL; } } shift = mantBits + 1 - bits; mp_init(&b); if (shift > 0) { mp_mul_2d(a, shift, &b); } else if (shift < 0) { mp_div_2d(a, -shift, &b, NULL); } else { mp_copy(a, &b); } mp_add_d(&b, 1, &b); mp_div_2d(&b, 1, &b, NULL); /* * Accumulate the result, one mp_digit at a time. */ r = 0.0; for (i=b.used-1 ; i>=0 ; --i) { r = ldexp(r, DIGIT_BIT) + b.dp[i]; } mp_clear(&b); /* * Scale the result to the correct number of bits. */ r = ldexp(r, bits - mantBits); /* * Return the result with the appropriate sign. */ if (a->sign == MP_ZPOS) { return r; } else { return -r; } } 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; return ((a->sign == MP_ZPOS) ? r : -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; }  /* *---------------------------------------------------------------------- * * SafeLdExp -- * * Do an 'ldexp' operation, but handle denormals gracefully. * * Results: * Returns the appropriately scaled value. * * On some platforms, 'ldexp' fails when presented with a number too * small to represent as a normalized double. This routine does 'ldexp' * in two steps for those numbers, to return correctly denormalized * values. * *---------------------------------------------------------------------- */ static double SafeLdExp( double fract, int expt) { int minexpt = DBL_MIN_EXP * log2FLT_RADIX; volatile double a, b, retval; if (expt < minexpt) { a = ldexp(fract, expt - mantBits - minexpt); b = ldexp(1.0, mantBits + minexpt); retval = a * b; } else { retval = ldexp(fract, expt); } return retval; }  /* *---------------------------------------------------------------------- * * TclFormatNaN -- * * Makes the string representation of a "Not a Number" * * Results: * None. * * Side effects: * Stores the string representation in the supplied buffer, which must be * at least TCL_DOUBLE_SPACE characters. * *---------------------------------------------------------------------- */ void TclFormatNaN( double value, /* The Not-a-Number to format. */ char *buffer) /* String representation. */ { #ifndef IEEE_FLOATING_POINT strcpy(buffer, "NaN"); return; #else union { double dv; Tcl_WideUInt iv; } bitwhack; bitwhack.dv = value; if (n770_fp) { bitwhack.iv = Nokia770Twiddle(bitwhack.iv); } if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63); *buffer++ = '-'; } *buffer++ = 'N'; *buffer++ = 'a'; *buffer++ = 'N'; bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1; if (bitwhack.iv != 0) { sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv); } else { *buffer = '\0'; } #endif /* IEEE_FLOATING_POINT */ }  /* *---------------------------------------------------------------------- * * Nokia770Twiddle -- * * Transpose the two words of a number for Nokia 770 floating * point handling. * *---------------------------------------------------------------------- */ static Tcl_WideUInt Nokia770Twiddle( Tcl_WideUInt w) /* Number to transpose */ { return (((w >> 32) & 0xffffffff) | (w << 32)); }  /* *---------------------------------------------------------------------- * * TclNokia770Doubles -- * * Transpose the two words of a number for Nokia 770 floating * point handling. * *---------------------------------------------------------------------- */ int TclNokia770Doubles(void) { return n770_fp; }  /* * Local Variables: * mode: c * c-basic-offset: 4 * fill-column: 78 * End: */

• ソフトウェア開発
• ホームページ制作
• Doxygen版
• STLport版