/* * 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. */ #include "tclInt.h" #include "tommath.h" #include /* * 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 /* * Rounding controls. (Thanks a lot, Intel!) */ #ifdef __i386 /* * 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__) 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 #define TCL_IEEE_DOUBLE_ROUNDING \ fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING; \ fpu_control_t oldRoundingMode; \ _FPU_GETCW(oldRoundingMode); \ _FPU_SETCW(roundTo53Bits) #define TCL_DEFAULT_DOUBLE_ROUNDING \ _FPU_SETCW(oldRoundingMode) /* * Sun ProC needs sunmath for rounding control on x86 like gcc above. */ #elif defined(__sun) #include #define TCL_IEEE_DOUBLE_ROUNDING \ ieee_flags("set","precision","double",NULL) #define TCL_DEFAULT_DOUBLE_ROUNDING \ ieee_flags("clear","precision",NULL,NULL) /* * Other platforms are assumed to always operate in full IEEE mode, so we make * the macros to go in and out of that mode do nothing. */ #else /* !__GNUC__ && !__sun */ #define TCL_IEEE_DOUBLE_ROUNDING ((void) 0) #define TCL_DEFAULT_DOUBLE_ROUNDING ((void) 0) #endif #else /* !__i386 */ #define TCL_IEEE_DOUBLE_ROUNDING ((void) 0) #define TCL_DEFAULT_DOUBLE_ROUNDING ((void) 0) #endif /* * MIPS floating-point units need special settings in control registers to use * gradual underflow as we expect. This fix is for the MIPSpro compiler. */ #if defined(__sgi) && defined(_COMPILER_VERSION) #include #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). */ /* Magic constants */ #define LOG10_2 0.3010299956639812 #define TWO_OVER_3LOG10 0.28952965460216784 #define LOG10_3HALVES_PLUS_FUDGE 0.1760912590558 /* * Definitions of the parts of an IEEE754-format floating point number. */ #define SIGN_BIT 0x80000000 /* Mask for the sign bit in the first word of * a double. */ #define EXP_MASK 0x7ff00000 /* Mask for the exponent field in the first * word of a double. */ #define EXP_SHIFT 20 /* Shift count to make the exponent an * integer. */ #define HIDDEN_BIT (((Tcl_WideUInt) 0x00100000) << 32) /* Hidden 1 bit for the significand. */ #define HI_ORDER_SIG_MASK 0x000fffff /* Mask for the high-order part of the * significand in the first word of a * double. */ #define SIG_MASK (((Tcl_WideUInt) HI_ORDER_SIG_MASK << 32) \ | 0xffffffff) /* Mask for the 52-bit significand. */ #define FP_PRECISION 53 /* Number of bits of significand plus the * hidden bit. */ #define EXPONENT_BIAS 0x3ff /* Bias of the exponent 0. */ /* * Derived quantities. */ #define TEN_PMAX 22 /* floor(FP_PRECISION*log(2)/log(5)) */ #define QUICK_MAX 14 /* floor((FP_PRECISION-1)*log(2)/log(10))-1 */ #define BLETCH 0x10 /* Highest power of two that is greater than * DBL_MAX_10_EXP, divided by 16. */ #define DIGIT_GROUP 8 /* floor(DIGIT_BIT*log(2)/log(10)) */ /* * Union used to dismantle floating point numbers. */ typedef union Double { struct { #ifdef WORDS_BIGENDIAN int word0; int word1; #else int word1; int word0; #endif } w; double d; Tcl_WideUInt q; } Double; 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 pow10vals[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 = 0.0; /* 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 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. */ /* * Table of powers of 5 that are small enough to fit in an mp_digit. */ static const mp_digit dpow5[13] = { 1, 5, 25, 125, 625, 3125, 15625, 78125, 390625, 1953125, 9765625, 48828125, 244140625 }; /* * Table of powers: pow5_13[n] = 5**(13*2**(n+1)) */ static mp_int pow5_13[5]; /* Table of powers: 5**13, 5**26, 5**52, * 5**104, 5**208 */ static const double tens[] = { 1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 }; static const int itens [] = { 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 }; static const double bigtens[] = { 1e016, 1e032, 1e064, 1e128, 1e256 }; #define N_BIGTENS 5 static const int log2pow5[27] = { 01, 3, 5, 7, 10, 12, 14, 17, 19, 21, 24, 26, 28, 31, 33, 35, 38, 40, 42, 45, 47, 49, 52, 54, 56, 59, 61 }; #define N_LOG2POW5 27 static const Tcl_WideUInt wuipow5[27] = { (Tcl_WideUInt) 1, /* 5**0 */ (Tcl_WideUInt) 5, (Tcl_WideUInt) 25, (Tcl_WideUInt) 125, (Tcl_WideUInt) 625, (Tcl_WideUInt) 3125, /* 5**5 */ (Tcl_WideUInt) 3125*5, (Tcl_WideUInt) 3125*25, (Tcl_WideUInt) 3125*125, (Tcl_WideUInt) 3125*625, (Tcl_WideUInt) 3125*3125, /* 5**10 */ (Tcl_WideUInt) 3125*3125*5, (Tcl_WideUInt) 3125*3125*25, (Tcl_WideUInt) 3125*3125*125, (Tcl_WideUInt) 3125*3125*625, (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */ (Tcl_WideUInt) 3125*3125*3125*5, (Tcl_WideUInt) 3125*3125*3125*25, (Tcl_WideUInt) 3125*3125*3125*125, (Tcl_WideUInt) 3125*3125*3125*625, (Tcl_WideUInt) 3125*3125*3125*3125, /* 5**20 */ (Tcl_WideUInt) 3125*3125*3125*3125*5, (Tcl_WideUInt) 3125*3125*3125*3125*25, (Tcl_WideUInt) 3125*3125*3125*3125*125, (Tcl_WideUInt) 3125*3125*3125*3125*625, (Tcl_WideUInt) 3125*3125*3125*3125*3125, /* 5**25 */ (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */ }; /* * Static functions defined in this file. */ static int AccumulateDecimalDigit(unsigned, int, Tcl_WideUInt *, mp_int *, int); static double MakeHighPrecisionDouble(int signum, mp_int *significand, int nSigDigs, int exponent); static double MakeLowPrecisionDouble(int signum, Tcl_WideUInt significand, int nSigDigs, int exponent); #ifdef IEEE_FLOATING_POINT static double MakeNaN(int signum, Tcl_WideUInt tag); #endif static double RefineApproximation(double approx, mp_int *exactSignificand, int exponent); static void MulPow5(mp_int *, unsigned, mp_int *); static int NormalizeRightward(Tcl_WideUInt *); static int RequiredPrecision(Tcl_WideUInt); static void DoubleToExpAndSig(double, Tcl_WideUInt *, int *, int *); static void TakeAbsoluteValue(Double *, int *); static char * FormatInfAndNaN(Double *, int *, char **); static char * FormatZero(int *, char **); static int ApproximateLog10(Tcl_WideUInt, int, int); static int BetterLog10(double, int, int *); static void ComputeScale(int, int, int *, int *, int *, int *); static void SetPrecisionLimits(int, int, int *, int *, int *, int *); static char * BumpUp(char *, char *, int *); static int AdjustRange(double *, int); static char * ShorteningQuickFormat(double, int, int, double, char *, int *); static char * StrictQuickFormat(double, int, int, double, char *, int *); static char * QuickConversion(double, int, int, int, int, int, int, int *, char **); static void CastOutPowersOf2(int *, int *, int *); static char * ShorteningInt64Conversion(Double *, int, Tcl_WideUInt, int, int, int, int, int, int, int, int, int, int, int, int *, char **); static char * StrictInt64Conversion(Double *, int, Tcl_WideUInt, int, int, int, int, int, int, int, int, int *, char **); static int ShouldBankerRoundUpPowD(mp_int *, int, int); static int ShouldBankerRoundUpToNextPowD(mp_int *, mp_int *, int, int, int, mp_int *); static char * ShorteningBignumConversionPowD(Double *dPtr, int convType, Tcl_WideUInt bw, int b2, int b5, int m2plus, int m2minus, int m5, int sd, int k, int len, int ilim, int ilim1, int *decpt, char **endPtr); static char * StrictBignumConversionPowD(Double *dPtr, int convType, Tcl_WideUInt bw, int b2, int b5, int sd, int k, int len, int ilim, int ilim1, int *decpt, char **endPtr); static int ShouldBankerRoundUp(mp_int *, mp_int *, int); static int ShouldBankerRoundUpToNext(mp_int *, mp_int *, mp_int *, int, int, mp_int *); static char * ShorteningBignumConversion(Double *dPtr, int convType, Tcl_WideUInt bw, int b2, int m2plus, int m2minus, int s2, int s5, int k, int len, int ilim, int ilim1, int *decpt, char **endPtr); static char * StrictBignumConversion(Double *dPtr, int convType, Tcl_WideUInt bw, int b2, int s2, int s5, int k, int len, int ilim, int ilim1, int *decpt, char **endPtr); static double BignumToBiasedFrExp(const mp_int *big, int *machexp); static double Pow10TimesFrExp(int exponent, double fraction, int *machexp); static double SafeLdExp(double fraction, int exponent); #ifdef IEEE_FLOATING_POINT static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w); #endif /* *---------------------------------------------------------------------- * * 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, 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 */ #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 (TclIsSpaceProc(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_BINARY_ONLY) { goto zerob; } 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, b, 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' and 'b'. */ 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 (flags & TCL_PARSE_BINARY_ONLY) { goto zerob; } if (c == 'o' || c == 'O') { state = ZERO_O; break; } goto decimal; 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; } 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: zerob: 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. */ decimal: 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 (TclIsSpaceProc(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'; } else { goto endgame; } numSigDigs++; 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 && TclIsSpaceProc(*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 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: #ifdef IEEE_FLOATING_POINT case sN: case sNA: case sNANPAREN: case sNANHEX: Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'", acceptState, bytes); #endif 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 TCL_WIDE_INT_IS_LONG 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 TCL_WIDE_INT_IS_LONG 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; #ifdef IEEE_FLOATING_POINT case sNAN: case sNANFINISH: objPtr->internalRep.doubleValue = MakeNaN(signum,significandWide); objPtr->typePtr = &tclDoubleType; break; #endif 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 = Tcl_ObjPrintf("expected %s but got \"", expected); Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, ""); Tcl_AppendToObj(msg, "\"", -1); Tcl_SetObjResult(interp, msg); Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", NULL); } } /* * 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. */ TCL_IEEE_DOUBLE_ROUNDING; /* * 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 * pow10vals[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 * pow10vals[diff]); retval = factor * pow10vals[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 / pow10vals[-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. */ TCL_DEFAULT_DOUBLE_ROUNDING; 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. */ TCL_IEEE_DOUBLE_ROUNDING; /* * 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 (tiny == 0.0) { tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits); } 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. */ TCL_DEFAULT_DOUBLE_ROUNDING; 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>= 1; ++r; } if (p != result) { mp_copy(p, result); } } /* *---------------------------------------------------------------------- * * NormalizeRightward -- * * Shifts a number rightward until it is odd (that is, until the least * significant bit is nonzero. * * Results: * Returns the number of bit positions by which the number was shifted. * * Side effects: * Shifts the number in place; *wPtr is replaced by the shifted number. * *---------------------------------------------------------------------- */ inline static int NormalizeRightward( Tcl_WideUInt *wPtr) /* INOUT: Number to shift. */ { int rv = 0; Tcl_WideUInt w = *wPtr; if (!(w & (Tcl_WideUInt) 0xffffffff)) { w >>= 32; rv += 32; } if (!(w & (Tcl_WideUInt) 0xffff)) { w >>= 16; rv += 16; } if (!(w & (Tcl_WideUInt) 0xff)) { w >>= 8; rv += 8; } if (!(w & (Tcl_WideUInt) 0xf)) { w >>= 4; rv += 4; } if (!(w & 0x3)) { w >>= 2; rv += 2; } if (!(w & 0x1)) { w >>= 1; ++rv; } *wPtr = w; return rv; } /* *---------------------------------------------------------------------- * * RequiredPrecision -- * * Determines the number of bits needed to hold an intger. * * Results: * Returns the position of the most significant bit (0 - 63). Returns 0 * if the number is zero. * *---------------------------------------------------------------------- */ static int RequiredPrecision( Tcl_WideUInt w) /* Number to interrogate. */ { int rv; unsigned long wi; if (w & ((Tcl_WideUInt) 0xffffffff << 32)) { wi = (unsigned long) (w >> 32); rv = 32; } else { wi = (unsigned long) w; rv = 0; } if (wi & 0xffff0000) { wi >>= 16; rv += 16; } if (wi & 0xff00) { wi >>= 8; rv += 8; } if (wi & 0xf0) { wi >>= 4; rv += 4; } if (wi & 0xc) { wi >>= 2; rv += 2; } if (wi & 0x2) { wi >>= 1; ++rv; } if (wi & 0x1) { ++rv; } return rv; } /* *---------------------------------------------------------------------- * * DoubleToExpAndSig -- * * Separates a 'double' into exponent and significand. * * Side effects: * Stores the significand in '*significand' and the exponent in '*expon' * so that dv == significand * 2.0**expon, and significand is odd. Also * stores the position of the leftmost 1-bit in 'significand' in 'bits'. * *---------------------------------------------------------------------- */ inline static void DoubleToExpAndSig( double dv, /* Number to convert. */ Tcl_WideUInt *significand, /* OUTPUT: Significand of the number. */ int *expon, /* OUTPUT: Exponent to multiply the number * by. */ int *bits) /* OUTPUT: Number of significant bits. */ { Double d; /* Number being converted. */ Tcl_WideUInt z; /* Significand under construction. */ int de; /* Exponent of the number. */ int k; /* Bit count. */ d.d = dv; /* * Extract exponent and significand. */ de = (d.w.word0 & EXP_MASK) >> EXP_SHIFT; z = d.q & SIG_MASK; if (de != 0) { z |= HIDDEN_BIT; k = NormalizeRightward(&z); *bits = FP_PRECISION - k; *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1); } else { k = NormalizeRightward(&z); *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1) + 1; *bits = RequiredPrecision(z); } *significand = z; } /* *---------------------------------------------------------------------- * * TakeAbsoluteValue -- * * Takes the absolute value of a 'double' including 0, Inf and NaN * * Side effects: * The 'double' in *d is replaced with its absolute value. The signum is * stored in 'sign': 1 for negative, 0 for nonnegative. * *---------------------------------------------------------------------- */ inline static void TakeAbsoluteValue( Double *d, /* Number to replace with absolute value. */ int *sign) /* Place to put the signum. */ { if (d->w.word0 & SIGN_BIT) { *sign = 1; d->w.word0 &= ~SIGN_BIT; } else { *sign = 0; } } /* *---------------------------------------------------------------------- * * FormatInfAndNaN -- * * Bailout for formatting infinities and Not-A-Number. * * Results: * Returns one of the strings 'Infinity' and 'NaN'. The string returned * must be freed by the caller using 'ckfree'. * * Side effects: * Stores 9999 in *decpt, and sets '*endPtr' to designate the terminating * NUL byte of the string if 'endPtr' is not NULL. * *---------------------------------------------------------------------- */ inline static char * FormatInfAndNaN( Double *d, /* Exceptional number to format. */ int *decpt, /* Decimal point to set to a bogus value. */ char **endPtr) /* Pointer to the end of the formatted data */ { char *retval; *decpt = 9999; if (!(d->w.word1) && !(d->w.word0 & HI_ORDER_SIG_MASK)) { retval = ckalloc(9); strcpy(retval, "Infinity"); if (endPtr) { *endPtr = retval + 8; } } else { retval = ckalloc(4); strcpy(retval, "NaN"); if (endPtr) { *endPtr = retval + 3; } } return retval; } /* *---------------------------------------------------------------------- * * FormatZero -- * * Bailout to format a zero floating-point number. * * Results: * Returns the constant string "0" * * Side effects: * Stores 1 in '*decpt' and puts a pointer to the NUL byte terminating * the string in '*endPtr' if 'endPtr' is not NULL. * *---------------------------------------------------------------------- */ inline static char * FormatZero( int *decpt, /* Location of the decimal point. */ char **endPtr) /* Pointer to the end of the formatted data */ { char *retval = ckalloc(2); strcpy(retval, "0"); if (endPtr) { *endPtr = retval+1; } *decpt = 0; return retval; } /* *---------------------------------------------------------------------- * * ApproximateLog10 -- * * Computes a two-term Taylor series approximation to the common log of a * number, and computes the number's binary log. * * Results: * Return an approximation to floor(log10(bw*2**be)) that is either exact * or 1 too high. * *---------------------------------------------------------------------- */ inline static int ApproximateLog10( Tcl_WideUInt bw, /* Integer significand of the number. */ int be, /* Power of two to scale bw. */ int bbits) /* Number of bits of precision in bw. */ { int i; /* Log base 2 of the number. */ int k; /* Floor(Log base 10 of the number) */ double ds; /* Mantissa of the number. */ Double d2; /* * Compute i and d2 such that d = d2*2**i, and 1 < d2 < 2. * Compute an approximation to log10(d), * log10(d) ~ log10(2) * i + log10(1.5) * + (significand-1.5)/(1.5 * log(10)) */ d2.q = bw << (FP_PRECISION - bbits) & SIG_MASK; d2.w.word0 |= (EXPONENT_BIAS) << EXP_SHIFT; i = be + bbits - 1; ds = (d2.d - 1.5) * TWO_OVER_3LOG10 + LOG10_3HALVES_PLUS_FUDGE + LOG10_2 * i; k = (int) ds; if (k > ds) { --k; } return k; } /* *---------------------------------------------------------------------- * * BetterLog10 -- * * Improves the result of ApproximateLog10 for numbers in the range * 1 .. 10**(TEN_PMAX)-1 * * Side effects: * Sets k_check to 0 if the new result is known to be exact, and to 1 if * it may still be one too high. * * Results: * Returns the improved approximation to log10(d). * *---------------------------------------------------------------------- */ inline static int BetterLog10( double d, /* Original number to format. */ int k, /* Characteristic(Log base 10) of the * number. */ int *k_check) /* Flag == 1 if k is inexact. */ { /* * Performance hack. If k is in the range 0..TEN_PMAX, then we can use a * powers-of-ten table to check it. */ if (k >= 0 && k <= TEN_PMAX) { if (d < tens[k]) { k--; } *k_check = 0; } else { *k_check = 1; } return k; } /* *---------------------------------------------------------------------- * * ComputeScale -- * * Prepares to format a floating-point number as decimal. * * Parameters: * floor(log10*x) is k (or possibly k-1). floor(log2(x) is i. The * significand of x requires bbits bits to represent. * * Results: * Determines integers b2, b5, s2, s5 so that sig*2**b2*5**b5/2**s2*2**s5 * exactly represents the value of the x/10**k. This value will lie in * the range [1 .. 10), and allows for computing successive digits by * multiplying sig%10 by 10. * *---------------------------------------------------------------------- */ inline static void ComputeScale( int be, /* Exponent part of number: d = bw * 2**be. */ int k, /* Characteristic of log10(number). */ int *b2, /* OUTPUT: Power of 2 in the numerator. */ int *b5, /* OUTPUT: Power of 5 in the numerator. */ int *s2, /* OUTPUT: Power of 2 in the denominator. */ int *s5) /* OUTPUT: Power of 5 in the denominator. */ { /* * Scale numerator and denominator powers of 2 so that the input binary * number is the ratio of integers. */ if (be <= 0) { *b2 = 0; *s2 = -be; } else { *b2 = be; *s2 = 0; } /* * Scale numerator and denominator so that the output decimal number is * the ratio of integers. */ if (k >= 0) { *b5 = 0; *s5 = k; *s2 += k; } else { *b2 -= k; *b5 = -k; *s5 = 0; } } /* *---------------------------------------------------------------------- * * SetPrecisionLimits -- * * Determines how many digits of significance should be computed (and, * hence, how much memory need be allocated) for formatting a floating * point number. * * Given that 'k' is floor(log10(x)): * if 'shortest' format is used, there will be at most 18 digits in the * result. * if 'F' format is used, there will be at most 'ndigits' + k + 1 digits * if 'E' format is used, there will be exactly 'ndigits' digits. * * Side effects: * Adjusts '*ndigitsPtr' to have a valid value. Stores the maximum memory * allocation needed in *iPtr. Sets '*iLimPtr' to the limiting number of * digits to convert if k has been guessed correctly, and '*iLim1Ptr' to * the limiting number of digits to convert if k has been guessed to be * one too high. * *---------------------------------------------------------------------- */ inline static void SetPrecisionLimits( int convType, /* Type of conversion: TCL_DD_SHORTEST, * TCL_DD_STEELE0, TCL_DD_E_FMT, * TCL_DD_F_FMT. */ int k, /* Floor(log10(number to convert)) */ int *ndigitsPtr, /* IN/OUT: Number of digits requested (will be * adjusted if needed). */ int *iPtr, /* OUT: Maximum number of digits to return. */ int *iLimPtr, /* OUT: Number of digits of significance if * the bignum method is used.*/ int *iLim1Ptr) /* OUT: Number of digits of significance if * the quick method is used. */ { switch (convType) { case TCL_DD_SHORTEST0: case TCL_DD_STEELE0: *iLimPtr = *iLim1Ptr = -1; *iPtr = 18; *ndigitsPtr = 0; break; case TCL_DD_E_FORMAT: if (*ndigitsPtr <= 0) { *ndigitsPtr = 1; } *iLimPtr = *iLim1Ptr = *iPtr = *ndigitsPtr; break; case TCL_DD_F_FORMAT: *iPtr = *ndigitsPtr + k + 1; *iLimPtr = *iPtr; *iLim1Ptr = *iPtr - 1; if (*iPtr <= 0) { *iPtr = 1; } break; default: *iPtr = -1; *iLimPtr = -1; *iLim1Ptr = -1; Tcl_Panic("impossible conversion type in TclDoubleDigits"); } } /* *---------------------------------------------------------------------- * * BumpUp -- * * Increases a string of digits ending in a series of nines to designate * the next higher number. xxxxb9999... -> xxxx(b+1)0000... * * Results: * Returns a pointer to the end of the adjusted string. * * Side effects: * In the case that the string consists solely of '999999', sets it to * "1" and moves the decimal point (*kPtr) one place to the right. * *---------------------------------------------------------------------- */ inline static char * BumpUp( char *s, /* Cursor pointing one past the end of the * string. */ char *retval, /* Start of the string of digits. */ int *kPtr) /* Position of the decimal point. */ { while (*--s == '9') { if (s == retval) { ++(*kPtr); *s = '1'; return s+1; } } ++*s; ++s; return s; } /* *---------------------------------------------------------------------- * * AdjustRange -- * * Rescales a 'double' in preparation for formatting it using the 'quick' * double-to-string method. * * Results: * Returns the precision that has been lost in the prescaling as a count * of units in the least significant place. * *---------------------------------------------------------------------- */ inline static int AdjustRange( double *dPtr, /* INOUT: Number to adjust. */ int k) /* IN: floor(log10(d)) */ { int ieps; /* Number of roundoff errors that have * accumulated. */ double d = *dPtr; /* Number to adjust. */ double ds; int i, j, j1; ieps = 2; if (k > 0) { /* * The number must be reduced to bring it into range. */ ds = tens[k & 0xf]; j = k >> 4; if (j & BLETCH) { j &= (BLETCH-1); d /= bigtens[N_BIGTENS - 1]; ieps++; } i = 0; for (; j != 0; j>>=1) { if (j & 1) { ds *= bigtens[i]; ++ieps; } ++i; } d /= ds; } else if ((j1 = -k) != 0) { /* * The number must be increased to bring it into range. */ d *= tens[j1 & 0xf]; i = 0; for (j = j1>>4; j; j>>=1) { if (j & 1) { ieps++; d *= bigtens[i]; } ++i; } } *dPtr = d; return ieps; } /* *---------------------------------------------------------------------- * * ShorteningQuickFormat -- * * Returns a 'quick' format of a double precision number to a string of * digits, preferring a shorter string of digits if the shorter string is * still within 1/2 ulp of the number. * * Results: * Returns the string of digits. Returns NULL if the 'quick' method fails * and the bignum method must be used. * * Side effects: * Stores the position of the decimal point at '*kPtr'. * *---------------------------------------------------------------------- */ inline static char * ShorteningQuickFormat( double d, /* Number to convert. */ int k, /* floor(log10(d)) */ int ilim, /* Number of significant digits to return. */ double eps, /* Estimated roundoff error. */ char *retval, /* Buffer to receive the digit string. */ int *kPtr) /* Pointer to stash the position of the * decimal point. */ { char *s = retval; /* Cursor in the return value. */ int digit; /* Current digit. */ int i; eps = 0.5 / tens[ilim-1] - eps; i = 0; for (;;) { /* * Convert a digit. */ digit = (int) d; d -= digit; *s++ = '0' + digit; /* * Truncate the conversion if the string of digits is within 1/2 ulp * of the actual value. */ if (d < eps) { *kPtr = k; return s; } if ((1. - d) < eps) { *kPtr = k; return BumpUp(s, retval, kPtr); } /* * Bail out if the conversion fails to converge to a sufficiently * precise value. */ if (++i >= ilim) { return NULL; } /* * Bring the next digit to the integer part. */ eps *= 10; d *= 10.0; } } /* *---------------------------------------------------------------------- * * StrictQuickFormat -- * * Convert a double precision number of a string of a precise number of * digits, using the 'quick' double precision method. * * Results: * Returns the digit string, or NULL if the bignum method must be used to * do the formatting. * * Side effects: * Stores the position of the decimal point in '*kPtr'. * *---------------------------------------------------------------------- */ inline static char * StrictQuickFormat( double d, /* Number to convert. */ int k, /* floor(log10(d)) */ int ilim, /* Number of significant digits to return. */ double eps, /* Estimated roundoff error. */ char *retval, /* Start of the digit string. */ int *kPtr) /* Pointer to stash the position of the * decimal point. */ { char *s = retval; /* Cursor in the return value. */ int digit; /* Current digit of the answer. */ int i; eps *= tens[ilim-1]; i = 1; for (;;) { /* * Extract a digit. */ digit = (int) d; d -= digit; if (d == 0.0) { ilim = i; } *s++ = '0' + digit; /* * When the given digit count is reached, handle trailing strings of 0 * and 9. */ if (i == ilim) { if (d > 0.5 + eps) { *kPtr = k; return BumpUp(s, retval, kPtr); } else if (d < 0.5 - eps) { while (*--s == '0') { /* do nothing */ } s++; *kPtr = k; return s; } else { return NULL; } } /* * Advance to the next digit. */ ++i; d *= 10.0; } } /* *---------------------------------------------------------------------- * * QuickConversion -- * * Converts a floating point number the 'quick' way, when only a limited * number of digits is required and floating point arithmetic can * therefore be used for the intermediate results. * * Results: * Returns the converted string, or NULL if the bignum method must be * used. * *---------------------------------------------------------------------- */ inline static char * QuickConversion( double e, /* Number to format. */ int k, /* floor(log10(d)), approximately. */ int k_check, /* 0 if k is exact, 1 if it may be too high */ int flags, /* Flags passed to dtoa: * TCL_DD_SHORTEN_FLAG */ int len, /* Length of the return value. */ int ilim, /* Number of digits to store. */ int ilim1, /* Number of digits to store if we misguessed * k. */ int *decpt, /* OUTPUT: Location of the decimal point. */ char **endPtr) /* OUTPUT: Pointer to the terminal null * byte. */ { int ieps; /* Number of 1-ulp roundoff errors that have * accumulated in the calculation. */ Double eps; /* Estimated roundoff error. */ char *retval; /* Returned string. */ char *end; /* Pointer to the terminal null byte in the * returned string. */ volatile double d; /* Workaround for a bug in mingw gcc 3.4.5 */ /* * Bring d into the range [1 .. 10). */ ieps = AdjustRange(&e, k); d = e; /* * If the guessed value of k didn't get d into range, adjust it by one. If * that leaves us outside the range in which quick format is accurate, * bail out. */ if (k_check && d < 1. && ilim > 0) { if (ilim1 < 0) { return NULL; } ilim = ilim1; --k; d *= 10.0; ++ieps; } /* * Compute estimated roundoff error. */ eps.d = ieps * d + 7.; eps.w.word0 -= (FP_PRECISION-1) << EXP_SHIFT; /* * Handle the peculiar case where the result has no significant digits. */ retval = ckalloc(len + 1); if (ilim == 0) { d -= 5.; if (d > eps.d) { *retval = '1'; *decpt = k; return retval; } else if (d < -eps.d) { *decpt = k; return retval; } else { ckfree(retval); return NULL; } } /* * Format the digit string. */ if (flags & TCL_DD_SHORTEN_FLAG) { end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt); } else { end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt); } if (end == NULL) { ckfree(retval); return NULL; } *end = '\0'; if (endPtr != NULL) { *endPtr = end; } return retval; } /* *---------------------------------------------------------------------- * * CastOutPowersOf2 -- * * Adjust the factors 'b2', 'm2', and 's2' to cast out common powers of 2 * from numerator and denominator in preparation for the 'bignum' method * of floating point conversion. * *---------------------------------------------------------------------- */ inline static void CastOutPowersOf2( int *b2, /* Power of 2 to multiply the significand. */ int *m2, /* Power of 2 to multiply 1/2 ulp. */ int *s2) /* Power of 2 to multiply the common * denominator. */ { int i; if (*m2 > 0 && *s2 > 0) { /* Find the smallest power of 2 in the * numerator. */ if (*m2 < *s2) { /* Find the lowest common denominator. */ i = *m2; } else { i = *s2; } *b2 -= i; /* Reduce to lowest terms. */ *m2 -= i; *s2 -= i; } } /* *---------------------------------------------------------------------- * * ShorteningInt64Conversion -- * * Converts a double-precision number to the shortest string of digits * that reconverts exactly to the given number, or to 'ilim' digits if * that will yield a shorter result. The numerator and denominator in * David Gay's conversion algorithm are known to fit in Tcl_WideUInt, * giving considerably faster arithmetic than mp_int's. * * Results: * Returns the string of significant decimal digits, in newly allocated * memory * * Side effects: * Stores the location of the decimal point in '*decpt' and the location * of the terminal null byte in '*endPtr'. * *---------------------------------------------------------------------- */ inline static char * ShorteningInt64Conversion( Double *dPtr, /* Original number to convert. */ int convType, /* Type of conversion (shortest, Steele, * E format, F format). */ Tcl_WideUInt bw, /* Integer significand. */ int b2, int b5, /* Scale factor for the significand in the * numerator. */ int m2plus, int m2minus, int m5, /* Scale factors for 1/2 ulp in the numerator * (will be different if bw == 1. */ int s2, int s5, /* Scale factors for the denominator. */ int k, /* Number of output digits before the decimal * point. */ int len, /* Number of digits to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Position of the terminal '\0' at * the end of the returned string. */ { char *retval = ckalloc(len + 1); /* Output buffer. */ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2; /* Numerator of the fraction being * converted. */ Tcl_WideUInt S = wuipow5[s5] << s2; /* Denominator of the fraction being * converted. */ Tcl_WideUInt mplus, mminus; /* Ranges for testing whether the result is * within roundoff of being exact. */ int digit; /* Current output digit. */ char *s = retval; /* Cursor in the output buffer. */ int i; /* Current position in the output buffer. */ /* * Adjust if the logarithm was guessed wrong. */ if (b < S) { b = 10 * b; ++m2plus; ++m2minus; ++m5; ilim = ilim1; --k; } /* * Compute roundoff ranges. */ mplus = wuipow5[m5] << m2plus; mminus = wuipow5[m5] << m2minus; /* * Loop through the digits. */ i = 1; for (;;) { digit = (int)(b / S); if (digit > 10) { Tcl_Panic("wrong digit!"); } b = b % S; /* * Does the current digit put us on the low side of the exact value * but within within roundoff of being exact? */ if (b < mplus || (b == mplus && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { /* * Make sure we shouldn't be rounding *up* instead, in case the * next number above is closer. */ if (2 * b > S || (2 * b == S && (digit & 1) != 0)) { ++digit; if (digit == 10) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } } /* * Stash the current digit. */ *s++ = '0' + digit; break; } /* * Does one plus the current digit put us within roundoff of the * number? */ if (b > S - mminus || (b == S - mminus && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { if (digit == 9) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } ++digit; *s++ = '0' + digit; break; } /* * Have we converted all the requested digits? */ *s++ = '0' + digit; if (i == ilim) { if (2*b > S || (2*b == S && (digit & 1) != 0)) { s = BumpUp(s, retval, &k); } break; } /* * Advance to the next digit. */ b = 10 * b; mplus = 10 * mplus; mminus = 10 * mminus; ++i; } /* * Endgame - store the location of the decimal point and the end of the * string. */ *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * StrictInt64Conversion -- * * Converts a double-precision number to a fixed-length string of 'ilim' * digits that reconverts exactly to the given number. ('ilim' should be * replaced with 'ilim1' in the case where log10(d) has been * overestimated). The numerator and denominator in David Gay's * conversion algorithm are known to fit in Tcl_WideUInt, giving * considerably faster arithmetic than mp_int's. * * Results: * Returns the string of significant decimal digits, in newly allocated * memory * * Side effects: * Stores the location of the decimal point in '*decpt' and the location * of the terminal null byte in '*endPtr'. * *---------------------------------------------------------------------- */ inline static char * StrictInt64Conversion( Double *dPtr, /* Original number to convert. */ int convType, /* Type of conversion (shortest, Steele, * E format, F format). */ Tcl_WideUInt bw, /* Integer significand. */ int b2, int b5, /* Scale factor for the significand in the * numerator. */ int s2, int s5, /* Scale factors for the denominator. */ int k, /* Number of output digits before the decimal * point. */ int len, /* Number of digits to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Position of the terminal '\0' at * the end of the returned string. */ { char *retval = ckalloc(len + 1); /* Output buffer. */ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2; /* Numerator of the fraction being * converted. */ Tcl_WideUInt S = wuipow5[s5] << s2; /* Denominator of the fraction being * converted. */ int digit; /* Current output digit. */ char *s = retval; /* Cursor in the output buffer. */ int i; /* Current position in the output buffer. */ /* * Adjust if the logarithm was guessed wrong. */ if (b < S) { b = 10 * b; ilim = ilim1; --k; } /* * Loop through the digits. */ i = 1; for (;;) { digit = (int)(b / S); if (digit > 10) { Tcl_Panic("wrong digit!"); } b = b % S; /* * Have we converted all the requested digits? */ *s++ = '0' + digit; if (i == ilim) { if (2*b > S || (2*b == S && (digit & 1) != 0)) { s = BumpUp(s, retval, &k); } else { while (*--s == '0') { /* do nothing */ } ++s; } break; } /* * Advance to the next digit. */ b = 10 * b; ++i; } /* * Endgame - store the location of the decimal point and the end of the * string. */ *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * ShouldBankerRoundUpPowD -- * * Test whether bankers' rounding should round a digit up. Assumption is * made that the denominator of the fraction being tested is a power of * 2**DIGIT_BIT. * * Results: * Returns 1 iff the fraction is more than 1/2, or if the fraction is * exactly 1/2 and the digit is odd. * *---------------------------------------------------------------------- */ inline static int ShouldBankerRoundUpPowD( mp_int *b, /* Numerator of the fraction. */ int sd, /* Denominator is 2**(sd*DIGIT_BIT). */ int isodd) /* 1 if the digit is odd, 0 if even. */ { int i; static const mp_digit topbit = 1 << (DIGIT_BIT - 1); if (b->used < sd || (b->dp[sd-1] & topbit) == 0) { return 0; } if (b->dp[sd-1] != topbit) { return 1; } for (i = sd-2; i >= 0; --i) { if (b->dp[i] != 0) { return 1; } } return isodd; } /* *---------------------------------------------------------------------- * * ShouldBankerRoundUpToNextPowD -- * * Tests whether bankers' rounding will round down in the "denominator is * a power of 2**MP_DIGIT" case. * * Results: * Returns 1 if the rounding will be performed - which increases the * digit by one - and 0 otherwise. * *---------------------------------------------------------------------- */ inline static int ShouldBankerRoundUpToNextPowD( mp_int *b, /* Numerator of the fraction. */ mp_int *m, /* Numerator of the rounding tolerance. */ int sd, /* Common denominator is 2**(sd*DIGIT_BIT). */ int convType, /* Conversion type: STEELE defeats * round-to-even (not sure why one wants to do * this; I copied it from Gay). FIXME */ int isodd, /* 1 if the integer significand is odd. */ mp_int *temp) /* Work area for the calculation. */ { int i; /* * Compare B and S-m - which is the same as comparing B+m and S - which we * do by computing b+m and doing a bitwhack compare against * 2**(DIGIT_BIT*sd) */ mp_add(b, m, temp); if (temp->used <= sd) { /* Too few digits to be > s */ return 0; } if (temp->used > sd+1 || temp->dp[sd] > 1) { /* >= 2s */ return 1; } for (i = sd-1; i >= 0; --i) { /* Check for ==s */ if (temp->dp[i] != 0) { /* > s */ return 1; } } if (convType == TCL_DD_STEELE0) { /* Biased rounding. */ return 0; } return isodd; } /* *---------------------------------------------------------------------- * * ShorteningBignumConversionPowD -- * * Converts a double-precision number to the shortest string of digits * that reconverts exactly to the given number, or to 'ilim' digits if * that will yield a shorter result. The denominator in David Gay's * conversion algorithm is known to be a power of 2**DIGIT_BIT, and hence * the division in the main loop may be replaced by a digit shift and * mask. * * Results: * Returns the string of significant decimal digits, in newly allocated * memory * * Side effects: * Stores the location of the decimal point in '*decpt' and the location * of the terminal null byte in '*endPtr'. * *---------------------------------------------------------------------- */ inline static char * ShorteningBignumConversionPowD( Double *dPtr, /* Original number to convert. */ int convType, /* Type of conversion (shortest, Steele, * E format, F format). */ Tcl_WideUInt bw, /* Integer significand. */ int b2, int b5, /* Scale factor for the significand in the * numerator. */ int m2plus, int m2minus, int m5, /* Scale factors for 1/2 ulp in the numerator * (will be different if bw == 1). */ int sd, /* Scale factor for the denominator. */ int k, /* Number of output digits before the decimal * point. */ int len, /* Number of digits to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Position of the terminal '\0' at * the end of the returned string. */ { char *retval = ckalloc(len + 1); /* Output buffer. */ mp_int b; /* Numerator of the fraction being * converted. */ mp_int mplus, mminus; /* Bounds for roundoff. */ mp_digit digit; /* Current output digit. */ char *s = retval; /* Cursor in the output buffer. */ int i; /* Index in the output buffer. */ mp_int temp; int r1; /* * b = bw * 2**b2 * 5**b5 * mminus = 5**m5 */ TclBNInitBignumFromWideUInt(&b, bw); mp_init_set_int(&mminus, 1); MulPow5(&b, b5, &b); mp_mul_2d(&b, b2, &b); /* * Adjust if the logarithm was guessed wrong. */ if (b.used <= sd) { mp_mul_d(&b, 10, &b); ++m2plus; ++m2minus; ++m5; ilim = ilim1; --k; } /* * mminus = 5**m5 * 2**m2minus * mplus = 5**m5 * 2**m2plus */ mp_mul_2d(&mminus, m2minus, &mminus); MulPow5(&mminus, m5, &mminus); if (m2plus > m2minus) { mp_init_copy(&mplus, &mminus); mp_mul_2d(&mplus, m2plus-m2minus, &mplus); } mp_init(&temp); /* * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT) * by mp_digit extraction. */ i = 0; for (;;) { if (b.used <= sd) { digit = 0; } else { digit = b.dp[sd]; if (b.used > sd+1 || digit >= 10) { Tcl_Panic("wrong digit!"); } --b.used; mp_clamp(&b); } /* * Does the current digit put us on the low side of the exact value * but within within roundoff of being exact? */ r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus); if (r1 == MP_LT || (r1 == MP_EQ && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { /* * Make sure we shouldn't be rounding *up* instead, in case the * next number above is closer. */ if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { ++digit; if (digit == 10) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } } /* * Stash the last digit. */ *s++ = '0' + digit; break; } /* * Does one plus the current digit put us within roundoff of the * number? */ if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd, convType, dPtr->w.word1 & 1, &temp)) { if (digit == 9) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } ++digit; *s++ = '0' + digit; break; } /* * Have we converted all the requested digits? */ *s++ = '0' + digit; if (i == ilim) { if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { s = BumpUp(s, retval, &k); } break; } /* * Advance to the next digit. */ mp_mul_d(&b, 10, &b); mp_mul_d(&mminus, 10, &mminus); if (m2plus > m2minus) { mp_mul_2d(&mminus, m2plus-m2minus, &mplus); } ++i; } /* * Endgame - store the location of the decimal point and the end of the * string. */ if (m2plus > m2minus) { mp_clear(&mplus); } mp_clear_multi(&b, &mminus, &temp, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * StrictBignumConversionPowD -- * * Converts a double-precision number to a fixed-lengt string of 'ilim' * digits (or 'ilim1' if log10(d) has been overestimated). The * denominator in David Gay's conversion algorithm is known to be a power * of 2**DIGIT_BIT, and hence the division in the main loop may be * replaced by a digit shift and mask. * * Results: * Returns the string of significant decimal digits, in newly allocated * memory. * * Side effects: * Stores the location of the decimal point in '*decpt' and the location * of the terminal null byte in '*endPtr'. * *---------------------------------------------------------------------- */ inline static char * StrictBignumConversionPowD( Double *dPtr, /* Original number to convert. */ int convType, /* Type of conversion (shortest, Steele, * E format, F format). */ Tcl_WideUInt bw, /* Integer significand. */ int b2, int b5, /* Scale factor for the significand in the * numerator. */ int sd, /* Scale factor for the denominator. */ int k, /* Number of output digits before the decimal * point. */ int len, /* Number of digits to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Position of the terminal '\0' at * the end of the returned string. */ { char *retval = ckalloc(len + 1); /* Output buffer. */ mp_int b; /* Numerator of the fraction being * converted. */ mp_digit digit; /* Current output digit. */ char *s = retval; /* Cursor in the output buffer. */ int i; /* Index in the output buffer. */ mp_int temp; /* * b = bw * 2**b2 * 5**b5 */ TclBNInitBignumFromWideUInt(&b, bw); MulPow5(&b, b5, &b); mp_mul_2d(&b, b2, &b); /* * Adjust if the logarithm was guessed wrong. */ if (b.used <= sd) { mp_mul_d(&b, 10, &b); ilim = ilim1; --k; } mp_init(&temp); /* * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT) * by mp_digit extraction. */ i = 1; for (;;) { if (b.used <= sd) { digit = 0; } else { digit = b.dp[sd]; if (b.used > sd+1 || digit >= 10) { Tcl_Panic("wrong digit!"); } --b.used; mp_clamp(&b); } /* * Have we converted all the requested digits? */ *s++ = '0' + digit; if (i == ilim) { if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { s = BumpUp(s, retval, &k); } while (*--s == '0') { /* do nothing */ } ++s; break; } /* * Advance to the next digit. */ mp_mul_d(&b, 10, &b); ++i; } /* * Endgame - store the location of the decimal point and the end of the * string. */ mp_clear_multi(&b, &temp, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * ShouldBankerRoundUp -- * * Tests whether a digit should be rounded up or down when finishing * bignum-based floating point conversion. * * Results: * Returns 1 if the number needs to be rounded up, 0 otherwise. * *---------------------------------------------------------------------- */ inline static int ShouldBankerRoundUp( mp_int *twor, /* 2x the remainder from thd division that * produced the last digit. */ mp_int *S, /* Denominator. */ int isodd) /* Flag == 1 if the last digit is odd. */ { int r = mp_cmp_mag(twor, S); switch (r) { case MP_LT: return 0; case MP_EQ: return isodd; case MP_GT: return 1; } Tcl_Panic("in ShouldBankerRoundUp, trichotomy fails!"); return 0; } /* *---------------------------------------------------------------------- * * ShouldBankerRoundUpToNext -- * * Tests whether the remainder is great enough to force rounding to the * next higher digit. * * Results: * Returns 1 if the number should be rounded up, 0 otherwise. * *---------------------------------------------------------------------- */ inline static int ShouldBankerRoundUpToNext( mp_int *b, /* Remainder from the division that produced * the last digit. */ mp_int *m, /* Numerator of the rounding tolerance. */ mp_int *S, /* Denominator. */ int convType, /* Conversion type: STEELE0 defeats * round-to-even. (Not sure why one would want * this; I coped it from Gay). FIXME */ int isodd, /* 1 if the integer significand is odd. */ mp_int *temp) /* Work area needed for the calculation. */ { int r; /* * Compare b and S-m: this is the same as comparing B+m and S. */ mp_add(b, m, temp); r = mp_cmp_mag(temp, S); switch(r) { case MP_LT: return 0; case MP_EQ: if (convType == TCL_DD_STEELE0) { return 0; } else { return isodd; } case MP_GT: return 1; } Tcl_Panic("in ShouldBankerRoundUpToNext, trichotomy fails!"); return 0; } /* *---------------------------------------------------------------------- * * ShorteningBignumConversion -- * * Convert a floating point number to a variable-length digit string * using the multiprecision method. * * Results: * Returns the string of digits. * * Side effects: * Stores the position of the decimal point in *decpt. Stores a pointer * to the end of the number in *endPtr. * *---------------------------------------------------------------------- */ inline static char * ShorteningBignumConversion( Double *dPtr, /* Original number being converted. */ int convType, /* Conversion type. */ Tcl_WideUInt bw, /* Integer significand and exponent. */ int b2, /* Scale factor for the significand. */ int m2plus, int m2minus, /* Scale factors for 1/2 ulp in numerator. */ int s2, int s5, /* Scale factors for denominator. */ int k, /* Guessed position of the decimal point. */ int len, /* Size of the digit buffer to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Pointer to the end of the number */ { char *retval = ckalloc(len+1); /* Buffer of digits to return. */ char *s = retval; /* Cursor in the return value. */ mp_int b; /* Numerator of the result. */ mp_int mminus; /* 1/2 ulp below the result. */ mp_int mplus; /* 1/2 ulp above the result. */ mp_int S; /* Denominator of the result. */ mp_int dig; /* Current digit of the result. */ int digit; /* Current digit of the result. */ mp_int temp; /* Work area. */ int minit = 1; /* Fudge factor for when we misguess k. */ int i; int r1; /* * b = bw * 2**b2 * 5**b5 * S = 2**s2 * 5*s5 */ TclBNInitBignumFromWideUInt(&b, bw); mp_mul_2d(&b, b2, &b); mp_init_set_int(&S, 1); MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S); /* * Handle the case where we guess the position of the decimal point wrong. */ if (mp_cmp_mag(&b, &S) == MP_LT) { mp_mul_d(&b, 10, &b); minit = 10; ilim =ilim1; --k; } /* * mminus = 2**m2minus * 5**m5 */ mp_init_set_int(&mminus, minit); mp_mul_2d(&mminus, m2minus, &mminus); if (m2plus > m2minus) { mp_init_copy(&mplus, &mminus); mp_mul_2d(&mplus, m2plus-m2minus, &mplus); } mp_init(&temp); /* * Loop through the digits. */ mp_init(&dig); i = 1; for (;;) { mp_div(&b, &S, &dig, &b); if (dig.used > 1 || dig.dp[0] >= 10) { Tcl_Panic("wrong digit!"); } digit = dig.dp[0]; /* * Does the current digit leave us with a remainder small enough to * round to it? */ r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus); if (r1 == MP_LT || (r1 == MP_EQ && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { ++digit; if (digit == 10) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } } *s++ = '0' + digit; break; } /* * Does the current digit leave us with a remainder large enough to * commit to rounding up to the next higher digit? */ if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType, dPtr->w.word1 & 1, &temp)) { ++digit; if (digit == 10) { *s++ = '9'; s = BumpUp(s, retval, &k); break; } *s++ = '0' + digit; break; } /* * Have we converted all the requested digits? */ *s++ = '0' + digit; if (i == ilim) { mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); } break; } /* * Advance to the next digit. */ if (s5 > 0) { /* * Can possibly shorten the denominator. */ mp_mul_2d(&b, 1, &b); mp_mul_2d(&mminus, 1, &mminus); if (m2plus > m2minus) { mp_mul_2d(&mplus, 1, &mplus); } mp_div_d(&S, 5, &S, NULL); --s5; /* * IDEA: It might possibly be a win to fall back to int64 * arithmetic here if S < 2**64/10. But it's a win only for * a fairly narrow range of magnitudes so perhaps not worth * bothering. We already know that we shorten the * denominator by at least 1 mp_digit, perhaps 2, as we do * the conversion for 17 digits of significance. * Possible savings: * 10**26 1 trip through loop before fallback possible * 10**27 1 trip * 10**28 2 trips * 10**29 3 trips * 10**30 4 trips * 10**31 5 trips * 10**32 6 trips * 10**33 7 trips * 10**34 8 trips * 10**35 9 trips * 10**36 10 trips * 10**37 11 trips * 10**38 12 trips * 10**39 13 trips * 10**40 14 trips * 10**41 15 trips * 10**42 16 trips * thereafter no gain. */ } else { mp_mul_d(&b, 10, &b); mp_mul_d(&mminus, 10, &mminus); if (m2plus > m2minus) { mp_mul_2d(&mplus, 10, &mplus); } } ++i; } /* * Endgame - store the location of the decimal point and the end of the * string. */ if (m2plus > m2minus) { mp_clear(&mplus); } mp_clear_multi(&b, &mminus, &temp, &dig, &S, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * StrictBignumConversion -- * * Convert a floating point number to a fixed-length digit string using * the multiprecision method. * * Results: * Returns the string of digits. * * Side effects: * Stores the position of the decimal point in *decpt. Stores a pointer * to the end of the number in *endPtr. * *---------------------------------------------------------------------- */ inline static char * StrictBignumConversion( Double *dPtr, /* Original number being converted. */ int convType, /* Conversion type. */ Tcl_WideUInt bw, /* Integer significand and exponent. */ int b2, /* Scale factor for the significand. */ int s2, int s5, /* Scale factors for denominator. */ int k, /* Guessed position of the decimal point. */ int len, /* Size of the digit buffer to allocate. */ int ilim, /* Number of digits to convert if b >= s */ int ilim1, /* Number of digits to convert if b < s */ int *decpt, /* OUTPUT: Position of the decimal point. */ char **endPtr) /* OUTPUT: Pointer to the end of the number */ { char *retval = ckalloc(len+1); /* Buffer of digits to return. */ char *s = retval; /* Cursor in the return value. */ mp_int b; /* Numerator of the result. */ mp_int S; /* Denominator of the result. */ mp_int dig; /* Current digit of the result. */ int digit; /* Current digit of the result. */ mp_int temp; /* Work area. */ int g; /* Size of the current digit ground. */ int i, j; /* * b = bw * 2**b2 * 5**b5 * S = 2**s2 * 5*s5 */ mp_init_multi(&temp, &dig, NULL); TclBNInitBignumFromWideUInt(&b, bw); mp_mul_2d(&b, b2, &b); mp_init_set_int(&S, 1); MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S); /* * Handle the case where we guess the position of the decimal point wrong. */ if (mp_cmp_mag(&b, &S) == MP_LT) { mp_mul_d(&b, 10, &b); ilim =ilim1; --k; } /* * Convert the leading digit. */ i = 0; mp_div(&b, &S, &dig, &b); if (dig.used > 1 || dig.dp[0] >= 10) { Tcl_Panic("wrong digit!"); } digit = dig.dp[0]; /* * Is a single digit all that was requested? */ *s++ = '0' + digit; if (++i >= ilim) { mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); } } else { for (;;) { /* * Shift by a group of digits. */ g = ilim - i; if (g > DIGIT_GROUP) { g = DIGIT_GROUP; } if (s5 >= g) { mp_div_d(&S, dpow5[g], &S, NULL); s5 -= g; } else if (s5 > 0) { mp_div_d(&S, dpow5[s5], &S, NULL); mp_mul_d(&b, dpow5[g - s5], &b); s5 = 0; } else { mp_mul_d(&b, dpow5[g], &b); } mp_mul_2d(&b, g, &b); /* * As with the shortening bignum conversion, it's possible at this * point that we will have reduced the denominator to less than * 2**64/10, at which point it would be possible to fall back to * to int64 arithmetic. But the potential payoff is tremendously * less - unless we're working in F format - because we know that * three groups of digits will always suffice for %#.17e, the * longest format that doesn't introduce empty precision. * * Extract the next group of digits. */ mp_div(&b, &S, &dig, &b); if (dig.used > 1) { Tcl_Panic("wrong digit!"); } digit = dig.dp[0]; for (j = g-1; j >= 0; --j) { int t = itens[j]; *s++ = digit / t + '0'; digit %= t; } i += g; /* * Have we converted all the requested digits? */ if (i == ilim) { mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); } break; } } } while (*--s == '0') { /* do nothing */ } ++s; /* * Endgame - store the location of the decimal point and the end of the * string. */ mp_clear_multi(&b, &S, &temp, &dig, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; } /* *---------------------------------------------------------------------- * * TclDoubleDigits -- * * Core of Tcl's conversion of double-precision floating point numbers to * decimal. * * Results: * Returns a newly-allocated string of digits. * * Side effects: * Sets *decpt to the index of the character in the string before the * place that the decimal point should go. If 'endPtr' is not NULL, sets * endPtr to point to the terminating '\0' byte of the string. Sets *sign * to 1 if a minus sign should be printed with the number, or 0 if a plus * sign (or no sign) should appear. * * This function is a service routine that produces the string of digits for * floating-point-to-decimal conversion. It can do a number of things * according to the 'flags' argument. Valid values for 'flags' include: * TCL_DD_SHORTEST - This is the default for floating point conversion if * ::tcl_precision is 0. It constructs the shortest string of * digits that will reconvert to the given number when scanned. * For floating point numbers that are exactly between two * decimal numbers, it resolves using the 'round to even' rule. * With this value, the 'ndigits' parameter is ignored. * TCL_DD_STEELE - This value is not recommended and may be removed in * the future. It follows the conversion algorithm outlined in * "How to Print Floating-Point Numbers Accurately" by Guy * L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, * pp. 112-126]. This rule has the effect of rendering 1e23 as * 9.9999999999999999e22 - which is a 'better' approximation in * the sense that it will reconvert correctly even if a * subsequent input conversion is 'round up' or 'round down' * rather than 'round to nearest', but is surprising otherwise. * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e format * conversion (or for default floating->string if tcl_precision * is not 0). It constructs a string of at most 'ndigits' digits, * choosing the one that is closest to the given number (and * resolving ties with 'round to even'). It is allowed to return * fewer than 'ndigits' if the number converts exactly; if the * TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG is supplied instead, it * also returns fewer digits if the shorter string will still * reconvert without loss to the given input number. In any case, * strings of trailing zeroes are suppressed. * TCL_DD_F_FORMAT - This value is used to prepare numbers for %f format * conversion. It requests that conversion proceed until * 'ndigits' digits after the decimal point have been converted. * It is possible for this format to result in a zero-length * string if the number is sufficiently small. Again, it is * permissible for TCL_DD_F_FORMAT to return fewer digits for a * number that converts exactly, and changing the argument to * TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow the routine * also to return fewer digits if the shorter string will still * reconvert without loss to the given input number. Strings of * trailing zeroes are suppressed. * * To any of these flags may be OR'ed TCL_DD_NO_QUICK; this flag requires * all calculations to be done in exact arithmetic. Normally, E and F * format with fewer than about 14 digits will be done with a quick * floating point approximation and fall back on the exact arithmetic * only if the input number is close enough to the midpoint between two * decimal strings that more precision is needed to resolve which string * is correct. * * The value stored in the 'decpt' argument on return may be negative * (indicating that the decimal point falls to the left of the string) or * greater than the length of the string. In addition, the value -9999 is used * as a sentinel to indicate that the string is one of the special values * "Infinity" and "NaN", and that no decimal point should be inserted. * *---------------------------------------------------------------------- */ char * TclDoubleDigits( double dv, /* Number to convert. */ int ndigits, /* Number of digits requested. */ int flags, /* Conversion flags. */ int *decpt, /* OUTPUT: Position of the decimal point. */ int *sign, /* OUTPUT: 1 if the result is negative. */ char **endPtr) /* OUTPUT: If not NULL, receives a pointer to * one character beyond the end of the * returned string. */ { int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK); /* Type of conversion being performed: * TCL_DD_SHORTEST0, TCL_DD_STEELE0, * TCL_DD_E_FORMAT, or TCL_DD_F_FORMAT. */ Double d; /* Union for deconstructing doubles. */ Tcl_WideUInt bw; /* Integer significand. */ int be; /* Power of 2 by which b must be multiplied */ int bbits; /* Number of bits needed to represent b. */ int denorm; /* Flag == 1 iff the input number was * denormalized. */ int k; /* Estimate of floor(log10(d)). */ int k_check; /* Flag == 1 if d is near enough to a power of * ten that k must be checked. */ int b2, b5, s2, s5; /* Powers of 2 and 5 in the numerator and * denominator of intermediate results. */ int ilim = -1, ilim1 = -1; /* Number of digits to convert, and number to * convert if log10(d) has been * overestimated. */ char *retval; /* Return value from this function. */ int i = -1; /* * Put the input number into a union for bit-whacking. */ d.d = dv; /* * Handle the cases of negative numbers (by taking the absolute value: * this includes -Inf and -NaN!), infinity, Not a Number, and zero. */ TakeAbsoluteValue(&d, sign); if ((d.w.word0 & EXP_MASK) == EXP_MASK) { return FormatInfAndNaN(&d, decpt, endPtr); } if (d.d == 0.0) { return FormatZero(decpt, endPtr); } /* * Unpack the floating point into a wide integer and an exponent. * Determine the number of bits that the big integer requires, and compute * a quick approximation (which may be one too high) of ceil(log10(d.d)). */ denorm = ((d.w.word0 & EXP_MASK) == 0); DoubleToExpAndSig(d.d, &bw, &be, &bbits); k = ApproximateLog10(bw, be, bbits); k = BetterLog10(d.d, k, &k_check); /* At this point, we have: * d is the number to convert. * bw are significand and exponent: d == bw*2**be, * bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits * k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0 if we * know that k is exactly ceil(log10(d)) and 1 if we need to check. * We want a rational number * r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5), * with b2, b5, s2, s5 >= 0. Note that the most significant decimal * digit is floor(r) and that successive digits can be obtained by * setting r <- 10*floor(r) (or b <= 10 * (b % S)). Find appropriate * b2, b5, s2, s5. */ ComputeScale(be, k, &b2, &b5, &s2, &s5); /* * Correct an incorrect caller-supplied 'ndigits'. Also determine: * i = The maximum number of decimal digits that will be returned in the * formatted string. This is k + 1 + ndigits for F format, 18 for * shortest and Steele, and ndigits for E format. * ilim = The number of significant digits to convert if k has been * guessed correctly. This is -1 for shortest and Steele (which * stop when all significance has been lost), 'ndigits' for E * format, and 'k + 1 + ndigits' for F format. * ilim1 = The minimum number of significant digits to convert if k has * been guessed 1 too high. This, too, is -1 for shortest and * Steele, and 'ndigits' for E format, but it's 'ndigits-1' for F * format. */ SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1); /* * Try to do low-precision conversion in floating point rather than * resorting to expensive multiprecision arithmetic. */ if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) { retval = QuickConversion(d.d, k, k_check, flags, i, ilim, ilim1, decpt, endPtr); if (retval != NULL) { return retval; } } /* * For shortening conversions, determine the upper and lower bounds for * the remainder at which we can stop. * m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the high * side, and * m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the low * side. * We may need to increase s2 to put m2plus, m2minus, b2 over a common * denominator. */ if (flags & TCL_DD_SHORTEN_FLAG) { int m2minus = b2; int m2plus; int m5 = b5; int len = i; /* * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5) is 1/2 unit * in the least significant place of the floating point number. */ if (denorm) { i = be + EXPONENT_BIAS + (FP_PRECISION-1); } else { i = 1 + FP_PRECISION - bbits; } b2 += i; s2 += i; /* * Reduce the fractions to lowest terms, since the above calculation * may have left excess powers of 2 in numerator and denominator. */ CastOutPowersOf2(&b2, &m2minus, &s2); /* * In the special case where bw==1, the nearest floating point number * to it on the low side is 1/4 ulp below it. Adjust accordingly. */ m2plus = m2minus; if (!denorm && bw == 1) { ++b2; ++s2; ++m2plus; } if (s5+1 < N_LOG2POW5 && s2+1 + log2pow5[s5+1] <= 64) { /* * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit word, * then all our intermediate calculations can be done using exact * 64-bit arithmetic with no need for expensive multiprecision * operations. (This will be true for all numbers in the range * [1.0e-3 .. 1.0e+24]). */ return ShorteningInt64Conversion(&d, convType, bw, b2, b5, m2plus, m2minus, m5, s2, s5, k, len, ilim, ilim1, decpt, endPtr); } else if (s5 == 0) { /* * The denominator is a power of 2, so we can replace division by * digit shifts. First we round up s2 to a multiple of DIGIT_BIT, * and adjust m2 and b2 accordingly. Then we launch into a version * of the comparison that's specialized for the 'power of mp_digit * in the denominator' case. */ if (s2 % DIGIT_BIT != 0) { int delta = DIGIT_BIT - (s2 % DIGIT_BIT); b2 += delta; m2plus += delta; m2minus += delta; s2 += delta; } return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5, m2plus, m2minus, m5, s2/DIGIT_BIT, k, len, ilim, ilim1, decpt, endPtr); } else { /* * Alas, there's no helpful special case; use full-up bignum * arithmetic for the conversion. */ return ShorteningBignumConversion(&d, convType, bw, b2, m2plus, m2minus, s2, s5, k, len, ilim, ilim1, decpt, endPtr); } } else { /* * Non-shortening conversion. */ int len = i; /* * Reduce numerator and denominator to lowest terms. */ if (b2 >= s2 && s2 > 0) { b2 -= s2; s2 = 0; } else if (s2 >= b2 && b2 > 0) { s2 -= b2; b2 = 0; } if (s5+1 < N_LOG2POW5 && s2+1 + log2pow5[s5+1] <= 64) { /* * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit word, * then all our intermediate calculations can be done using exact * 64-bit arithmetic with no need for expensive multiprecision * operations. */ return StrictInt64Conversion(&d, convType, bw, b2, b5, s2, s5, k, len, ilim, ilim1, decpt, endPtr); } else if (s5 == 0) { /* * The denominator is a power of 2, so we can replace division by * digit shifts. First we round up s2 to a multiple of DIGIT_BIT, * and adjust m2 and b2 accordingly. Then we launch into a version * of the comparison that's specialized for the 'power of mp_digit * in the denominator' case. */ if (s2 % DIGIT_BIT != 0) { int delta = DIGIT_BIT - (s2 % DIGIT_BIT); b2 += delta; s2 += delta; } return StrictBignumConversionPowD(&d, convType, bw, b2, b5, s2/DIGIT_BIT, k, len, ilim, ilim1, decpt, endPtr); } else { /* * There are no helpful special cases, but at least we know in * advance how many digits we will convert. We can run the * conversion in steps of DIGIT_GROUP digits, so as to have many * fewer mp_int divisions. */ return StrictBignumConversion(&d, convType, bw, b2, s2, s5, k, len, ilim, ilim1, decpt, endPtr); } } } /* *---------------------------------------------------------------------- * * 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 #if defined(__sgi) && defined(_COMPILER_VERSION) union fpc_csr mipsCR; mipsCR.fc_word = get_fpc_csr(); mipsCR.fc_struct.flush = 0; set_fpc_csr(mipsCR.fc_word); #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 = 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) { pow10vals[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); } mp_init_set_int(pow5_13, 1220703125); for (i = 1; i < 5; ++i) { mp_init(pow5_13 + i); mp_sqr(pow5_13 + i - 1, pow5_13 + i); } /* * 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. */ 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.)); 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; ckfree(pow10_wide); for (i=0; i<9; ++i) { mp_clear(pow5 + i); } for (i=0; i < 5; ++i) { mp_clear(pow5_13 + 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( const mp_int *a) /* Integer to convert. */ { mp_int b; int bits, shift, i, lsb; double r; /* * We need a 'mantBits'-bit significand. Determine what shift will * give us that. */ 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 - bits; /* * If shift > 0, shift the significand left by the requisite number of * bits. If shift == 0, the significand is already exactly 'mantBits' * in length. If shift < 0, we will need to shift the significand right * by the requisite number of bits, and round it. If the '1-shift' * least significant bits are 0, but the 'shift'th bit is nonzero, * then the significand lies exactly between two values and must be * 'rounded to even'. */ mp_init(&b); if (shift == 0) { mp_copy(a, &b); } else if (shift > 0) { mp_mul_2d(a, shift, &b); } else if (shift < 0) { lsb = mp_cnt_lsb(a); if (lsb == -1-shift) { /* * Round to even */ mp_div_2d(a, -shift, &b, NULL); if (mp_isodd(&b)) { if (b.sign == MP_ZPOS) { mp_add_d(&b, 1, &b); } else { mp_sub_d(&b, 1, &b); } } } else { /* * Ordinary rounding */ mp_div_2d(a, -1-shift, &b, NULL); if (b.sign == MP_ZPOS) { mp_add_d(&b, 1, &b); } else { mp_sub_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; } } /* *---------------------------------------------------------------------- * * TclCeil -- * * Computes the smallest floating point number that is at least the * mp_int argument. * * Results: * Returns the floating point number. * *---------------------------------------------------------------------- */ double TclCeil( const 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; } /* *---------------------------------------------------------------------- * * TclFloor -- * * Computes the largest floating point number less than or equal to the * mp_int argument. * * Results: * Returns the floating point value. * *---------------------------------------------------------------------- */ double TclFloor( const 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( const 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 * pow10vals[exponent&0xf], &j); expt += j; for (i=4; i<9; ++i) { if (exponent & (1<> 32) & 0xffffffff) | (w << 32)); } #endif /* *---------------------------------------------------------------------- * * 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: */