diff options
Diffstat (limited to 'generic/tclStrToD.c')
| -rwxr-xr-x | generic/tclStrToD.c | 5018 |
1 files changed, 0 insertions, 5018 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c deleted file mode 100755 index 332cfca..0000000 --- a/generic/tclStrToD.c +++ /dev/null @@ -1,5018 +0,0 @@ -/* - * 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 <math.h> - -/* - * Define KILL_OCTAL to suppress interpretation of numbers with leading zero - * as octal. (Ceterum censeo: numeros octonarios delendos esse.) - */ - -#undef KILL_OCTAL - -/* - * This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754 - * floating point; of these, only IEEE-754 can represent NaN. IEEE-754 can be - * uniquely determined by radix and by the widths of significand and exponent. - */ - -#if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024) -# define IEEE_FLOATING_POINT -#endif - -/* - * 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 <sunmath.h> -#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 <sys/fpu.h> -#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 int mantDIGIT; /* Number of mp_digit's needed to hold the - * significand of a double. */ -static const double pow_10_2_n[] = { /* Inexact higher powers of ten. */ - 1.0, - 100.0, - 10000.0, - 1.0e+8, - 1.0e+16, - 1.0e+32, - 1.0e+64, - 1.0e+128, - 1.0e+256 -}; - -static int n770_fp; /* Flag is 1 on Nokia N770 floating point. - * Nokia's floating point has the words - * reversed: if big-endian is 7654 3210, - * and little-endian is 0123 4567, - * then Nokia's FP is 4567 0123; - * little-endian within the 32-bit words but - * big-endian between them. */ - -/* - * 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, BAD_OCTAL, DECIMAL, - LEADING_RADIX_POINT, FRACTION, - EXPONENT_START, EXPONENT_SIGNUM, EXPONENT, - sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY -#ifdef IEEE_FLOATING_POINT - , sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH -#endif - } state = INITIAL; - enum State acceptState = INITIAL; - - int signum = 0; /* Sign of the number being parsed. */ - Tcl_WideUInt significandWide = 0; - /* Significand of the number being parsed (if - * no overflow). */ - mp_int significandBig; /* Significand of the number being parsed (if - * it overflows significandWide). */ - int significandOverflow = 0;/* Flag==1 iff significandBig is used. */ - Tcl_WideUInt octalSignificandWide = 0; - /* Significand of an octal number; needed - * because we don't know whether a number with - * a leading zero is octal or decimal until - * we've scanned forward to a '.' or 'e'. */ - mp_int octalSignificandBig; /* Significand of octal number once - * octalSignificandWide overflows. */ - int octalSignificandOverflow = 0; - /* Flag==1 if octalSignificandBig is used. */ - int numSigDigs = 0; /* Number of significant digits in the decimal - * significand. */ - int numTrailZeros = 0; /* Number of trailing zeroes at the current - * point in the parse. */ - int numDigitsAfterDp = 0; /* Number of digits scanned after the decimal - * point. */ - int exponentSignum = 0; /* Signum of the exponent of a floating point - * number. */ - long exponent = 0; /* Exponent of a floating point number. */ - const char *p; /* Pointer to next character to scan. */ - size_t len; /* Number of characters remaining after p. */ - const char *acceptPoint; /* Pointer to position after last character in - * an acceptable number. */ - size_t acceptLen; /* Number of characters following that - * point. */ - int status = TCL_OK; /* Status to return to caller. */ - char d = 0; /* Last hexadecimal digit scanned; initialized - * to avoid a compiler warning. */ - int shift = 0; /* Amount to shift when accumulating binary */ - int explicitOctal = 0; - -#define ALL_BITS (~(Tcl_WideUInt)0) -#define MOST_BITS (ALL_BITS >> 1) - - /* - * Initialize bytes to start of the object's string rep if the caller - * didn't pass anything else. - */ - - if (bytes == NULL) { - bytes = TclGetString(objPtr); - } - - p = bytes; - len = numBytes; - acceptPoint = p; - acceptLen = len; - while (1) { - char c = len ? *p : '\0'; - switch (state) { - - case INITIAL: - /* - * Initial state. Acceptable characters are +, -, digits, period, - * I, N, and whitespace. - */ - - if (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') { - explicitOctal = 1; - state = ZERO_O; - break; - } -#ifdef KILL_OCTAL - goto decimal; -#endif - /* FALLTHROUGH */ - - case OCTAL: - /* - * Scanned an optional + or -, followed by a string of octal - * digits. Acceptable inputs are more digits, period, or E. If 8 - * or 9 is encountered, commit to floating point. - */ - - acceptState = state; - acceptPoint = p; - acceptLen = len; - /* FALLTHROUGH */ - case ZERO_O: - zeroo: - if (c == '0') { - numTrailZeros++; - state = OCTAL; - break; - } else if (c >= '1' && c <= '7') { - if (objPtr != NULL) { - shift = 3 * (numTrailZeros + 1); - significandOverflow = AccumulateDecimalDigit( - (unsigned)(c-'0'), numTrailZeros, - &significandWide, &significandBig, - significandOverflow); - - if (!octalSignificandOverflow) { - /* - * Shifting by more bits than are in the value being - * shifted is at least de facto nonportable. Check for - * too large shifts first. - */ - - if ((octalSignificandWide != 0) - && (((size_t)shift >= - CHAR_BIT*sizeof(Tcl_WideUInt)) - || (octalSignificandWide > - (~(Tcl_WideUInt)0 >> shift)))) { - octalSignificandOverflow = 1; - TclBNInitBignumFromWideUInt(&octalSignificandBig, - octalSignificandWide); - } - } - if (!octalSignificandOverflow) { - octalSignificandWide = - (octalSignificandWide << shift) + (c - '0'); - } else { - mp_mul_2d(&octalSignificandBig, shift, - &octalSignificandBig); - mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'), - &octalSignificandBig); - } - } - if (numSigDigs != 0) { - numSigDigs += numTrailZeros+1; - } else { - numSigDigs = 1; - } - numTrailZeros = 0; - state = OCTAL; - break; - } - /* FALLTHROUGH */ - - case BAD_OCTAL: - if (explicitOctal) { - /* - * No forgiveness for bad digits in explicitly octal numbers. - */ - - goto endgame; - } - if (flags & TCL_PARSE_INTEGER_ONLY) { - /* - * No seeking floating point when parsing only integer. - */ - - goto endgame; - } -#ifndef KILL_OCTAL - - /* - * Scanned a number with a leading zero that contains an 8, 9, - * radix point or E. This is an invalid octal number, but might - * still be floating point. - */ - - if (c == '0') { - numTrailZeros++; - state = BAD_OCTAL; - break; - } else if (isdigit(UCHAR(c))) { - if (objPtr != NULL) { - significandOverflow = AccumulateDecimalDigit( - (unsigned)(c-'0'), numTrailZeros, - &significandWide, &significandBig, - significandOverflow); - } - if (numSigDigs != 0) { - numSigDigs += (numTrailZeros + 1); - } else { - numSigDigs = 1; - } - numTrailZeros = 0; - state = BAD_OCTAL; - break; - } else if (c == '.') { - state = FRACTION; - break; - } else if (c == 'E' || c == 'e') { - state = EXPONENT_START; - break; - } -#endif - goto endgame; - - /* - * Scanned 0x. If state is HEXADECIMAL, scanned at least one - * character following the 0x. The only acceptable inputs are - * hexadecimal digits. - */ - - case HEXADECIMAL: - acceptState = state; - acceptPoint = p; - acceptLen = len; - /* FALLTHROUGH */ - - case ZERO_X: - zerox: - if (c == '0') { - numTrailZeros++; - state = HEXADECIMAL; - break; - } else if (isdigit(UCHAR(c))) { - d = (c-'0'); - } else if (c >= 'A' && c <= 'F') { - d = (c-'A'+10); - } else if (c >= 'a' && c <= 'f') { - d = (c-'a'+10); - } else { - goto endgame; - } - if (objPtr != NULL) { - shift = 4 * (numTrailZeros + 1); - if (!significandOverflow) { - /* - * Shifting by more bits than are in the value being - * shifted is at least de facto nonportable. Check for too - * large shifts first. - */ - - if (significandWide != 0 && - ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || - significandWide > (~(Tcl_WideUInt)0 >> shift))) { - significandOverflow = 1; - TclBNInitBignumFromWideUInt(&significandBig, - significandWide); - } - } - if (!significandOverflow) { - significandWide = (significandWide << shift) + d; - } else { - mp_mul_2d(&significandBig, shift, &significandBig); - mp_add_d(&significandBig, (mp_digit) d, &significandBig); - } - } - numTrailZeros = 0; - state = HEXADECIMAL; - break; - - case BINARY: - acceptState = state; - acceptPoint = p; - acceptLen = len; - case ZERO_B: - 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. - */ - -#ifdef KILL_OCTAL - decimal: -#endif - acceptState = state; - acceptPoint = p; - acceptLen = len; - if (c == '0') { - numTrailZeros++; - state = DECIMAL; - break; - } else if (isdigit(UCHAR(c))) { - if (objPtr != NULL) { - significandOverflow = AccumulateDecimalDigit( - (unsigned)(c - '0'), numTrailZeros, - &significandWide, &significandBig, - significandOverflow); - } - numSigDigs += numTrailZeros+1; - numTrailZeros = 0; - state = DECIMAL; - break; - } else if (flags & TCL_PARSE_INTEGER_ONLY) { - goto endgame; - } else if (c == '.') { - state = FRACTION; - break; - } else if (c == 'E' || c == 'e') { - state = EXPONENT_START; - break; - } - goto endgame; - - /* - * Found a decimal point. If no digits have yet been scanned, E is - * not allowed; otherwise, it introduces the exponent. If at least - * one digit has been found, we have a possible complete number. - */ - - case FRACTION: - acceptState = state; - acceptPoint = p; - acceptLen = len; - if (c == 'E' || c=='e') { - state = EXPONENT_START; - break; - } - /* FALLTHROUGH */ - - case LEADING_RADIX_POINT: - if (c == '0') { - numDigitsAfterDp++; - numTrailZeros++; - state = FRACTION; - break; - } else if (isdigit(UCHAR(c))) { - numDigitsAfterDp++; - if (objPtr != NULL) { - significandOverflow = AccumulateDecimalDigit( - (unsigned)(c-'0'), numTrailZeros, - &significandWide, &significandBig, - significandOverflow); - } - if (numSigDigs != 0) { - numSigDigs += numTrailZeros+1; - } else { - numSigDigs = 1; - } - numTrailZeros = 0; - state = FRACTION; - break; - } - goto endgame; - - case EXPONENT_START: - /* - * Scanned the E at the start of an exponent. Make sure a legal - * character follows before using the C library strtol routine, - * which allows whitespace. - */ - - if (c == '+') { - state = EXPONENT_SIGNUM; - break; - } else if (c == '-') { - exponentSignum = 1; - state = EXPONENT_SIGNUM; - break; - } - /* FALLTHROUGH */ - - case EXPONENT_SIGNUM: - /* - * Found the E at the start of the exponent, followed by a sign - * character. - */ - - if (isdigit(UCHAR(c))) { - exponent = c - '0'; - state = EXPONENT; - break; - } - goto endgame; - - case EXPONENT: - /* - * Found an exponent with at least one digit. Accumulate it, - * making sure to hard-pin it to LONG_MAX on overflow. - */ - - acceptState = state; - acceptPoint = p; - acceptLen = len; - if (isdigit(UCHAR(c))) { - if (exponent < (LONG_MAX - 9) / 10) { - exponent = 10 * exponent + (c - '0'); - } else { - exponent = LONG_MAX; - } - state = EXPONENT; - break; - } - goto endgame; - - /* - * Parse out INFINITY by simply spelling it out. INF is accepted - * as an abbreviation; other prefices are not. - */ - - case sI: - if (c == 'n' || c == 'N') { - state = sIN; - break; - } - goto endgame; - case sIN: - if (c == 'f' || c == 'F') { - state = sINF; - break; - } - goto endgame; - case sINF: - acceptState = state; - acceptPoint = p; - acceptLen = len; - if (c == 'i' || c == 'I') { - state = sINFI; - break; - } - goto endgame; - case sINFI: - if (c == 'n' || c == 'N') { - state = sINFIN; - break; - } - goto endgame; - case sINFIN: - if (c == 'i' || c == 'I') { - state = sINFINI; - break; - } - goto endgame; - case sINFINI: - if (c == 't' || c == 'T') { - state = sINFINIT; - break; - } - goto endgame; - case sINFINIT: - if (c == 'y' || c == 'Y') { - state = sINFINITY; - break; - } - goto endgame; - - /* - * Parse NaN's. - */ -#ifdef IEEE_FLOATING_POINT - case sN: - if (c == 'a' || c == 'A') { - state = sNA; - break; - } - goto endgame; - case sNA: - if (c == 'n' || c == 'N') { - state = sNAN; - break; - } - goto endgame; - case sNAN: - acceptState = state; - acceptPoint = p; - acceptLen = len; - if (c == '(') { - state = sNANPAREN; - break; - } - goto endgame; - - /* - * Parse NaN(hexdigits) - */ - case sNANHEX: - if (c == ')') { - state = sNANFINISH; - break; - } - /* FALLTHROUGH */ - case sNANPAREN: - if (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 BAD_OCTAL: - case ZERO_X: - case ZERO_O: - case ZERO_B: - case LEADING_RADIX_POINT: - case EXPONENT_START: - case EXPONENT_SIGNUM: - case sI: - case sIN: - case sINFI: - case sINFIN: - case sINFINI: - case sINFINIT: -#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 NO_WIDE_TYPE - if (octalSignificandWide <= (MOST_BITS + signum)) { - objPtr->typePtr = &tclWideIntType; - if (signum) { - objPtr->internalRep.wideValue = - - (Tcl_WideInt) octalSignificandWide; - } else { - objPtr->internalRep.wideValue = - (Tcl_WideInt) octalSignificandWide; - } - break; - } -#endif - TclBNInitBignumFromWideUInt(&octalSignificandBig, - octalSignificandWide); - octalSignificandOverflow = 1; - } else { - objPtr->typePtr = &tclIntType; - if (signum) { - objPtr->internalRep.longValue = - - (long) octalSignificandWide; - } else { - objPtr->internalRep.longValue = - (long) octalSignificandWide; - } - } - } - if (octalSignificandOverflow) { - if (signum) { - mp_neg(&octalSignificandBig, &octalSignificandBig); - } - TclSetBignumIntRep(objPtr, &octalSignificandBig); - } - break; - - case ZERO: - case DECIMAL: - significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1, - &significandWide, &significandBig, significandOverflow); - if (!significandOverflow && (significandWide > MOST_BITS+signum)){ - significandOverflow = 1; - TclBNInitBignumFromWideUInt(&significandBig, significandWide); - } - returnInteger: - if (!significandOverflow) { - if (significandWide > - (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { -#ifndef NO_WIDE_TYPE - if (significandWide <= MOST_BITS+signum) { - objPtr->typePtr = &tclWideIntType; - if (signum) { - objPtr->internalRep.wideValue = - - (Tcl_WideInt) significandWide; - } else { - objPtr->internalRep.wideValue = - (Tcl_WideInt) significandWide; - } - break; - } -#endif - TclBNInitBignumFromWideUInt(&significandBig, - significandWide); - significandOverflow = 1; - } else { - objPtr->typePtr = &tclIntType; - if (signum) { - objPtr->internalRep.longValue = - - (long) significandWide; - } else { - objPtr->internalRep.longValue = - (long) significandWide; - } - } - } - if (significandOverflow) { - if (signum) { - mp_neg(&significandBig, &significandBig); - } - TclSetBignumIntRep(objPtr, &significandBig); - } - break; - - case FRACTION: - case EXPONENT: - - /* - * Here, we're parsing a floating-point number. 'significandWide' - * or 'significandBig' contains the exact significand, according - * to whether 'significandOverflow' is set. The desired floating - * point value is significand * 10**k, where - * k = numTrailZeros+exponent-numDigitsAfterDp. - */ - - objPtr->typePtr = &tclDoubleType; - if (exponentSignum) { - exponent = -exponent; - } - if (!significandOverflow) { - objPtr->internalRep.doubleValue = MakeLowPrecisionDouble( - signum, significandWide, numSigDigs, - numTrailZeros + exponent - numDigitsAfterDp); - } else { - objPtr->internalRep.doubleValue = MakeHighPrecisionDouble( - signum, &significandBig, numSigDigs, - numTrailZeros + exponent - numDigitsAfterDp); - } - break; - - case sINF: - case sINFINITY: - if (signum) { - objPtr->internalRep.doubleValue = -HUGE_VAL; - } else { - objPtr->internalRep.doubleValue = HUGE_VAL; - } - objPtr->typePtr = &tclDoubleType; - break; - -#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); - if (state == BAD_OCTAL) { - Tcl_AppendToObj(msg, " (looks like invalid octal number)", -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<minincr) { - quot = minincr; - } - - mp_clear(&twoMd); - mp_clear(&twoMv); - - return approxResult + quot; -} - -/* - *---------------------------------------------------------------------- - * - * MultPow5 -- - * - * Multiply a bignum by a power of 5. - * - * Side effects: - * Stores base*5**n in result. - * - *---------------------------------------------------------------------- - */ - -inline static void -MulPow5( - mp_int *base, /* Number to multiply. */ - unsigned n, /* Power of 5 to multiply by. */ - mp_int *result) /* Place to store the result. */ -{ - mp_int *p = base; - int n13 = n / 13; - int r = n % 13; - - if (r != 0) { - mp_mul_d(p, dpow5[r], result); - p = result; - } - r = 0; - while (n13 != 0) { - if (n13 & 1) { - mp_mul(p, pow5_13+r, result); - p = result; - } - n13 >>= 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.)); - mantDIGIT = (mantBits + DIGIT_BIT-1) / DIGIT_BIT; - log10_DIGIT_MAX = (int) floor(DIGIT_BIT * log(2.) / log(10.)); - - /* - * Nokia 770's software-emulated floating point is "middle endian": the - * bytes within a 32-bit word are little-endian (like the native - * integers), but the two words of a 'double' are presented most - * significant word first. - */ - -#ifdef IEEE_FLOATING_POINT - bitwhack.dv = 1.000000238418579; - /* 3ff0 0000 4000 0000 */ - if ((bitwhack.iv >> 32) == 0x3ff00000) { - n770_fp = 0; - } else if ((bitwhack.iv & 0xffffffff) == 0x3ff00000) { - n770_fp = 1; - } else { - Tcl_Panic("unknown floating point word order on this machine"); - } -#endif -} - -/* - *---------------------------------------------------------------------- - * - * TclFinalizeDoubleConversion -- - * - * Cleans up this file on exit. - * - * Results: - * None - * - * Side effects: - * Memory allocated by TclInitDoubleConversion is freed. - * - *---------------------------------------------------------------------- - */ - -void -TclFinalizeDoubleConversion(void) -{ - int i; - - 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<<i)) { - retval = frexp(retval * pow_10_2_n[i], &j); - expt += j; - } - } - } else if (exponent < 0) { - /* - * Divide by 10**-exponent. - */ - - retval = frexp(retval / pow10vals[(-exponent) & 0xf], &j); - expt += j; - for (i=4; i<9; ++i) { - if ((-exponent) & (1<<i)) { - retval = frexp(retval / pow_10_2_n[i], &j); - expt += j; - } - } - } - - *machexp = expt; - return retval; -} - -/* - *---------------------------------------------------------------------- - * - * SafeLdExp -- - * - * Do an 'ldexp' operation, but handle denormals gracefully. - * - * Results: - * Returns the appropriately scaled value. - * - * On some platforms, 'ldexp' fails when presented with a number too - * small to represent as a normalized double. This routine does 'ldexp' - * in two steps for those numbers, to return correctly denormalized - * values. - * - *---------------------------------------------------------------------- - */ - -static double -SafeLdExp( - double fract, - int expt) -{ - int minexpt = DBL_MIN_EXP * log2FLT_RADIX; - volatile double a, b, retval; - - if (expt < minexpt) { - a = ldexp(fract, expt - mantBits - minexpt); - b = ldexp(1.0, mantBits + minexpt); - retval = a * b; - } else { - retval = ldexp(fract, expt); - } - return retval; -} - -/* - *---------------------------------------------------------------------- - * - * TclFormatNaN -- - * - * Makes the string representation of a "Not a Number" - * - * Results: - * None. - * - * Side effects: - * Stores the string representation in the supplied buffer, which must be - * at least TCL_DOUBLE_SPACE characters. - * - *---------------------------------------------------------------------- - */ - -void -TclFormatNaN( - double value, /* The Not-a-Number to format. */ - char *buffer) /* String representation. */ -{ -#ifndef IEEE_FLOATING_POINT - strcpy(buffer, "NaN"); - return; -#else - union { - double dv; - Tcl_WideUInt iv; - } bitwhack; - - bitwhack.dv = value; - if (n770_fp) { - bitwhack.iv = Nokia770Twiddle(bitwhack.iv); - } - if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { - bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63); - *buffer++ = '-'; - } - *buffer++ = 'N'; - *buffer++ = 'a'; - *buffer++ = 'N'; - bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1; - if (bitwhack.iv != 0) { - sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv); - } else { - *buffer = '\0'; - } -#endif /* IEEE_FLOATING_POINT */ -} - -/* - *---------------------------------------------------------------------- - * - * Nokia770Twiddle -- - * - * Transpose the two words of a number for Nokia 770 floating point - * handling. - * - *---------------------------------------------------------------------- - */ -#ifdef IEEE_FLOATING_POINT -static Tcl_WideUInt -Nokia770Twiddle( - Tcl_WideUInt w) /* Number to transpose. */ -{ - return (((w >> 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: - */ |