summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
Diffstat (limited to 'generic/tclStrToD.c')
-rwxr-xr-xgeneric/tclStrToD.c4962
1 files changed, 4962 insertions, 0 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
new file mode 100755
index 0000000..d5c6e9c
--- /dev/null
+++ b/generic/tclStrToD.c
@@ -0,0 +1,4962 @@
+/*
+ *----------------------------------------------------------------------
+ *
+ * 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
+
+/*
+ * gcc on x86 needs access to rounding controls, because of a questionable
+ * feature where it retains intermediate results as IEEE 'long double' values
+ * somewhat unpredictably. It is tempting to include fpu_control.h, but that
+ * file exists only on Linux; it is missing on Cygwin and MinGW. Most gcc-isms
+ * and ix86-isms are factored out here.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
+#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
+#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
+# define FPU_IEEE_ROUNDING 0x027f
+# define ADJUST_FPU_CONTROL_WORD
+#endif
+
+/* Sun ProC needs sunmath for rounding control on x86 like gcc above.
+ *
+ *
+ */
+#if defined(__sun) && defined(__i386) && !defined(__GNUC__)
+#include <sunmath.h>
+#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 Tcl_WideUInt wtens[] = {
+ 1, 10, 100, 1000, 10000, 100000, 1000000,
+ (Tcl_WideUInt) 1000000*10, (Tcl_WideUInt) 1000000*100,
+ (Tcl_WideUInt) 1000000*1000, (Tcl_WideUInt) 1000000*10000,
+ (Tcl_WideUInt) 1000000*100000, (Tcl_WideUInt) 1000000*1000000,
+ (Tcl_WideUInt) 1000000*1000000*10, (Tcl_WideUInt) 1000000*1000000*100,
+ (Tcl_WideUInt) 1000000*1000000*1000,(Tcl_WideUInt) 1000000*1000000*10000
+
+};
+
+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);
+static double MakeNaN(int signum, Tcl_WideUInt tag);
+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(mp_int *big, int *machexp);
+static double Pow10TimesFrExp(int exponent, double fraction,
+ int *machexp);
+static double SafeLdExp(double fraction, int exponent);
+static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclParseNumber --
+ *
+ * Scans bytes, interpreted as characters in Tcl's internal encoding, and
+ * parses the longest prefix that is the string representation of a
+ * number in a format recognized by Tcl.
+ *
+ * The arguments bytes, numBytes, and objPtr are the inputs which
+ * determine the string to be parsed. If bytes is non-NULL, it points to
+ * the first byte to be scanned. If bytes is NULL, then objPtr must be
+ * non-NULL, and the string representation of objPtr will be scanned
+ * (generated first, if necessary). The numBytes argument determines the
+ * number of bytes to be scanned. If numBytes is negative, the first NUL
+ * byte encountered will terminate the scan. If numBytes is non-negative,
+ * then no more than numBytes bytes will be scanned.
+ *
+ * The argument flags is an input that controls the numeric formats
+ * recognized by the parser. The flag bits are:
+ *
+ * - TCL_PARSE_INTEGER_ONLY: accept only integer values; reject
+ * strings that denote floating point values (or accept only the
+ * leading portion of them that are integer values).
+ * - TCL_PARSE_SCAN_PREFIXES: ignore the prefixes 0b and 0o that are
+ * not part of the [scan] command's vocabulary. Use only in
+ * combination with TCL_PARSE_INTEGER_ONLY.
+ * - TCL_PARSE_OCTAL_ONLY: parse only in the octal format, whether
+ * or not a prefix is present that would lead to octal parsing.
+ * Use only in combination with TCL_PARSE_INTEGER_ONLY.
+ * - TCL_PARSE_HEXADECIMAL_ONLY: parse only in the hexadecimal format,
+ * whether or not a prefix is present that would lead to
+ * hexadecimal parsing. Use only in combination with
+ * TCL_PARSE_INTEGER_ONLY.
+ * - TCL_PARSE_DECIMAL_ONLY: parse only in the decimal format, no
+ * matter whether a 0 prefix would normally force a different
+ * base.
+ * - TCL_PARSE_NO_WHITESPACE: reject any leading/trailing whitespace
+ *
+ * The arguments interp and expected are inputs that control error
+ * message generation. If interp is NULL, no error message will be
+ * generated. If interp is non-NULL, then expected must also be non-NULL.
+ * When TCL_ERROR is returned, an error message will be left in the
+ * result of interp, and the expected argument will appear in the error
+ * message as the thing TclParseNumber expected, but failed to find in
+ * the string.
+ *
+ * The arguments objPtr and endPtrPtr as well as the return code are the
+ * outputs.
+ *
+ * When the parser cannot find any prefix of the string that matches a
+ * format it is looking for, TCL_ERROR is returned and an error message
+ * may be generated and returned as described above. The contents of
+ * objPtr will not be changed. If endPtrPtr is non-NULL, a pointer to the
+ * character in the string that terminated the scan will be written to
+ * *endPtrPtr.
+ *
+ * When the parser determines that the entire string matches a format it
+ * is looking for, TCL_OK is returned, and if objPtr is non-NULL, then
+ * the internal rep and Tcl_ObjType of objPtr are set to the "canonical"
+ * numeric value that matches the scanned string. If endPtrPtr is not
+ * NULL, a pointer to the end of the string will be written to *endPtrPtr
+ * (that is, either bytes+numBytes or a pointer to a terminating NUL
+ * byte).
+ *
+ * When the parser determines that a partial string matches a format it
+ * is looking for, the value of endPtrPtr determines what happens:
+ *
+ * - If endPtrPtr is NULL, then TCL_ERROR is returned, with error message
+ * generation as above.
+ *
+ * - If endPtrPtr is non-NULL, then TCL_OK is returned and objPtr
+ * internals are set as above. Also, a pointer to the first
+ * character following the parsed numeric string is written to
+ * *endPtrPtr.
+ *
+ * In some cases where the string being scanned is the string rep of
+ * objPtr, this routine can leave objPtr in an inconsistent state where
+ * its string rep and its internal rep do not agree. In these cases the
+ * internal rep will be in agreement with only some substring of the
+ * string rep. This might happen if the caller passes in a non-NULL bytes
+ * value that points somewhere into the string rep. It might happen if
+ * the caller passes in a numBytes value that limits the scan to only a
+ * prefix of the string rep. Or it might happen if a non-NULL value of
+ * endPtrPtr permits a TCL_OK return from only a partial string match. It
+ * is the responsibility of the caller to detect and correct such
+ * inconsistencies when they can and do arise.
+ *
+ * Results:
+ * Returns a standard Tcl result.
+ *
+ * Side effects:
+ * The string representaton of objPtr may be generated.
+ *
+ * The internal representation and Tcl_ObjType of objPtr may be changed.
+ * This may involve allocation and/or freeing of memory.
+ *
+ *----------------------------------------------------------------------
+ */
+
+int
+TclParseNumber(
+ Tcl_Interp *interp, /* Used for error reporting. May be NULL. */
+ Tcl_Obj *objPtr, /* Object to receive the internal rep. */
+ const char *expected, /* Description of the type of number the
+ * caller expects to be able to parse
+ * ("integer", "boolean value", etc.). */
+ const char *bytes, /* Pointer to the start of the string to
+ * scan. */
+ int numBytes, /* Maximum number of bytes to scan, see
+ * above. */
+ const char **endPtrPtr, /* Place to store pointer to the character
+ * that terminated the scan. */
+ int flags) /* Flags governing the parse. */
+{
+ enum State {
+ INITIAL, SIGNUM, ZERO, ZERO_X,
+ ZERO_O, ZERO_B, BINARY,
+ HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL,
+ LEADING_RADIX_POINT, FRACTION,
+ EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
+ sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
+#ifdef IEEE_FLOATING_POINT
+ , sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
+#endif
+ } state = INITIAL;
+ enum State acceptState = INITIAL;
+
+ int signum = 0; /* Sign of the number being parsed */
+ Tcl_WideUInt significandWide = 0;
+ /* Significand of the number being parsed (if
+ * no overflow) */
+ mp_int significandBig; /* Significand of the number being parsed (if
+ * it overflows significandWide) */
+ int significandOverflow = 0;/* Flag==1 iff significandBig is used */
+ Tcl_WideUInt octalSignificandWide = 0;
+ /* Significand of an octal number; needed
+ * because we don't know whether a number with
+ * a leading zero is octal or decimal until
+ * we've scanned forward to a '.' or 'e' */
+ mp_int octalSignificandBig; /* Significand of octal number once
+ * octalSignificandWide overflows */
+ int octalSignificandOverflow = 0;
+ /* Flag==1 if octalSignificandBig is used */
+ int numSigDigs = 0; /* Number of significant digits in the decimal
+ * significand */
+ int numTrailZeros = 0; /* Number of trailing zeroes at the current
+ * point in the parse. */
+ int numDigitsAfterDp = 0; /* Number of digits scanned after the decimal
+ * point */
+ int exponentSignum = 0; /* Signum of the exponent of a floating point
+ * number */
+ long exponent = 0; /* Exponent of a floating point number */
+ const char *p; /* Pointer to next character to scan */
+ size_t len; /* Number of characters remaining after p */
+ const char *acceptPoint; /* Pointer to position after last character in
+ * an acceptable number */
+ size_t acceptLen; /* Number of characters following that
+ * point. */
+ int status = TCL_OK; /* Status to return to caller */
+ char d = 0; /* Last hexadecimal digit scanned; initialized
+ * to avoid a compiler warning. */
+ int shift = 0; /* Amount to shift when accumulating binary */
+ int explicitOctal = 0;
+
+#define ALL_BITS (~(Tcl_WideUInt)0)
+#define MOST_BITS (ALL_BITS >> 1)
+
+ /*
+ * Initialize bytes to start of the object's string rep if the caller
+ * didn't pass anything else.
+ */
+
+ if (bytes == NULL) {
+ bytes = TclGetString(objPtr);
+ }
+
+ p = bytes;
+ len = numBytes;
+ acceptPoint = p;
+ acceptLen = len;
+ while (1) {
+ char c = len ? *p : '\0';
+ switch (state) {
+
+ case INITIAL:
+ /*
+ * Initial state. Acceptable characters are +, -, digits, period,
+ * I, N, and whitespace.
+ */
+
+ if (isspace(UCHAR(c))) {
+ if (flags & TCL_PARSE_NO_WHITESPACE) {
+ goto endgame;
+ }
+ break;
+ } else if (c == '+') {
+ state = SIGNUM;
+ break;
+ } else if (c == '-') {
+ signum = 1;
+ state = SIGNUM;
+ break;
+ }
+ /* FALLTHROUGH */
+
+ case SIGNUM:
+ /*
+ * Scanned a leading + or -. Acceptable characters are digits,
+ * period, I, and N.
+ */
+
+ if (c == '0') {
+ if (flags & TCL_PARSE_DECIMAL_ONLY) {
+ state = DECIMAL;
+ } else {
+ state = ZERO;
+ }
+ break;
+ } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
+ goto zerox;
+ } else if (flags & TCL_PARSE_OCTAL_ONLY) {
+ goto zeroo;
+ } else if (isdigit(UCHAR(c))) {
+ significandWide = c - '0';
+ numSigDigs = 1;
+ state = DECIMAL;
+ break;
+ } else if (flags & TCL_PARSE_INTEGER_ONLY) {
+ goto endgame;
+ } else if (c == '.') {
+ state = LEADING_RADIX_POINT;
+ break;
+ } else if (c == 'I' || c == 'i') {
+ state = sI;
+ break;
+#ifdef IEEE_FLOATING_POINT
+ } else if (c == 'N' || c == 'n') {
+ state = sN;
+ break;
+#endif
+ }
+ goto endgame;
+
+ case ZERO:
+ /*
+ * Scanned a leading zero (perhaps with a + or -). Acceptable
+ * inputs are digits, period, X, 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 (c == 'o' || c == 'O') {
+ explicitOctal = 1;
+ state = ZERO_O;
+ break;
+ }
+#ifdef KILL_OCTAL
+ goto decimal;
+#endif
+ /* FALLTHROUGH */
+
+ case OCTAL:
+ /*
+ * Scanned an optional + or -, followed by a string of octal
+ * digits. Acceptable inputs are more digits, period, or E. If 8
+ * or 9 is encountered, commit to floating point.
+ */
+
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ /* FALLTHROUGH */
+ case ZERO_O:
+ zeroo:
+ if (c == '0') {
+ numTrailZeros++;
+ state = OCTAL;
+ break;
+ } else if (c >= '1' && c <= '7') {
+ if (objPtr != NULL) {
+ shift = 3 * (numTrailZeros + 1);
+ significandOverflow = AccumulateDecimalDigit(
+ (unsigned)(c-'0'), numTrailZeros,
+ &significandWide, &significandBig,
+ significandOverflow);
+
+ if (!octalSignificandOverflow) {
+ /*
+ * Shifting by more bits than are in the value being
+ * shifted is at least de facto nonportable. Check for
+ * too large shifts first.
+ */
+
+ if ((octalSignificandWide != 0)
+ && (((size_t)shift >=
+ CHAR_BIT*sizeof(Tcl_WideUInt))
+ || (octalSignificandWide >
+ (~(Tcl_WideUInt)0 >> shift)))) {
+ octalSignificandOverflow = 1;
+ TclBNInitBignumFromWideUInt(&octalSignificandBig,
+ octalSignificandWide);
+ }
+ }
+ if (!octalSignificandOverflow) {
+ octalSignificandWide =
+ (octalSignificandWide << shift) + (c - '0');
+ } else {
+ mp_mul_2d(&octalSignificandBig, shift,
+ &octalSignificandBig);
+ mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
+ &octalSignificandBig);
+ }
+ }
+ if (numSigDigs != 0) {
+ numSigDigs += numTrailZeros+1;
+ } else {
+ numSigDigs = 1;
+ }
+ numTrailZeros = 0;
+ state = OCTAL;
+ break;
+ }
+ /* FALLTHROUGH */
+
+ case BAD_OCTAL:
+ if (explicitOctal) {
+ /*
+ * No forgiveness for bad digits in explicitly octal numbers.
+ */
+
+ goto endgame;
+ }
+ if (flags & TCL_PARSE_INTEGER_ONLY) {
+ /*
+ * No seeking floating point when parsing only integer.
+ */
+
+ goto endgame;
+ }
+#ifndef KILL_OCTAL
+
+ /*
+ * Scanned a number with a leading zero that contains an 8, 9,
+ * radix point or E. This is an invalid octal number, but might
+ * still be floating point.
+ */
+
+ if (c == '0') {
+ numTrailZeros++;
+ state = BAD_OCTAL;
+ break;
+ } else if (isdigit(UCHAR(c))) {
+ if (objPtr != NULL) {
+ significandOverflow = AccumulateDecimalDigit(
+ (unsigned)(c-'0'), numTrailZeros,
+ &significandWide, &significandBig,
+ significandOverflow);
+ }
+ if (numSigDigs != 0) {
+ numSigDigs += (numTrailZeros + 1);
+ } else {
+ numSigDigs = 1;
+ }
+ numTrailZeros = 0;
+ state = BAD_OCTAL;
+ break;
+ } else if (c == '.') {
+ state = FRACTION;
+ break;
+ } else if (c == 'E' || c == 'e') {
+ state = EXPONENT_START;
+ break;
+ }
+#endif
+ goto endgame;
+
+ /*
+ * Scanned 0x. If state is HEXADECIMAL, scanned at least one
+ * character following the 0x. The only acceptable inputs are
+ * hexadecimal digits.
+ */
+
+ case HEXADECIMAL:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ /* FALLTHROUGH */
+
+ case ZERO_X:
+ zerox:
+ if (c == '0') {
+ numTrailZeros++;
+ state = HEXADECIMAL;
+ break;
+ } else if (isdigit(UCHAR(c))) {
+ d = (c-'0');
+ } else if (c >= 'A' && c <= 'F') {
+ d = (c-'A'+10);
+ } else if (c >= 'a' && c <= 'f') {
+ d = (c-'a'+10);
+ } else {
+ goto endgame;
+ }
+ if (objPtr != NULL) {
+ shift = 4 * (numTrailZeros + 1);
+ if (!significandOverflow) {
+ /*
+ * Shifting by more bits than are in the value being
+ * shifted is at least de facto nonportable. Check for too
+ * large shifts first.
+ */
+
+ if (significandWide != 0 &&
+ ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
+ significandWide > (~(Tcl_WideUInt)0 >> shift))) {
+ significandOverflow = 1;
+ TclBNInitBignumFromWideUInt(&significandBig,
+ significandWide);
+ }
+ }
+ if (!significandOverflow) {
+ significandWide = (significandWide << shift) + d;
+ } else {
+ mp_mul_2d(&significandBig, shift, &significandBig);
+ mp_add_d(&significandBig, (mp_digit) d, &significandBig);
+ }
+ }
+ numTrailZeros = 0;
+ state = HEXADECIMAL;
+ break;
+
+ case BINARY:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ case ZERO_B:
+ if (c == '0') {
+ numTrailZeros++;
+ state = BINARY;
+ break;
+ } else if (c != '1') {
+ goto endgame;
+ }
+ if (objPtr != NULL) {
+ shift = numTrailZeros + 1;
+ if (!significandOverflow) {
+ /*
+ * Shifting by more bits than are in the value being
+ * shifted is at least de facto nonportable. Check for too
+ * large shifts first.
+ */
+
+ if (significandWide != 0 &&
+ ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
+ significandWide > (~(Tcl_WideUInt)0 >> shift))) {
+ significandOverflow = 1;
+ TclBNInitBignumFromWideUInt(&significandBig,
+ significandWide);
+ }
+ }
+ if (!significandOverflow) {
+ significandWide = (significandWide << shift) + 1;
+ } else {
+ mp_mul_2d(&significandBig, shift, &significandBig);
+ mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
+ }
+ }
+ numTrailZeros = 0;
+ state = BINARY;
+ break;
+
+ case DECIMAL:
+ /*
+ * Scanned an optional + or - followed by a string of decimal
+ * digits.
+ */
+
+#ifdef KILL_OCTAL
+ decimal:
+#endif
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ if (c == '0') {
+ numTrailZeros++;
+ state = DECIMAL;
+ break;
+ } else if (isdigit(UCHAR(c))) {
+ if (objPtr != NULL) {
+ significandOverflow = AccumulateDecimalDigit(
+ (unsigned)(c - '0'), numTrailZeros,
+ &significandWide, &significandBig,
+ significandOverflow);
+ }
+ numSigDigs += numTrailZeros+1;
+ numTrailZeros = 0;
+ state = DECIMAL;
+ break;
+ } else if (flags & TCL_PARSE_INTEGER_ONLY) {
+ goto endgame;
+ } else if (c == '.') {
+ state = FRACTION;
+ break;
+ } else if (c == 'E' || c == 'e') {
+ state = EXPONENT_START;
+ break;
+ }
+ goto endgame;
+
+ /*
+ * Found a decimal point. If no digits have yet been scanned, E is
+ * not allowed; otherwise, it introduces the exponent. If at least
+ * one digit has been found, we have a possible complete number.
+ */
+
+ case FRACTION:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ if (c == 'E' || c=='e') {
+ state = EXPONENT_START;
+ break;
+ }
+ /* FALLTHROUGH */
+
+ case LEADING_RADIX_POINT:
+ if (c == '0') {
+ numDigitsAfterDp++;
+ numTrailZeros++;
+ state = FRACTION;
+ break;
+ } else if (isdigit(UCHAR(c))) {
+ numDigitsAfterDp++;
+ if (objPtr != NULL) {
+ significandOverflow = AccumulateDecimalDigit(
+ (unsigned)(c-'0'), numTrailZeros,
+ &significandWide, &significandBig,
+ significandOverflow);
+ }
+ if (numSigDigs != 0) {
+ numSigDigs += numTrailZeros+1;
+ } else {
+ numSigDigs = 1;
+ }
+ numTrailZeros = 0;
+ state = FRACTION;
+ break;
+ }
+ goto endgame;
+
+ case EXPONENT_START:
+ /*
+ * Scanned the E at the start of an exponent. Make sure a legal
+ * character follows before using the C library strtol routine,
+ * which allows whitespace.
+ */
+
+ if (c == '+') {
+ state = EXPONENT_SIGNUM;
+ break;
+ } else if (c == '-') {
+ exponentSignum = 1;
+ state = EXPONENT_SIGNUM;
+ break;
+ }
+ /* FALLTHROUGH */
+
+ case EXPONENT_SIGNUM:
+ /*
+ * Found the E at the start of the exponent, followed by a sign
+ * character.
+ */
+
+ if (isdigit(UCHAR(c))) {
+ exponent = c - '0';
+ state = EXPONENT;
+ break;
+ }
+ goto endgame;
+
+ case EXPONENT:
+ /*
+ * Found an exponent with at least one digit. Accumulate it,
+ * making sure to hard-pin it to LONG_MAX on overflow.
+ */
+
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ if (isdigit(UCHAR(c))) {
+ if (exponent < (LONG_MAX - 9) / 10) {
+ exponent = 10 * exponent + (c - '0');
+ } else {
+ exponent = LONG_MAX;
+ }
+ state = EXPONENT;
+ break;
+ }
+ goto endgame;
+
+ /*
+ * Parse out INFINITY by simply spelling it out. INF is accepted
+ * as an abbreviation; other prefices are not.
+ */
+
+ case sI:
+ if (c == 'n' || c == 'N') {
+ state = sIN;
+ break;
+ }
+ goto endgame;
+ case sIN:
+ if (c == 'f' || c == 'F') {
+ state = sINF;
+ break;
+ }
+ goto endgame;
+ case sINF:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ if (c == 'i' || c == 'I') {
+ state = sINFI;
+ break;
+ }
+ goto endgame;
+ case sINFI:
+ if (c == 'n' || c == 'N') {
+ state = sINFIN;
+ break;
+ }
+ goto endgame;
+ case sINFIN:
+ if (c == 'i' || c == 'I') {
+ state = sINFINI;
+ break;
+ }
+ goto endgame;
+ case sINFINI:
+ if (c == 't' || c == 'T') {
+ state = sINFINIT;
+ break;
+ }
+ goto endgame;
+ case sINFINIT:
+ if (c == 'y' || c == 'Y') {
+ state = sINFINITY;
+ break;
+ }
+ goto endgame;
+
+ /*
+ * Parse NaN's.
+ */
+#ifdef IEEE_FLOATING_POINT
+ case sN:
+ if (c == 'a' || c == 'A') {
+ state = sNA;
+ break;
+ }
+ goto endgame;
+ case sNA:
+ if (c == 'n' || c == 'N') {
+ state = sNAN;
+ break;
+ }
+ goto endgame;
+ case sNAN:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ if (c == '(') {
+ state = sNANPAREN;
+ break;
+ }
+ goto endgame;
+
+ /*
+ * Parse NaN(hexdigits)
+ */
+ case sNANHEX:
+ if (c == ')') {
+ state = sNANFINISH;
+ break;
+ }
+ /* FALLTHROUGH */
+ case sNANPAREN:
+ if (isspace(UCHAR(c))) {
+ break;
+ }
+ if (numSigDigs < 13) {
+ if (c >= '0' && c <= '9') {
+ d = c - '0';
+ } else if (c >= 'a' && c <= 'f') {
+ d = 10 + c - 'a';
+ } else if (c >= 'A' && c <= 'F') {
+ d = 10 + c - 'A';
+ }
+ significandWide = (significandWide << 4) + d;
+ state = sNANHEX;
+ break;
+ }
+ goto endgame;
+ case sNANFINISH:
+#endif
+
+ case sINFINITY:
+ acceptState = state;
+ acceptPoint = p;
+ acceptLen = len;
+ goto endgame;
+ }
+ p++;
+ len--;
+ }
+
+ endgame:
+ if (acceptState == INITIAL) {
+ /*
+ * No numeric string at all found.
+ */
+
+ status = TCL_ERROR;
+ if (endPtrPtr != NULL) {
+ *endPtrPtr = p;
+ }
+ } else {
+ /*
+ * Back up to the last accepting state in the lexer.
+ */
+
+ p = acceptPoint;
+ len = acceptLen;
+ if (!(flags & TCL_PARSE_NO_WHITESPACE)) {
+ /*
+ * Accept trailing whitespace.
+ */
+
+ while (len != 0 && isspace(UCHAR(*p))) {
+ p++;
+ len--;
+ }
+ }
+ if (endPtrPtr == NULL) {
+ if ((len != 0) && ((numBytes > 0) || (*p != '\0'))) {
+ status = TCL_ERROR;
+ }
+ } else {
+ *endPtrPtr = p;
+ }
+ }
+
+ /*
+ * Generate and store the appropriate internal rep.
+ */
+
+ if (status == TCL_OK && objPtr != NULL) {
+ TclFreeIntRep(objPtr);
+ switch (acceptState) {
+ case SIGNUM:
+ case BAD_OCTAL:
+ case ZERO_X:
+ case ZERO_O:
+ case ZERO_B:
+ case LEADING_RADIX_POINT:
+ case EXPONENT_START:
+ case EXPONENT_SIGNUM:
+ case sI:
+ case sIN:
+ case sINFI:
+ case sINFIN:
+ case sINFINI:
+ case sINFINIT:
+#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;
+
+ TclNewLiteralStringObj(msg, "expected ");
+ Tcl_AppendToObj(msg, expected, -1);
+ Tcl_AppendToObj(msg, " but got \"", -1);
+ Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, "");
+ Tcl_AppendToObj(msg, "\"", -1);
+ if (state == BAD_OCTAL) {
+ Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1);
+ }
+ Tcl_SetObjResult(interp, msg);
+ 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.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+ fpu_control_t roundTo53Bits = 0x027f;
+ fpu_control_t oldRoundingMode;
+ _FPU_GETCW(oldRoundingMode);
+ _FPU_SETCW(roundTo53Bits);
+#endif
+#if defined(__sun) && defined(__i386) && !defined(__GNUC__)
+ ieee_flags("set","precision","double",NULL);
+#endif
+
+ /*
+ * Test for the easy cases.
+ */
+
+ if (numSigDigs <= DBL_DIG) {
+ if (exponent >= 0) {
+ if (exponent <= mmaxpow) {
+ /*
+ * The significand is an exact integer, and so is
+ * 10**exponent. The product will be correct to within 1/2 ulp
+ * without special handling.
+ */
+
+ retval = (double)(Tcl_WideInt)significand * 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.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+ _FPU_SETCW(oldRoundingMode);
+#endif
+#if defined(__sun) && defined(__i386) && !defined(__GNUC__)
+ ieee_flags("clear","precision",NULL,NULL);
+#endif
+
+ return retval;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * MakeHighPrecisionDouble --
+ *
+ * Makes the double precision number, signum*significand*10**exponent.
+ *
+ * Results:
+ * Returns the constructed number.
+ *
+ * MakeHighPrecisionDouble is used when arbitrary-precision arithmetic is
+ * needed to ensure correct rounding. It begins by calculating a
+ * low-precision approximation to the desired number, and then refines
+ * the answer in high precision.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+MakeHighPrecisionDouble(
+ int signum, /* 1=negative, 0=nonnegative */
+ mp_int *significand, /* Exact significand of the number */
+ int numSigDigs, /* Number of significant digits */
+ int exponent) /* Power of 10 by which to multiply */
+{
+ double retval;
+ int machexp; /* Machine exponent of a power of 10 */
+
+ /*
+ * With gcc on x86, the floating point rounding mode is double-extended.
+ * This causes the result of double-precision calculations to be rounded
+ * twice: once to the precision of double-extended and then again to the
+ * precision of double. Double-rounding introduces gratuitous errors of 1
+ * ulp, so we need to change rounding mode to 53-bits.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+ fpu_control_t roundTo53Bits = 0x027f;
+ fpu_control_t oldRoundingMode;
+ _FPU_GETCW(oldRoundingMode);
+ _FPU_SETCW(roundTo53Bits);
+#endif
+#if defined(__sun) && defined(__i386) && !defined(__GNUC__)
+ ieee_flags("set","precision","double",NULL);
+#endif
+
+ /*
+ * Quick checks for over/underflow.
+ */
+
+ if (numSigDigs+exponent-1 > maxDigits) {
+ retval = HUGE_VAL;
+ goto returnValue;
+ }
+ if (numSigDigs+exponent-1 < minDigits) {
+ retval = 0;
+ goto returnValue;
+ }
+
+ /*
+ * Develop a first approximation to the significand. It is tempting simply
+ * to force bignum to double, but that will overflow on input numbers like
+ * 1.[string repeat 0 1000]1; while this is a not terribly likely
+ * scenario, we still have to deal with it. Use fraction and exponent
+ * instead. Once we have the significand, multiply by 10**exponent. Test
+ * for overflow. Convert back to a double, and test for underflow.
+ */
+
+ retval = BignumToBiasedFrExp(significand, &machexp);
+ retval = Pow10TimesFrExp(exponent, retval, &machexp);
+ if (machexp > DBL_MAX_EXP*log2FLT_RADIX) {
+ retval = HUGE_VAL;
+ goto returnValue;
+ }
+ retval = SafeLdExp(retval, machexp);
+ if (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.
+ */
+
+#if defined(__GNUC__) && defined(__i386)
+ _FPU_SETCW(oldRoundingMode);
+#endif
+#if defined(__sun) && defined(__i386) && !defined(__GNUC__)
+ ieee_flags("clear","precision",NULL,NULL);
+#endif
+ return retval;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * MakeNaN --
+ *
+ * Makes a "Not a Number" given a set of bits to put in the tag bits
+ *
+ * Note that a signalling NaN is never returned.
+ *
+ *----------------------------------------------------------------------
+ */
+
+#ifdef IEEE_FLOATING_POINT
+static double
+MakeNaN(
+ int signum, /* Sign bit (1=negative, 0=nonnegative */
+ Tcl_WideUInt tags) /* Tag bits to put in the NaN */
+{
+ union {
+ Tcl_WideUInt iv;
+ double dv;
+ } theNaN;
+
+ theNaN.iv = tags;
+ theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
+ if (signum) {
+ theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48;
+ } else {
+ theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
+ }
+ if (n770_fp) {
+ theNaN.iv = Nokia770Twiddle(theNaN.iv);
+ }
+ return theNaN.dv;
+}
+#endif
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * RefineApproximation --
+ *
+ * Given a poor approximation to a floating point number, returns a
+ * better one. (The better approximation is correct to within 1 ulp, and
+ * is entirely correct if the poor approximation is correct to 1 ulp.)
+ *
+ * Results:
+ * Returns the improved result.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+RefineApproximation(
+ double approxResult, /* Approximate result of conversion */
+ mp_int *exactSignificand, /* Integer significand */
+ int exponent) /* Power of 10 to multiply by significand */
+{
+ int M2, M5; /* Powers of 2 and of 5 needed to put the
+ * decimal and binary numbers over a common
+ * denominator. */
+ double significand; /* Sigificand of the binary number */
+ int binExponent; /* Exponent of the binary number */
+ int msb; /* Most significant bit position of an
+ * intermediate result */
+ int nDigits; /* Number of mp_digit's in an intermediate
+ * result */
+ mp_int twoMv; /* Approx binary value expressed as an exact
+ * integer scaled by the multiplier 2M */
+ mp_int twoMd; /* Exact decimal value expressed as an exact
+ * integer scaled by the multiplier 2M */
+ int scale; /* Scale factor for M */
+ int multiplier; /* Power of two to scale M */
+ double num, den; /* Numerator and denominator of the correction
+ * term */
+ double quot; /* Correction term */
+ double minincr; /* Lower bound on the absolute value of the
+ * correction term. */
+ int i;
+
+ /*
+ * The first approximation is always low. If we find that it's HUGE_VAL,
+ * we're done.
+ */
+
+ if (approxResult == HUGE_VAL) {
+ return approxResult;
+ }
+
+ /*
+ * Find a common denominator for the decimal and binary fractions. The
+ * common denominator will be 2**M2 + 5**M5.
+ */
+
+ significand = frexp(approxResult, &binExponent);
+ i = mantBits - binExponent;
+ if (i < 0) {
+ M2 = 0;
+ } else {
+ M2 = i;
+ }
+ if (exponent > 0) {
+ M5 = 0;
+ } else {
+ M5 = -exponent;
+ if ((M5-1) > M2) {
+ M2 = M5-1;
+ }
+ }
+
+ /*
+ * The floating point number is significand*2**binExponent. Compute the
+ * large integer significand*2**(binExponent+M2+1). The 2**-1 bit of the
+ * significand (the most significant) corresponds to the
+ * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold
+ * that quantity, then convert the significand to a large integer, scaled
+ * appropriately. Then multiply by the appropriate power of 5.
+ */
+
+ msb = binExponent + M2; /* 1008 */
+ nDigits = msb / DIGIT_BIT + 1;
+ mp_init_size(&twoMv, nDigits);
+ i = (msb % DIGIT_BIT + 1);
+ twoMv.used = nDigits;
+ significand *= SafeLdExp(1.0, i);
+ while (--nDigits >= 0) {
+ twoMv.dp[nDigits] = (mp_digit) significand;
+ significand -= (mp_digit) significand;
+ significand = SafeLdExp(significand, DIGIT_BIT);
+ }
+ for (i = 0; i <= 8; ++i) {
+ if (M5 & (1 << i)) {
+ mp_mul(&twoMv, pow5+i, &twoMv);
+ }
+ }
+
+ /*
+ * Collect the decimal significand as a high precision integer. The least
+ * significant bit corresponds to bit M2+exponent+1 so it will need to be
+ * shifted left by that many bits after being multiplied by
+ * 5**(M5+exponent).
+ */
+
+ mp_init_copy(&twoMd, exactSignificand);
+ for (i=0; i<=8; ++i) {
+ if ((M5+exponent) & (1 << i)) {
+ mp_mul(&twoMd, pow5+i, &twoMd);
+ }
+ }
+ mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
+ mp_sub(&twoMd, &twoMv, &twoMd);
+
+ /*
+ * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
+ * term. Because 2M may well overflow a double, we need to scale the
+ * denominator by a factor of 2**binExponent-mantBits
+ */
+
+ scale = binExponent - mantBits - 1;
+
+ mp_set(&twoMv, 1);
+ for (i=0; i<=8; ++i) {
+ if (M5 & (1 << i)) {
+ mp_mul(&twoMv, pow5+i, &twoMv);
+ }
+ }
+ multiplier = M2 + scale + 1;
+ if (multiplier > 0) {
+ mp_mul_2d(&twoMv, multiplier, &twoMv);
+ } else if (multiplier < 0) {
+ mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
+ }
+
+ /*
+ * If the result is less than unity, the error is less than 1/2 unit in
+ * the last place, so there's no correction to make.
+ */
+
+ if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) {
+ mp_clear(&twoMd);
+ mp_clear(&twoMv);
+ return approxResult;
+ }
+
+ /*
+ * Convert the numerator and denominator of the corrector term accurately
+ * to floating point numbers.
+ */
+
+ num = TclBignumToDouble(&twoMd);
+ den = TclBignumToDouble(&twoMv);
+
+ quot = SafeLdExp(num/den, scale);
+ minincr = SafeLdExp(1.0, binExponent-mantBits);
+
+ if (quot<0. && quot>-minincr) {
+ quot = -minincr;
+ } else if (quot>0. && quot<minincr) {
+ quot = minincr;
+ }
+
+ mp_clear(&twoMd);
+ mp_clear(&twoMv);
+
+ return approxResult + quot;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * 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;
+}
+
+/*
+ *-----------------------------------------------------------------------------0
+ *
+ * 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'.
+ *
+ * Side effects:
+ * Stores 9999 in *decpt, and sets '*endPtr' to designate the
+ * terminating NUL byte of the string if 'endPtr' is not NULL.
+ *
+ * The string returned must be freed by the caller using 'ckfree'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+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 d, /* 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
+ * musguessed 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 */
+
+ /*
+ * Bring d into the range [1 .. 10)
+ */
+ ieps = AdjustRange(&d, k);
+
+ /*
+ * 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 denominatorr */
+ 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;
+ const static 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);
+ } else {
+ 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, 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 groun */
+ int i, j;
+
+ /*
+ * 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);
+ ilim =ilim1;
+ --k;
+ }
+ mp_init(&temp);
+
+ /* Convert the leading digit */
+
+ mp_init(&dig);
+ 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);
+ } else {
+ while (*--s == '0') {
+ /* do nothing */
+ }
+ ++s;
+ }
+ break;
+ }
+ }
+ }
+ /*
+ * 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;
+
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * 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 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
+ * 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)) {
+ if ((retval = QuickConversion(d.d, k, k_check, flags,
+ i, ilim, ilim1,
+ decpt, endPtr)) != 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 = (Tcl_WideUInt *)
+ ckalloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
+ u = 1;
+ for (i = 0; i < maxpow10_wide; ++i) {
+ pow10_wide[i] = u;
+ u *= 10;
+ }
+ pow10_wide[i] = u;
+
+ /*
+ * Determine how many bits of precision a double has, and how many
+ * decimal digits that represents.
+ */
+
+ if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
+ Tcl_Panic("This code doesn't work on a decimal machine!");
+ }
+ log2FLT_RADIX--;
+ mantBits = DBL_MANT_DIG * log2FLT_RADIX;
+ d = 1.0;
+
+ /*
+ * Initialize a table of powers of ten that can be exactly represented
+ * in a double.
+ */
+
+ x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
+ if (x < MAXPOW) {
+ mmaxpow = x;
+ } else {
+ mmaxpow = MAXPOW;
+ }
+ for (i=0 ; i<=mmaxpow ; ++i) {
+ 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((char *) pow10_wide);
+ for (i=0; i<9; ++i) {
+ mp_clear(pow5 + i);
+ }
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * Tcl_InitBignumFromDouble --
+ *
+ * Extracts the integer part of a double and converts it to an arbitrary
+ * precision integer.
+ *
+ * Results:
+ * None.
+ *
+ * Side effects:
+ * Initializes the bignum supplied, and stores the converted number in
+ * it.
+ *
+ *----------------------------------------------------------------------
+ */
+
+int
+Tcl_InitBignumFromDouble(
+ Tcl_Interp *interp, /* For error message */
+ double d, /* Number to convert */
+ mp_int *b) /* Place to store the result */
+{
+ double fract;
+ int expt;
+
+ /*
+ * Infinite values can't convert to bignum.
+ */
+
+ if (TclIsInfinite(d)) {
+ if (interp != NULL) {
+ const char *s = "integer value too large to represent";
+
+ Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1));
+ Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL);
+ }
+ return TCL_ERROR;
+ }
+
+ fract = frexp(d,&expt);
+ if (expt <= 0) {
+ mp_init(b);
+ mp_zero(b);
+ } else {
+ Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
+ int shift = expt - mantBits;
+
+ TclBNInitBignumFromWideInt(b, w);
+ if (shift < 0) {
+ mp_div_2d(b, -shift, b, NULL);
+ } else if (shift > 0) {
+ mp_mul_2d(b, shift, b);
+ }
+ }
+ return TCL_OK;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclBignumToDouble --
+ *
+ * Convert an arbitrary-precision integer to a native floating point
+ * number.
+ *
+ * Results:
+ * Returns the converted number. Sets errno to ERANGE if the number is
+ * too large to convert.
+ *
+ *----------------------------------------------------------------------
+ */
+
+double
+TclBignumToDouble(
+ mp_int *a) /* Integer to convert. */
+{
+ mp_int b;
+ int bits, shift, i;
+ double r;
+
+ /*
+ * Determine how many bits we need, and extract that many from the input.
+ * Round to nearest unit in the last place.
+ */
+
+ bits = mp_count_bits(a);
+ if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
+ errno = ERANGE;
+ if (a->sign == MP_ZPOS) {
+ return HUGE_VAL;
+ } else {
+ return -HUGE_VAL;
+ }
+ }
+ shift = mantBits + 1 - bits;
+ mp_init(&b);
+ if (shift > 0) {
+ mp_mul_2d(a, shift, &b);
+ } else if (shift < 0) {
+ mp_div_2d(a, -shift, &b, NULL);
+ } else {
+ mp_copy(a, &b);
+ }
+ mp_add_d(&b, 1, &b);
+ mp_div_2d(&b, 1, &b, NULL);
+
+ /*
+ * Accumulate the result, one mp_digit at a time.
+ */
+
+ r = 0.0;
+ for (i=b.used-1 ; i>=0 ; --i) {
+ r = ldexp(r, DIGIT_BIT) + b.dp[i];
+ }
+ mp_clear(&b);
+
+ /*
+ * Scale the result to the correct number of bits.
+ */
+
+ r = ldexp(r, bits - mantBits);
+
+ /*
+ * Return the result with the appropriate sign.
+ */
+
+ if (a->sign == MP_ZPOS) {
+ return r;
+ } else {
+ return -r;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclCeil --
+ *
+ * Computes the smallest floating point number that is at least the
+ * mp_int argument.
+ *
+ * Results:
+ * Returns the floating point number.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+double
+TclCeil(
+ mp_int *a) /* Integer to convert. */
+{
+ double r = 0.0;
+ mp_int b;
+
+ mp_init(&b);
+ if (mp_cmp_d(a, 0) == MP_LT) {
+ mp_neg(a, &b);
+ r = -TclFloor(&b);
+ } else {
+ int bits = mp_count_bits(a);
+
+ if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
+ r = HUGE_VAL;
+ } else {
+ int i, exact = 1, shift = mantBits - bits;
+
+ if (shift > 0) {
+ mp_mul_2d(a, shift, &b);
+ } else if (shift < 0) {
+ mp_int d;
+ mp_init(&d);
+ mp_div_2d(a, -shift, &b, &d);
+ exact = mp_iszero(&d);
+ mp_clear(&d);
+ } else {
+ mp_copy(a, &b);
+ }
+ if (!exact) {
+ mp_add_d(&b, 1, &b);
+ }
+ for (i=b.used-1 ; i>=0 ; --i) {
+ r = ldexp(r, DIGIT_BIT) + b.dp[i];
+ }
+ r = ldexp(r, bits - mantBits);
+ }
+ }
+ mp_clear(&b);
+ return r;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclFloor --
+ *
+ * Computes the largest floating point number less than or equal to
+ * the mp_int argument.
+ *
+ * Results:
+ * Returns the floating point value.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+double
+TclFloor(
+ mp_int *a) /* Integer to convert. */
+{
+ double r = 0.0;
+ mp_int b;
+
+ mp_init(&b);
+ if (mp_cmp_d(a, 0) == MP_LT) {
+ mp_neg(a, &b);
+ r = -TclCeil(&b);
+ } else {
+ int bits = mp_count_bits(a);
+
+ if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
+ r = DBL_MAX;
+ } else {
+ int i, shift = mantBits - bits;
+
+ if (shift > 0) {
+ mp_mul_2d(a, shift, &b);
+ } else if (shift < 0) {
+ mp_div_2d(a, -shift, &b, NULL);
+ } else {
+ mp_copy(a, &b);
+ }
+ for (i=b.used-1 ; i>=0 ; --i) {
+ r = ldexp(r, DIGIT_BIT) + b.dp[i];
+ }
+ r = ldexp(r, bits - mantBits);
+ }
+ }
+ mp_clear(&b);
+ return r;
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * BignumToBiasedFrExp --
+ *
+ * Convert an arbitrary-precision integer to a native floating point
+ * number in the range [0.5,1) times a power of two. NOTE: Intentionally
+ * converts to a number that's a few ulp too small, so that
+ * RefineApproximation will not overflow near the high end of the
+ * machine's arithmetic range.
+ *
+ * Results:
+ * Returns the converted number.
+ *
+ * Side effects:
+ * Stores the exponent of two in 'machexp'.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+BignumToBiasedFrExp(
+ mp_int *a, /* Integer to convert */
+ int *machexp) /* Power of two */
+{
+ mp_int b;
+ int bits;
+ int shift;
+ int i;
+ double r;
+
+ /*
+ * Determine how many bits we need, and extract that many from the input.
+ * Round to nearest unit in the last place.
+ */
+
+ bits = mp_count_bits(a);
+ shift = mantBits - 2 - bits;
+ mp_init(&b);
+ if (shift > 0) {
+ mp_mul_2d(a, shift, &b);
+ } else if (shift < 0) {
+ mp_div_2d(a, -shift, &b, NULL);
+ } else {
+ mp_copy(a, &b);
+ }
+
+ /*
+ * Accumulate the result, one mp_digit at a time.
+ */
+
+ r = 0.0;
+ for (i=b.used-1; i>=0; --i) {
+ r = ldexp(r, DIGIT_BIT) + b.dp[i];
+ }
+ mp_clear(&b);
+
+ /*
+ * Return the result with the appropriate sign.
+ */
+
+ *machexp = bits - mantBits + 2;
+ return ((a->sign == MP_ZPOS) ? r : -r);
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * Pow10TimesFrExp --
+ *
+ * Multiply a power of ten by a number expressed as fraction and
+ * exponent.
+ *
+ * Results:
+ * Returns the significand of the result.
+ *
+ * Side effects:
+ * Overwrites the 'machexp' parameter with the exponent of the result.
+ *
+ * Assumes that 'exponent' is such that 10**exponent would be a double, even
+ * though 'fraction*10**(machexp+exponent)' might overflow.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static double
+Pow10TimesFrExp(
+ int exponent, /* Power of 10 to multiply by */
+ double fraction, /* Significand of multiplicand */
+ int *machexp) /* On input, exponent of multiplicand. On
+ * output, exponent of result. */
+{
+ int i, j;
+ int expt = *machexp;
+ double retval = fraction;
+
+ if (exponent > 0) {
+ /*
+ * Multiply by 10**exponent
+ */
+
+ retval = frexp(retval * 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.
+ *
+ *----------------------------------------------------------------------
+ */
+
+static Tcl_WideUInt
+Nokia770Twiddle(
+ Tcl_WideUInt w) /* Number to transpose */
+{
+ return (((w >> 32) & 0xffffffff) | (w << 32));
+}
+
+/*
+ *----------------------------------------------------------------------
+ *
+ * TclNokia770Doubles --
+ *
+ * Transpose the two words of a number for Nokia 770 floating
+ * point handling.
+ *
+ *----------------------------------------------------------------------
+ */
+
+int
+TclNokia770Doubles(void)
+{
+ return n770_fp;
+}
+
+/*
+ * Local Variables:
+ * mode: c
+ * c-basic-offset: 4
+ * fill-column: 78
+ * End:
+ */
generated by cgit v0.12 at 2024-07-16 14:21:08 (GMT)