Please enter the commit message for your changes. Lines starting

[eos/base.git] / util / src / TclTk / tcl8.6.12 / generic / tclStrToD.c

/*
 * tclStrToD.c --
 *
 *      This file contains a collection of procedures for managing conversions
 *      to/from floating-point in Tcl. They include TclParseNumber, which
 *      parses numbers from strings; TclDoubleDigits, which formats numbers
 *      into strings of digits, and procedures for interconversion among
 *      'double' and 'mp_int' types.
 *
 * 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 <float.h>
#include <math.h>

#ifdef _WIN32
#define copysign _copysign
#endif

/*
 * Define KILL_OCTAL to suppress interpretation of numbers with leading zero
 * as octal. (Ceterum censeo: numeros octonarios delendos esse.)
 */

#undef  KILL_OCTAL

/*
 * This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754
 * floating point; of these, only IEEE-754 can represent NaN. IEEE-754 can be
 * uniquely determined by radix and by the widths of significand and exponent.
 */

#if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024)
#   define IEEE_FLOATING_POINT
#endif

/*
 * Rounding controls. (Thanks a lot, Intel!)
 */

#ifdef __i386
/*
 * gcc on x86 needs access to rounding controls, because of a questionable
 * feature where it retains intermediate results as IEEE 'long double' values
 * somewhat unpredictably. It is tempting to include fpu_control.h, but that
 * file exists only on Linux; it is missing on Cygwin and MinGW. Most gcc-isms
 * and ix86-isms are factored out here.
 */

#if defined(__GNUC__)
typedef unsigned int    fpu_control_t __attribute__ ((__mode__ (__HI__)));

#define _FPU_GETCW(cw)  __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
#define _FPU_SETCW(cw)  __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
#   define FPU_IEEE_ROUNDING    0x027F
#   define ADJUST_FPU_CONTROL_WORD
#define TCL_IEEE_DOUBLE_ROUNDING \
    fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING;    \
    fpu_control_t oldRoundingMode;                      \
    _FPU_GETCW(oldRoundingMode);                        \
    _FPU_SETCW(roundTo53Bits)
#define TCL_DEFAULT_DOUBLE_ROUNDING \
    _FPU_SETCW(oldRoundingMode)

/*
 * Sun ProC needs sunmath for rounding control on x86 like gcc above.
 */
#elif defined(__sun)
#include <sunmath.h>
#define TCL_IEEE_DOUBLE_ROUNDING \
    ieee_flags("set","precision","double",NULL)
#define TCL_DEFAULT_DOUBLE_ROUNDING \
    ieee_flags("clear","precision",NULL,NULL)

/*
 * Other platforms are assumed to always operate in full IEEE mode, so we make
 * the macros to go in and out of that mode do nothing.
 */

#else /* !__GNUC__ && !__sun */
#define TCL_IEEE_DOUBLE_ROUNDING        ((void) 0)
#define TCL_DEFAULT_DOUBLE_ROUNDING     ((void) 0)
#endif
#else /* !__i386 */
#define TCL_IEEE_DOUBLE_ROUNDING        ((void) 0)
#define TCL_DEFAULT_DOUBLE_ROUNDING     ((void) 0)
#endif

/*
 * MIPS floating-point units need special settings in control registers to use
 * gradual underflow as we expect.  This fix is for the MIPSpro compiler.
 */

#if defined(__sgi) && defined(_COMPILER_VERSION)
#include <sys/fpu.h>
#endif

/*
 * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN.
 * Everyone else uses 7ff8000000000000. (Why, HP, why?)
 */

#ifdef __hppa
#   define NAN_START    0x7FF4
#   define NAN_MASK     (((Tcl_WideUInt) 1) << 50)
#else
#   define NAN_START    0x7FF8
#   define NAN_MASK     (((Tcl_WideUInt) 1) << 51)
#endif

/*
 * Constants used by this file (most of which are only ever calculated at
 * runtime).
 */

/* Magic constants */

#define LOG10_2 0.3010299956639812
#define TWO_OVER_3LOG10 0.28952965460216784
#define LOG10_3HALVES_PLUS_FUDGE 0.1760912590558

/*
 * Definitions of the parts of an IEEE754-format floating point number.
 */

#define SIGN_BIT        0x80000000
                                /* Mask for the sign bit in the first word of
                                 * a double. */
#define EXP_MASK        0x7FF00000
                                /* Mask for the exponent field in the first
                                 * word of a double. */
#define EXP_SHIFT       20      /* Shift count to make the exponent an
                                 * integer. */
#define HIDDEN_BIT      (((Tcl_WideUInt) 0x00100000) << 32)
                                /* Hidden 1 bit for the significand. */
#define HI_ORDER_SIG_MASK 0x000FFFFF
                                /* Mask for the high-order part of the
                                 * significand in the first word of a
                                 * double. */
#define SIG_MASK        (((Tcl_WideUInt) HI_ORDER_SIG_MASK << 32) \
                        | 0xFFFFFFFF)
                                /* Mask for the 52-bit significand. */
#define FP_PRECISION    53      /* Number of bits of significand plus the
                                 * hidden bit. */
#define EXPONENT_BIAS   0x3FF   /* Bias of the exponent 0. */

/*
 * Derived quantities.
 */

#define TEN_PMAX        22      /* floor(FP_PRECISION*log(2)/log(5)) */
#define QUICK_MAX       14      /* floor((FP_PRECISION-1)*log(2)/log(10))-1 */
#define BLETCH          0x10    /* Highest power of two that is greater than
                                 * DBL_MAX_10_EXP, divided by 16. */
#define DIGIT_GROUP     8       /* floor(MP_DIGIT_BIT*log(2)/log(10)) */

/*
 * Union used to dismantle floating point numbers.
 */

typedef union Double {
    struct {
#ifdef WORDS_BIGENDIAN
        int word0;
        int word1;
#else
        int word1;
        int word0;
#endif
    } w;
    double d;
    Tcl_WideUInt q;
} Double;

static int maxpow10_wide;       /* The powers of ten that can be represented
                                 * exactly as wide integers. */
static Tcl_WideUInt *pow10_wide;
#define MAXPOW  22
static double pow10vals[MAXPOW+1];
                                /* The powers of ten that can be represented
                                 * exactly as IEEE754 doubles. */
static int mmaxpow;             /* Largest power of ten that can be
                                 * represented exactly in a 'double'. */
static int log10_DIGIT_MAX;     /* The number of decimal digits that fit in an
                                 * mp_digit. */
static int log2FLT_RADIX;       /* Logarithm of the floating point radix. */
static int mantBits;            /* Number of bits in a double's significand */
static mp_int pow5[9];          /* Table of powers of 5**(2**n), up to
                                 * 5**256 */
static double tiny = 0.0;       /* The smallest representable double. */
static int maxDigits;           /* The maximum number of digits to the left of
                                 * the decimal point of a double. */
static int minDigits;           /* The maximum number of digits to the right
                                 * of the decimal point in a double. */
static const double pow_10_2_n[] = {    /* Inexact higher powers of ten. */
    1.0,
    100.0,
    10000.0,
    1.0e+8,
    1.0e+16,
    1.0e+32,
    1.0e+64,
    1.0e+128,
    1.0e+256
};

static int n770_fp;             /* Flag is 1 on Nokia N770 floating point.
                                 * Nokia's floating point has the words
                                 * reversed: if big-endian is 7654 3210,
                                 * and little-endian is       0123 4567,
                                 * then Nokia's FP is         4567 0123;
                                 * little-endian within the 32-bit words but
                                 * big-endian between them. */

/*
 * Table of powers of 5 that are small enough to fit in an mp_digit.
 */

static const mp_digit dpow5[13] = {
               1,              5,             25,            125,
             625,           3125,          15625,          78125,
          390625,        1953125,        9765625,       48828125,
       244140625
};

/*
 * Table of powers: pow5_13[n] = 5**(13*2**(n+1))
 */

static mp_int pow5_13[5];       /* Table of powers: 5**13, 5**26, 5**52,
                                 * 5**104, 5**208 */
static const double tens[] = {
    1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09,
    1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
    1e20, 1e21, 1e22
};

static const int itens [] = {
    1,
    10,
    100,
    1000,
    10000,
    100000,
    1000000,
    10000000,
    100000000
};

static const double bigtens[] = {
    1e016, 1e032, 1e064, 1e128, 1e256
};
#define N_BIGTENS 5

static const int log2pow5[27] = {
    01,  3,  5,  7, 10, 12, 14, 17, 19, 21,
    24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
    47, 49, 52, 54, 56, 59, 61
};
#define N_LOG2POW5 27

static const Tcl_WideUInt wuipow5[27] = {
    (Tcl_WideUInt) 1,           /* 5**0 */
    (Tcl_WideUInt) 5,
    (Tcl_WideUInt) 25,
    (Tcl_WideUInt) 125,
    (Tcl_WideUInt) 625,
    (Tcl_WideUInt) 3125,        /* 5**5 */
    (Tcl_WideUInt) 3125*5,
    (Tcl_WideUInt) 3125*25,
    (Tcl_WideUInt) 3125*125,
    (Tcl_WideUInt) 3125*625,
    (Tcl_WideUInt) 3125*3125,   /* 5**10 */
    (Tcl_WideUInt) 3125*3125*5,
    (Tcl_WideUInt) 3125*3125*25,
    (Tcl_WideUInt) 3125*3125*125,
    (Tcl_WideUInt) 3125*3125*625,
    (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */
    (Tcl_WideUInt) 3125*3125*3125*5,
    (Tcl_WideUInt) 3125*3125*3125*25,
    (Tcl_WideUInt) 3125*3125*3125*125,
    (Tcl_WideUInt) 3125*3125*3125*625,
    (Tcl_WideUInt) 3125*3125*3125*3125, /* 5**20 */
    (Tcl_WideUInt) 3125*3125*3125*3125*5,
    (Tcl_WideUInt) 3125*3125*3125*3125*25,
    (Tcl_WideUInt) 3125*3125*3125*3125*125,
    (Tcl_WideUInt) 3125*3125*3125*3125*625,
    (Tcl_WideUInt) 3125*3125*3125*3125*3125,  /* 5**25 */
    (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */
};

/*
 * Static functions defined in this file.
 */

static int              AccumulateDecimalDigit(unsigned, int,
                            Tcl_WideUInt *, mp_int *, int);
static double           MakeHighPrecisionDouble(int signum,
                            mp_int *significand, int nSigDigs, long exponent);
static double           MakeLowPrecisionDouble(int signum,
                            Tcl_WideUInt significand, int nSigDigs,
                            long exponent);
#ifdef IEEE_FLOATING_POINT
static double           MakeNaN(int signum, Tcl_WideUInt tag);
#endif
static double           RefineApproximation(double approx,
                            mp_int *exactSignificand, int exponent);
static void             MulPow5(mp_int *, unsigned, mp_int *);
static int              NormalizeRightward(Tcl_WideUInt *);
static int              RequiredPrecision(Tcl_WideUInt);
static void             DoubleToExpAndSig(double, Tcl_WideUInt *, int *,
                            int *);
static void             TakeAbsoluteValue(Double *, int *);
static char *           FormatInfAndNaN(Double *, int *, char **);
static char *           FormatZero(int *, char **);
static int              ApproximateLog10(Tcl_WideUInt, int, int);
static int              BetterLog10(double, int, int *);
static void             ComputeScale(int, int, int *, int *, int *, int *);
static void             SetPrecisionLimits(int, int, int *, int *, int *,
                            int *);
static char *           BumpUp(char *, char *, int *);
static int              AdjustRange(double *, int);
static char *           ShorteningQuickFormat(double, int, int, double,
                            char *, int *);
static char *           StrictQuickFormat(double, int, int, double,
                            char *, int *);
static char *           QuickConversion(double, int, int, int, int, int, int,
                            int *, char **);
static void             CastOutPowersOf2(int *, int *, int *);
static char *           ShorteningInt64Conversion(Double *, int, Tcl_WideUInt,
                            int, int, int, int, int, int, int, int, int,
                            int, int, int *, char **);
static char *           StrictInt64Conversion(Double *, int, Tcl_WideUInt,
                            int, int, int, int, int, int,
                            int, int, int *, char **);
static int              ShouldBankerRoundUpPowD(mp_int *, int, int);
static int              ShouldBankerRoundUpToNextPowD(mp_int *, mp_int *,
                            int, int, int, mp_int *);
static char *           ShorteningBignumConversionPowD(Double *dPtr,
                            int convType, Tcl_WideUInt bw, int b2, int b5,
                            int m2plus, int m2minus, int m5,
                            int sd, int k, int len,
                            int ilim, int ilim1, int *decpt,
                            char **endPtr);
static char *           StrictBignumConversionPowD(Double *dPtr, int convType,
                            Tcl_WideUInt bw, int b2, int b5,
                            int sd, int k, int len,
                            int ilim, int ilim1, int *decpt,
                            char **endPtr);
static int              ShouldBankerRoundUp(mp_int *, mp_int *, int);
static int              ShouldBankerRoundUpToNext(mp_int *, mp_int *,
                            mp_int *, int, int, mp_int *);
static char *           ShorteningBignumConversion(Double *dPtr, int convType,
                            Tcl_WideUInt bw, int b2,
                            int m2plus, int m2minus,
                            int s2, int s5, int k, int len,
                            int ilim, int ilim1, int *decpt,
                            char **endPtr);
static char *           StrictBignumConversion(Double *dPtr, int convType,
                            Tcl_WideUInt bw, int b2,
                            int s2, int s5, int k, int len,
                            int ilim, int ilim1, int *decpt,
                            char **endPtr);
static double           BignumToBiasedFrExp(const mp_int *big, int *machexp);
static double           Pow10TimesFrExp(int exponent, double fraction,
                            int *machexp);
static double           SafeLdExp(double fraction, int exponent);
#ifdef IEEE_FLOATING_POINT
static Tcl_WideUInt     Nokia770Twiddle(Tcl_WideUInt w);
#endif
\f
/*
 *----------------------------------------------------------------------
 *
 * 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_BINARY_ONLY:        parse only in the binary format, whether
 *              or not a prefix is present that would lead to binary parsing.
 *              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 (TclIsSpaceProcM(c)) {
                if (flags & TCL_PARSE_NO_WHITESPACE) {
                    goto endgame;
                }
                break;
            } else if (c == '+') {
                state = SIGNUM;
                break;
            } else if (c == '-') {
                signum = 1;
                state = SIGNUM;
                break;
            }
            /* FALLTHROUGH */

        case SIGNUM:
            /*
             * Scanned a leading + or -. Acceptable characters are digits,
             * period, I, and N.
             */

            if (c == '0') {
                if (flags & TCL_PARSE_DECIMAL_ONLY) {
                    state = DECIMAL;
                } else {
                    state = ZERO;
                }
                break;
            } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
                goto zerox;
            } else if (flags & TCL_PARSE_BINARY_ONLY) {
                goto zerob;
            } else if (flags & TCL_PARSE_OCTAL_ONLY) {
                goto zeroo;
            } else if (isdigit(UCHAR(c))) {
                significandWide = c - '0';
                numSigDigs = 1;
                state = DECIMAL;
                break;
            } else if (flags & TCL_PARSE_INTEGER_ONLY) {
                goto endgame;
            } else if (c == '.') {
                state = LEADING_RADIX_POINT;
                break;
            } else if (c == 'I' || c == 'i') {
                state = sI;
                break;
#ifdef IEEE_FLOATING_POINT
            } else if (c == 'N' || c == 'n') {
                state = sN;
                break;
#endif
            }
            goto endgame;

        case ZERO:
            /*
             * Scanned a leading zero (perhaps with a + or -). Acceptable
             * inputs are digits, period, X, b, and E. If 8 or 9 is
             * encountered, the number can't be octal. This state and the
             * OCTAL state differ only in whether they recognize 'X' and 'b'.
             */

            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (c == 'x' || c == 'X') {
                if (flags & (TCL_PARSE_OCTAL_ONLY|TCL_PARSE_BINARY_ONLY)) {
                    goto endgame;
                }
                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') {
                if (flags & TCL_PARSE_OCTAL_ONLY) {
                    goto endgame;
                }
                state = ZERO_B;
                break;
            }
            if (flags & TCL_PARSE_BINARY_ONLY) {
                goto zerob;
            }
            if (c == 'o' || c == 'O') {
                explicitOctal = 1;
                state = ZERO_O;
                break;
            }
#ifdef KILL_OCTAL
            goto decimal;
#endif
            /* FALLTHROUGH */

        case OCTAL:
            /*
             * Scanned an optional + or -, followed by a string of octal
             * digits. Acceptable inputs are more digits, period, or E. If 8
             * or 9 is encountered, commit to floating point.
             */

            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            /* FALLTHROUGH */
        case ZERO_O:
        zeroo:
            if (c == '0') {
                numTrailZeros++;
                state = OCTAL;
                break;
            } else if (c >= '1' && c <= '7') {
                if (objPtr != NULL) {
                    shift = 3 * (numTrailZeros + 1);
                    significandOverflow = AccumulateDecimalDigit(
                            (unsigned)(c-'0'), numTrailZeros,
                            &significandWide, &significandBig,
                            significandOverflow);

                    if (!octalSignificandOverflow) {
                        /*
                         * Shifting by more bits than are in the value being
                         * shifted is at least de facto nonportable. Check for
                         * too large shifts first.
                         */

                        if ((octalSignificandWide != 0)
                                && (((size_t)shift >=
                                        CHAR_BIT*sizeof(Tcl_WideUInt))
                                || (octalSignificandWide >
                                        (~(Tcl_WideUInt)0 >> shift)))) {
                            octalSignificandOverflow = 1;
                            TclBNInitBignumFromWideUInt(&octalSignificandBig,
                                    octalSignificandWide);
                        }
                    }
                    if (!octalSignificandOverflow) {
                        octalSignificandWide =
                                (octalSignificandWide << shift) + (c - '0');
                    } else {
                        mp_mul_2d(&octalSignificandBig, shift,
                                &octalSignificandBig);
                        mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
                                &octalSignificandBig);
                    }
                }
                if (numSigDigs != 0) {
                    numSigDigs += numTrailZeros+1;
                } else {
                    numSigDigs = 1;
                }
                numTrailZeros = 0;
                state = OCTAL;
                break;
            }
            /* FALLTHROUGH */

        case BAD_OCTAL:
            if (explicitOctal) {
                /*
                 * No forgiveness for bad digits in explicitly octal numbers.
                 */

                goto endgame;
            }
            if (flags & TCL_PARSE_INTEGER_ONLY) {
                /*
                 * No seeking floating point when parsing only integer.
                 */

                goto endgame;
            }
#ifndef KILL_OCTAL

            /*
             * Scanned a number with a leading zero that contains an 8, 9,
             * radix point or E. This is an invalid octal number, but might
             * still be floating point.
             */

            if (c == '0') {
                numTrailZeros++;
                state = BAD_OCTAL;
                break;
            } else if (isdigit(UCHAR(c))) {
                if (objPtr != NULL) {
                    significandOverflow = AccumulateDecimalDigit(
                            (unsigned)(c-'0'), numTrailZeros,
                            &significandWide, &significandBig,
                            significandOverflow);
                }
                if (numSigDigs != 0) {
                    numSigDigs += (numTrailZeros + 1);
                } else {
                    numSigDigs = 1;
                }
                numTrailZeros = 0;
                state = BAD_OCTAL;
                break;
            } else if (c == '.') {
                state = FRACTION;
                break;
            } else if (c == 'E' || c == 'e') {
                state = EXPONENT_START;
                break;
            }
#endif
            goto endgame;

            /*
             * Scanned 0x. If state is HEXADECIMAL, scanned at least one
             * character following the 0x. The only acceptable inputs are
             * hexadecimal digits.
             */

        case HEXADECIMAL:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            /* FALLTHROUGH */

        case ZERO_X:
        zerox:
            if (c == '0') {
                numTrailZeros++;
                state = HEXADECIMAL;
                break;
            } else if (isdigit(UCHAR(c))) {
                d = (c-'0');
            } else if (c >= 'A' && c <= 'F') {
                d = (c-'A'+10);
            } else if (c >= 'a' && c <= 'f') {
                d = (c-'a'+10);
            } else {
                goto endgame;
            }
            if (objPtr != NULL) {
                shift = 4 * (numTrailZeros + 1);
                if (!significandOverflow) {
                    /*
                     * Shifting by more bits than are in the value being
                     * shifted is at least de facto nonportable. Check for too
                     * large shifts first.
                     */

                    if (significandWide != 0 &&
                            ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
                            significandWide > (~(Tcl_WideUInt)0 >> shift))) {
                        significandOverflow = 1;
                        TclBNInitBignumFromWideUInt(&significandBig,
                                significandWide);
                    }
                }
                if (!significandOverflow) {
                    significandWide = (significandWide << shift) + d;
                } else {
                    mp_mul_2d(&significandBig, shift, &significandBig);
                    mp_add_d(&significandBig, (mp_digit) d, &significandBig);
                }
            }
            numTrailZeros = 0;
            state = HEXADECIMAL;
            break;

        case BINARY:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            /* FALLTHRU */
        case ZERO_B:
        zerob:
            if (c == '0') {
                numTrailZeros++;
                state = BINARY;
                break;
            } else if (c != '1') {
                goto endgame;
            }
            if (objPtr != NULL) {
                shift = numTrailZeros + 1;
                if (!significandOverflow) {
                    /*
                     * Shifting by more bits than are in the value being
                     * shifted is at least de facto nonportable. Check for too
                     * large shifts first.
                     */

                    if (significandWide != 0 &&
                            ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
                            significandWide > (~(Tcl_WideUInt)0 >> shift))) {
                        significandOverflow = 1;
                        TclBNInitBignumFromWideUInt(&significandBig,
                                significandWide);
                    }
                }
                if (!significandOverflow) {
                    significandWide = (significandWide << shift) + 1;
                } else {
                    mp_mul_2d(&significandBig, shift, &significandBig);
                    mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
                }
            }
            numTrailZeros = 0;
            state = BINARY;
            break;

        case DECIMAL:
            /*
             * Scanned an optional + or - followed by a string of decimal
             * digits.
             */

#ifdef KILL_OCTAL
        decimal:
#endif
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (c == '0') {
                numTrailZeros++;
                state = DECIMAL;
                break;
            } else if (isdigit(UCHAR(c))) {
                if (objPtr != NULL) {
                    significandOverflow = AccumulateDecimalDigit(
                            (unsigned)(c - '0'), numTrailZeros,
                            &significandWide, &significandBig,
                            significandOverflow);
                }
                numSigDigs += numTrailZeros+1;
                numTrailZeros = 0;
                state = DECIMAL;
                break;
            } else if (flags & TCL_PARSE_INTEGER_ONLY) {
                goto endgame;
            } else if (c == '.') {
                state = FRACTION;
                break;
            } else if (c == 'E' || c == 'e') {
                state = EXPONENT_START;
                break;
            }
            goto endgame;

            /*
             * Found a decimal point. If no digits have yet been scanned, E is
             * not allowed; otherwise, it introduces the exponent. If at least
             * one digit has been found, we have a possible complete number.
             */

        case FRACTION:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (c == 'E' || c=='e') {
                state = EXPONENT_START;
                break;
            }
            /* FALLTHROUGH */

        case LEADING_RADIX_POINT:
            if (c == '0') {
                numDigitsAfterDp++;
                numTrailZeros++;
                state = FRACTION;
                break;
            } else if (isdigit(UCHAR(c))) {
                numDigitsAfterDp++;
                if (objPtr != NULL) {
                    significandOverflow = AccumulateDecimalDigit(
                            (unsigned)(c-'0'), numTrailZeros,
                            &significandWide, &significandBig,
                            significandOverflow);
                }
                if (numSigDigs != 0) {
                    numSigDigs += numTrailZeros+1;
                } else {
                    numSigDigs = 1;
                }
                numTrailZeros = 0;
                state = FRACTION;
                break;
            }
            goto endgame;

        case EXPONENT_START:
            /*
             * Scanned the E at the start of an exponent. Make sure a legal
             * character follows before using the C library strtol routine,
             * which allows whitespace.
             */

            if (c == '+') {
                state = EXPONENT_SIGNUM;
                break;
            } else if (c == '-') {
                exponentSignum = 1;
                state = EXPONENT_SIGNUM;
                break;
            }
            /* FALLTHROUGH */

        case EXPONENT_SIGNUM:
            /*
             * Found the E at the start of the exponent, followed by a sign
             * character.
             */

            if (isdigit(UCHAR(c))) {
                exponent = c - '0';
                state = EXPONENT;
                break;
            }
            goto endgame;

        case EXPONENT:
            /*
             * Found an exponent with at least one digit. Accumulate it,
             * making sure to hard-pin it to LONG_MAX on overflow.
             */

            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (isdigit(UCHAR(c))) {
                if (exponent < (LONG_MAX - 9) / 10) {
                    exponent = 10 * exponent + (c - '0');
                } else {
                    exponent = LONG_MAX;
                }
                state = EXPONENT;
                break;
            }
            goto endgame;

            /*
             * Parse out INFINITY by simply spelling it out. INF is accepted
             * as an abbreviation; other prefices are not.
             */

        case sI:
            if (c == 'n' || c == 'N') {
                state = sIN;
                break;
            }
            goto endgame;
        case sIN:
            if (c == 'f' || c == 'F') {
                state = sINF;
                break;
            }
            goto endgame;
        case sINF:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (c == 'i' || c == 'I') {
                state = sINFI;
                break;
            }
            goto endgame;
        case sINFI:
            if (c == 'n' || c == 'N') {
                state = sINFIN;
                break;
            }
            goto endgame;
        case sINFIN:
            if (c == 'i' || c == 'I') {
                state = sINFINI;
                break;
            }
            goto endgame;
        case sINFINI:
            if (c == 't' || c == 'T') {
                state = sINFINIT;
                break;
            }
            goto endgame;
        case sINFINIT:
            if (c == 'y' || c == 'Y') {
                state = sINFINITY;
                break;
            }
            goto endgame;

            /*
             * Parse NaN's.
             */
#ifdef IEEE_FLOATING_POINT
        case sN:
            if (c == 'a' || c == 'A') {
                state = sNA;
                break;
            }
            goto endgame;
        case sNA:
            if (c == 'n' || c == 'N') {
                state = sNAN;
                break;
            }
            goto endgame;
        case sNAN:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            if (c == '(') {
                state = sNANPAREN;
                break;
            }
            goto endgame;

            /*
             * Parse NaN(hexdigits)
             */
        case sNANHEX:
            if (c == ')') {
                state = sNANFINISH;
                break;
            }
            /* FALLTHROUGH */
        case sNANPAREN:
            if (TclIsSpaceProcM(c)) {
                break;
            }
            if (numSigDigs < 13) {
                if (c >= '0' && c <= '9') {
                    d = c - '0';
                } else if (c >= 'a' && c <= 'f') {
                    d = 10 + c - 'a';
                } else if (c >= 'A' && c <= 'F') {
                    d = 10 + c - 'A';
                } else {
                    goto endgame;
                }
                numSigDigs++;
                significandWide = (significandWide << 4) + d;
                state = sNANHEX;
                break;
            }
            goto endgame;
        case sNANFINISH:
#endif

        case sINFINITY:
            acceptState = state;
            acceptPoint = p;
            acceptLen = len;
            goto endgame;
        }
        p++;
        len--;
    }

  endgame:
    if (acceptState == INITIAL) {
        /*
         * No numeric string at all found.
         */

        status = TCL_ERROR;
        if (endPtrPtr != NULL) {
            *endPtrPtr = p;
        }
    } else {
        /*
         * Back up to the last accepting state in the lexer.
         */

        p = acceptPoint;
        len = acceptLen;
        if (!(flags & TCL_PARSE_NO_WHITESPACE)) {
            /*
             * Accept trailing whitespace.
             */

            while (len != 0 && TclIsSpaceProcM(*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:
#endif
            Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'",
                    acceptState, bytes);
        case BINARY:
            shift = numTrailZeros;
            if (!significandOverflow && significandWide != 0 &&
                    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
                    significandWide > (MOST_BITS + signum) >> shift)) {
                significandOverflow = 1;
                TclBNInitBignumFromWideUInt(&significandBig, significandWide);
            }
            if (shift) {
                if (!significandOverflow) {
                    significandWide <<= shift;
                } else {
                    mp_mul_2d(&significandBig, shift, &significandBig);
                }
            }
            goto returnInteger;

        case HEXADECIMAL:
            /*
             * Returning a hex integer. Final scaling step.
             */

            shift = 4 * numTrailZeros;
            if (!significandOverflow && significandWide !=0 &&
                    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
                    significandWide > (MOST_BITS + signum) >> shift)) {
                significandOverflow = 1;
                TclBNInitBignumFromWideUInt(&significandBig, significandWide);
            }
            if (shift) {
                if (!significandOverflow) {
                    significandWide <<= shift;
                } else {
                    mp_mul_2d(&significandBig, shift, &significandBig);
                }
            }
            goto returnInteger;

        case OCTAL:
            /*
             * Returning an octal integer. Final scaling step.
             */

            shift = 3 * numTrailZeros;
            if (!octalSignificandOverflow && octalSignificandWide != 0 &&
                    ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
                    octalSignificandWide > (MOST_BITS + signum) >> shift)) {
                octalSignificandOverflow = 1;
                TclBNInitBignumFromWideUInt(&octalSignificandBig,
                        octalSignificandWide);
            }
            if (shift) {
                if (!octalSignificandOverflow) {
                    octalSignificandWide <<= shift;
                } else {
                    mp_mul_2d(&octalSignificandBig, shift,
                            &octalSignificandBig);
                }
            }
            if (!octalSignificandOverflow) {
                if (octalSignificandWide >
                        (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
#ifndef TCL_WIDE_INT_IS_LONG
                    if (octalSignificandWide <= (MOST_BITS + signum)) {
                        objPtr->typePtr = &tclWideIntType;
                        if (signum) {
                            objPtr->internalRep.wideValue =
                                    - (Tcl_WideInt) octalSignificandWide;
                        } else {
                            objPtr->internalRep.wideValue =
                                    (Tcl_WideInt) octalSignificandWide;
                        }
                        break;
                    }
#endif
                    TclBNInitBignumFromWideUInt(&octalSignificandBig,
                            octalSignificandWide);
                    octalSignificandOverflow = 1;
                } else {
                    objPtr->typePtr = &tclIntType;
                    if (signum) {
                        objPtr->internalRep.longValue =
                                - (long) octalSignificandWide;
                    } else {
                        objPtr->internalRep.longValue =
                                (long) octalSignificandWide;
                    }
                }
            }
            if (octalSignificandOverflow) {
                if (signum) {
                    (void)mp_neg(&octalSignificandBig, &octalSignificandBig);
                }
                TclSetBignumInternalRep(objPtr, &octalSignificandBig);
            }
            break;

        case ZERO:
        case DECIMAL:
            significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
                    &significandWide, &significandBig, significandOverflow);
            if (!significandOverflow && (significandWide > MOST_BITS+signum)) {
                significandOverflow = 1;
                TclBNInitBignumFromWideUInt(&significandBig, significandWide);
            }
        returnInteger:
            if (!significandOverflow) {
                if (significandWide >
                        (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
#ifndef TCL_WIDE_INT_IS_LONG
                    if (significandWide <= MOST_BITS+signum) {
                        objPtr->typePtr = &tclWideIntType;
                        if (signum) {
                            objPtr->internalRep.wideValue =
                                    - (Tcl_WideInt) significandWide;
                        } else {
                            objPtr->internalRep.wideValue =
                                    (Tcl_WideInt) significandWide;
                        }
                        break;
                    }
#endif
                    TclBNInitBignumFromWideUInt(&significandBig,
                            significandWide);
                    significandOverflow = 1;
                } else {
                    objPtr->typePtr = &tclIntType;
                    if (signum) {
                        objPtr->internalRep.longValue =
                                - (long) significandWide;
                    } else {
                        objPtr->internalRep.longValue =
                                (long) significandWide;
                    }
                }
            }
            if (significandOverflow) {
                if (signum) {
                    (void)mp_neg(&significandBig, &significandBig);
                }
                TclSetBignumInternalRep(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) {
                /*
                 * At this point exponent>=0, so the following calculation
                 * cannot underflow.
                 */
                exponent = -exponent;
            }

            /*
             * Adjust the exponent for the number of trailing zeros that
             * have not been accumulated, and the number of digits after
             * the decimal point. Pin any overflow to LONG_MAX/LONG_MIN
             * respectively.
             */

            if (exponent >= 0) {
                if (exponent - numDigitsAfterDp > LONG_MAX - numTrailZeros) {
                    exponent = LONG_MAX;
                } else {
                    exponent = exponent - numDigitsAfterDp + numTrailZeros;
                }
            } else {
                if (exponent + numTrailZeros < LONG_MIN + numDigitsAfterDp) {
                    exponent = LONG_MIN;
                } else {
                    exponent = exponent + numTrailZeros - numDigitsAfterDp;
                }
            }

            /*
             * The desired number is now significandWide * 10**exponent
             * or significandBig * 10**exponent, depending on whether
             * the significand has overflowed a wide int.
             */
            if (!significandOverflow) {
                objPtr->internalRep.doubleValue = MakeLowPrecisionDouble(
                        signum, significandWide, numSigDigs, exponent);
            } else {
                objPtr->internalRep.doubleValue = MakeHighPrecisionDouble(
                        signum, &significandBig, numSigDigs, exponent);
            }
            break;

        case sINF:
        case sINFINITY:
            if (signum) {
                objPtr->internalRep.doubleValue = -HUGE_VAL;
            } else {
                objPtr->internalRep.doubleValue = HUGE_VAL;
            }
            objPtr->typePtr = &tclDoubleType;
            break;

#ifdef IEEE_FLOATING_POINT
        case sNAN:
        case sNANFINISH:
            objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide);
            objPtr->typePtr = &tclDoubleType;
            break;
#endif
        case INITIAL:
            /* This case only to silence compiler warning. */
            Tcl_Panic("TclParseNumber: state INITIAL can't happen here");
        }
    }

    /*
     * Format an error message when an invalid number is encountered.
     */

    if (status != TCL_OK) {
        if (interp != NULL) {
            Tcl_Obj *msg = Tcl_ObjPrintf("expected %s but got \"",
                    expected);

            Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, "");
            Tcl_AppendToObj(msg, "\"", -1);
            if (state == BAD_OCTAL) {
                Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1);
            }
            Tcl_SetObjResult(interp, msg);
            Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", NULL);
        }
    }

    /*
     * Free memory.
     */

    if (octalSignificandOverflow) {
        mp_clear(&octalSignificandBig);
    }
    if (significandOverflow) {
        mp_clear(&significandBig);
    }
    return status;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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 MP_DIGIT_BIT is at least 27. The
         * first multiplication, by up to 10**7, is done with a one-DIGIT
         * multiply (this presumes that MP_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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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 */
    long exponent)              /* Power of ten */
{
    double retval;              /* Value of the number. */
    mp_int significandBig;      /* Significand expressed as a bignum. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.
     */

    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Test for the easy cases.
     */

    if (significand == 0) {
        return copysign(0.0, -signum);
    }
    if (numSigDigs <= QUICK_MAX) {
        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 = QUICK_MAX - numSigDigs;

                if (exponent-diff <= mmaxpow) {
                    /*
                     * 10**exponent is not an exact integer, but
                     * 10**(exponent-diff) is exact, and so is
                     * significand*10**diff, so we can still compute the value
                     * with only one roundoff.
                     */

                    volatile double factor = (double)
                            ((Tcl_WideInt)significand * pow10vals[diff]);
                    retval = factor * pow10vals[exponent-diff];
                    goto returnValue;
                }
            }
        } else {
            if (exponent >= -mmaxpow) {
                /*
                 * 10**-exponent is an exact integer, and so is the
                 * significand. Compute the result by one division, again with
                 * only one rounding.
                 */

                retval = (double)
                        ((Tcl_WideInt)significand / pow10vals[-exponent]);
                goto returnValue;
            }
        }
    }

    /*
     * All the easy cases have failed. Promote ths significand to bignum and
     * call MakeHighPrecisionDouble to do it the hard way.
     */

    TclBNInitBignumFromWideUInt(&significandBig, significand);
    retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
            exponent);
    mp_clear(&significandBig);

    /*
     * Come here to return the computed value.
     */

  returnValue:
    if (signum) {
        retval = -retval;
    }

    /*
     * On gcc on x86, restore the floating point mode word.
     */

    TCL_DEFAULT_DOUBLE_ROUNDING;

    return retval;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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 */
    long exponent)              /* Power of 10 by which to multiply */
{
    double retval;
    int machexp;                /* Machine exponent of a power of 10. */

    /*
     * With gcc on x86, the floating point rounding mode is double-extended.
     * This causes the result of double-precision calculations to be rounded
     * twice: once to the precision of double-extended and then again to the
     * precision of double. Double-rounding introduces gratuitous errors of 1
     * ulp, so we need to change rounding mode to 53-bits.
     */

    TCL_IEEE_DOUBLE_ROUNDING;

    /*
     * Quick checks for zero, and over/underflow. Be careful to avoid
     * integer overflow when calculating with 'exponent'.
     */

    if (mp_iszero(significand)) {
        return copysign(0.0, -signum);
    }
    if (exponent >= 0 && exponent-1 > maxDigits-numSigDigs) {
        retval = HUGE_VAL;
        goto returnValue;
    } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) {
        retval = 0.0;
        goto returnValue;
    }

    /*
     * Develop a first approximation to the significand. It is tempting simply
     * to force bignum to double, but that will overflow on input numbers like
     * 1.[string repeat 0 1000]1; while this is a not terribly likely
     * scenario, we still have to deal with it. Use fraction and exponent
     * instead. Once we have the significand, multiply by 10**exponent. Test
     * for overflow. Convert back to a double, and test for underflow.
     */

    retval = BignumToBiasedFrExp(significand, &machexp);
    retval = Pow10TimesFrExp(exponent, retval, &machexp);
    if (machexp > DBL_MAX_EXP*log2FLT_RADIX) {
        retval = HUGE_VAL;
        goto returnValue;
    }
    retval = SafeLdExp(retval, machexp);
        if (tiny == 0.0) {
            tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
        }
    if (retval < tiny) {
        retval = tiny;
    }

    /*
     * Refine the result twice. (The second refinement should be necessary
     * only if the best approximation is a power of 2 minus 1/2 ulp).
     */

    retval = RefineApproximation(retval, significand, exponent);
    retval = RefineApproximation(retval, significand, exponent);

    /*
     * Come here to return the computed value.
     */

  returnValue:
    if (signum) {
        retval = -retval;
    }

    /*
     * On gcc on x86, restore the floating point mode word.
     */

    TCL_DEFAULT_DOUBLE_ROUNDING;

    return retval;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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
\f
/*
 *----------------------------------------------------------------------
 *
 * 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 roundToEven = 0;        /* Flag == TRUE if we need to invoke
                                 * "round to even" functionality */
    double rteSignificand;      /* Significand of the round-to-even result */
    int rteExponent;            /* Exponent of the round-to-even result */
    int shift;                  /* Shift count for converting numerator
                                 * and denominator of corrector to floating
                                 * point */
    Tcl_WideInt rteSigWide;     /* Wide integer version of the significand
                                 * for testing evenness */
    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;
    }
    significand = frexp(approxResult, &binExponent);

    /*
     * We are trying to compute a corrector term that, when added to the
     * approximate result, will yield close to the exact result.
     * The exact result is exactSignificand * 10**exponent.
     * The approximate result is significand * 2**binExponent
     * If exponent<0, we need to multiply the exact value by 10**-exponent
     * to make it an integer, plus another factor of 2 to decide on rounding.
     *  Similarly if binExponent<FP_PRECISION, we need
     * to multiply by 2**FP_PRECISION to make the approximate value an integer.
     *
     * Let M = 2**M2 * 5**M5 be the least common multiple of these two
     * multipliers.
     */

    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;
        }
    }

    /*
     * Compute twoMv as 2*M*v, where v is the approximate value.
     * This is done by bit-whacking to calculate 2**(M2+1)*significand,
     * and then multiplying by 5**M5.
     */

    msb = binExponent + M2;     /* 1008 */
    nDigits = msb / MP_DIGIT_BIT + 1;
    mp_init_size(&twoMv, nDigits);
    i = (msb % MP_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, MP_DIGIT_BIT);
    }
    for (i = 0; i <= 8; ++i) {
        if (M5 & (1 << i)) {
            mp_mul(&twoMv, pow5+i, &twoMv);
        }
    }

    /*
     * Compute twoMd as 2*M*d, where d is the exact value.
     * This is done by multiplying by 5**(M5+exponent) and then multiplying
     * by 2**(M5+exponent+1), which is, of couse, a left shift.
     */

    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);

    /*
     * Now let twoMd = twoMd - twoMv, the difference between the exact and
     * approximate values.
     */

    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. Place that factor
     * times 1/2 ULP into twoMd.
     */

    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);
    }

    /*
     * Will the eventual correction term be less than, equal to, or
     * greater than 1/2 ULP?
     */

    switch (mp_cmp_mag(&twoMd, &twoMv)) {
    case MP_LT:
        /*
         * If the error is less than 1/2 ULP, there's no correction to make.
         */
        mp_clear(&twoMd);
        mp_clear(&twoMv);
        return approxResult;
    case MP_EQ:
        /*
         * If the error is exactly 1/2 ULP, we need to round to even.
         */
        roundToEven = 1;
        break;
    case MP_GT:
        /*
         * We need to correct the result if the error exceeds 1/2 ULP.
         */
        break;
    }

    /*
     * If we're in the 'round to even' case, and the significand is already
     * even, we're done. Return the approximate result.
     */
    if (roundToEven) {
        rteSignificand = frexp(approxResult, &rteExponent);
        rteSigWide = (Tcl_WideInt) ldexp(rteSignificand, FP_PRECISION);
        if ((rteSigWide & 1) == 0) {
            mp_clear(&twoMd);
            mp_clear(&twoMv);
            return approxResult;
        }
    }

    /*
     * Reduce the numerator and denominator of the corrector term so that
     * they will fit in the floating point precision.
     */
    shift = mp_count_bits(&twoMv) - FP_PRECISION - 1;
    if (shift > 0) {
        mp_div_2d(&twoMv, shift, &twoMv, NULL);
        mp_div_2d(&twoMd, shift, &twoMd, NULL);
    }

    /*
     * 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * MultPow5 --
 *
 *      Multiply a bignum by a power of 5.
 *
 * Side effects:
 *      Stores base*5**n in result.
 *
 *----------------------------------------------------------------------
 */

static inline 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);
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * FormatInfAndNaN --
 *
 *      Bailout for formatting infinities and Not-A-Number.
 *
 * Results:
 *      Returns one of the strings 'Infinity' and 'NaN'.  The string returned
 *      must be freed by the caller using 'ckfree'.
 *
 * Side effects:
 *      Stores 9999 in *decpt, and sets '*endPtr' to designate the terminating
 *      NUL byte of the string if 'endPtr' is not NULL.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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).
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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");
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline char *
QuickConversion(
    double e,                   /* Number to format. */
    int k,                      /* floor(log10(d)), approximately. */
    int k_check,                /* 0 if k is exact, 1 if it may be too high */
    int flags,                  /* Flags passed to dtoa:
                                 *    TCL_DD_SHORTEN_FLAG */
    int len,                    /* Length of the return value. */
    int ilim,                   /* Number of digits to store. */
    int ilim1,                  /* Number of digits to store if we misguessed
                                 * k. */
    int *decpt,                 /* OUTPUT: Location of the decimal point. */
    char **endPtr)              /* OUTPUT: Pointer to the terminal null
                                 * byte. */
{
    int ieps;                   /* Number of 1-ulp roundoff errors that have
                                 * accumulated in the calculation. */
    Double eps;                 /* Estimated roundoff error. */
    char *retval;               /* Returned string. */
    char *end;                  /* Pointer to the terminal null byte in the
                                 * returned string. */
    volatile double d;          /* Workaround for a bug in mingw gcc 3.4.5 */

    /*
     * Bring d into the range [1 .. 10).
     */

    ieps = AdjustRange(&e, k);
    d = e;

    /*
     * If the guessed value of k didn't get d into range, adjust it by one. If
     * that leaves us outside the range in which quick format is accurate,
     * bail out.
     */

    if (k_check && d < 1. && ilim > 0) {
        if (ilim1 < 0) {
            return NULL;
        }
        ilim = ilim1;
        --k;
        d = 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 = 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline void
CastOutPowersOf2(
    int *b2,                    /* Power of 2 to multiply the significand. */
    int *m2,                    /* Power of 2 to multiply 1/2 ulp. */
    int *s2)                    /* Power of 2 to multiply the common
                                 * denominator. */
{
    int i;

    if (*m2 > 0 && *s2 > 0) {   /* Find the smallest power of 2 in the
                                 * numerator. */
        if (*m2 < *s2) {        /* Find the lowest common denominator. */
            i = *m2;
        } else {
            i = *s2;
        }
        *b2 -= i;               /* Reduce to lowest terms. */
        *m2 -= i;
        *s2 -= i;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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**MP_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.
 *
 *----------------------------------------------------------------------
 */

static inline int
ShouldBankerRoundUpPowD(
    mp_int *b,                  /* Numerator of the fraction. */
    int sd,                     /* Denominator is 2**(sd*MP_DIGIT_BIT). */
    int isodd)                  /* 1 if the digit is odd, 0 if even. */
{
    int i;
    static const mp_digit topbit = ((mp_digit)1) << (MP_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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline int
ShouldBankerRoundUpToNextPowD(
    mp_int *b,                  /* Numerator of the fraction. */
    mp_int *m,                  /* Numerator of the rounding tolerance. */
    int sd,                     /* Common denominator is 2**(sd*MP_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**(MP_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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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**MP_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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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(&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*MP_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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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**MP_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'.
 *
 *----------------------------------------------------------------------
 */

static inline 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*MP_DIGIT_BIT)
     * by mp_digit extraction.
     */

    i = 1;
    for (;;) {
        if (b.used <= sd) {
            digit = 0;
        } else {
            digit = b.dp[sd];
            if (b.used > sd+1 || digit >= 10) {
                Tcl_Panic("wrong digit!");
            }
            --b.used;
            mp_clamp(&b);
        }

        /*
         * Have we converted all the requested digits?
         */

        *s++ = '0' + digit;
        if (i == ilim) {
            if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
                s = BumpUp(s, retval, &k);
            }
            while (*--s == '0') {
                /* do nothing */
            }
            ++s;
            break;
        }

        /*
         * Advance to the next digit.
         */

        mp_mul_d(&b, 10, &b);
        ++i;
    }

    /*
     * Endgame - store the location of the decimal point and the end of the
     * string.
     */

    mp_clear_multi(&b, &temp, NULL);
    *s = '\0';
    *decpt = k;
    if (endPtr) {
        *endPtr = s;
    }
    return retval;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline 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(&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(&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_t
             *       arithmetic here if S < 2**64/10. But it's a win only for
             *       a fairly narrow range of magnitudes so perhaps not worth
             *       bothering.  We already know that we shorten the
             *       denominator by at least 1 mp_digit, perhaps 2, as we do
             *       the conversion for 17 digits of significance.
             * Possible savings:
             * 10**26   1 trip through loop before fallback possible
             * 10**27   1 trip
             * 10**28   2 trips
             * 10**29   3 trips
             * 10**30   4 trips
             * 10**31   5 trips
             * 10**32   6 trips
             * 10**33   7 trips
             * 10**34   8 trips
             * 10**35   9 trips
             * 10**36  10 trips
             * 10**37  11 trips
             * 10**38  12 trips
             * 10**39  13 trips
             * 10**40  14 trips
             * 10**41  15 trips
             * 10**42  16 trips
             * thereafter no gain.
             */
        } else {
            mp_mul_d(&b, 10, &b);
            mp_mul_d(&mminus, 10, &mminus);
            if (m2plus > m2minus) {
                mp_mul_2d(&mplus, 10, &mplus);
            }
        }

        ++i;
    }

    /*
     * Endgame - store the location of the decimal point and the end of the
     * string.
     */

    if (m2plus > m2minus) {
        mp_clear(&mplus);
    }
    mp_clear_multi(&b, &mminus, &temp, &dig, &S, NULL);
    *s = '\0';
    *decpt = k;
    if (endPtr) {
        *endPtr = s;
    }
    return retval;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.
 *
 *----------------------------------------------------------------------
 */

static inline char *
StrictBignumConversion(
    Double *dPtr,               /* Original number being converted. */
    int convType,               /* Conversion type. */
    Tcl_WideUInt bw,            /* Integer significand and exponent. */
    int b2,                     /* Scale factor for the significand. */
    int s2, int s5,             /* Scale factors for denominator. */
    int k,                      /* Guessed position of the decimal point. */
    int len,                    /* Size of the digit buffer to allocate. */
    int ilim,                   /* Number of digits to convert if b >= s */
    int ilim1,                  /* Number of digits to convert if b < s */
    int *decpt,                 /* OUTPUT: Position of the decimal point. */
    char **endPtr)              /* OUTPUT: Pointer to the end of the number */
{
    char *retval = ckalloc(len+1);
                                /* Buffer of digits to return. */
    char *s = retval;           /* Cursor in the return value. */
    mp_int b;                   /* Numerator of the result. */
    mp_int S;                   /* Denominator of the result. */
    mp_int dig;                 /* Current digit of the result. */
    int digit;                  /* Current digit of the result. */
    mp_int temp;                /* Work area. */
    int g;                      /* Size of the current digit ground. */
    int i, j;

    /*
     * b = bw * 2**b2 * 5**b5
     * S = 2**s2 * 5*s5
     */

    mp_init_multi(&temp, &dig, NULL);
    TclBNInitBignumFromWideUInt(&b, bw);
    mp_mul_2d(&b, b2, &b);
    mp_init_set(&S, 1);
    MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);

    /*
     * Handle the case where we guess the position of the decimal point wrong.
     */

    if (mp_cmp_mag(&b, &S) == MP_LT) {
        mp_mul_d(&b, 10, &b);
        ilim =ilim1;
        --k;
    }

    /*
     * Convert the leading digit.
     */

    i = 0;
    mp_div(&b, &S, &dig, &b);
    if (dig.used > 1 || dig.dp[0] >= 10) {
        Tcl_Panic("wrong digit!");
    }
    digit = dig.dp[0];

    /*
     * Is a single digit all that was requested?
     */

    *s++ = '0' + digit;
    if (++i >= ilim) {
        mp_mul_2d(&b, 1, &b);
        if (ShouldBankerRoundUp(&b, &S, digit&1)) {
            s = BumpUp(s, retval, &k);
        }
    } else {
        for (;;) {
            /*
             * Shift by a group of digits.
             */

            g = ilim - i;
            if (g > DIGIT_GROUP) {
                g = DIGIT_GROUP;
            }
            if (s5 >= g) {
                mp_div_d(&S, dpow5[g], &S, NULL);
                s5 -= g;
            } else if (s5 > 0) {
                mp_div_d(&S, dpow5[s5], &S, NULL);
                mp_mul_d(&b, dpow5[g - s5], &b);
                s5 = 0;
            } else {
                mp_mul_d(&b, dpow5[g], &b);
            }
            mp_mul_2d(&b, g, &b);

            /*
             * As with the shortening bignum conversion, it's possible at this
             * point that we will have reduced the denominator to less than
             * 2**64/10, at which point it would be possible to fall back to
             * to int64_t arithmetic. But the potential payoff is tremendously
             * less - unless we're working in F format - because we know that
             * three groups of digits will always suffice for %#.17e, the
             * longest format that doesn't introduce empty precision.
             *
             * Extract the next group of digits.
             */

            mp_div(&b, &S, &dig, &b);
            if (dig.used > 1) {
                Tcl_Panic("wrong digit!");
            }
            digit = dig.dp[0];
            for (j = g-1; j >= 0; --j) {
                int t = itens[j];

                *s++ = digit / t + '0';
                digit %= t;
            }
            i += g;

            /*
             * Have we converted all the requested digits?
             */

            if (i == ilim) {
                mp_mul_2d(&b, 1, &b);
                if (ShouldBankerRoundUp(&b, &S, digit&1)) {
                    s = BumpUp(s, retval, &k);
                }
                break;
            }
        }
    }
    while (*--s == '0') {
        /* do nothing */
    }
    ++s;

    /*
     * Endgame - store the location of the decimal point and the end of the
     * string.
     */

    mp_clear_multi(&b, &S, &temp, &dig, NULL);
    *s = '\0';
    *decpt = k;
    if (endPtr) {
        *endPtr = s;
    }
    return retval;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * TclDoubleDigits --
 *
 *      Core of Tcl's conversion of double-precision floating point numbers to
 *      decimal.
 *
 * Results:
 *      Returns a newly-allocated string of digits.
 *
 * Side effects:
 *      Sets *decpt to the index of the character in the string before the
 *      place that the decimal point should go. If 'endPtr' is not NULL, sets
 *      endPtr to point to the terminating '\0' byte of the string. Sets *sign
 *      to 1 if a minus sign should be printed with the number, or 0 if a plus
 *      sign (or no sign) should appear.
 *
 * This function is a service routine that produces the string of digits for
 * floating-point-to-decimal conversion. It can do a number of things
 * according to the 'flags' argument. Valid values for 'flags' include:
 *      TCL_DD_SHORTEST - This is the default for floating point conversion if
 *              ::tcl_precision is 0. It constructs the shortest string of
 *              digits that will reconvert to the given number when scanned.
 *              For floating point numbers that are exactly between two
 *              decimal numbers, it resolves using the 'round to even' rule.
 *              With this value, the 'ndigits' parameter is ignored.
 *      TCL_DD_STEELE - This value is not recommended and may be removed in
 *              the future. It follows the conversion algorithm outlined in
 *              "How to Print Floating-Point Numbers Accurately" by Guy
 *              L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90,
 *              pp. 112-126]. This rule has the effect of rendering 1e23 as
 *              9.9999999999999999e22 - which is a 'better' approximation in
 *              the sense that it will reconvert correctly even if a
 *              subsequent input conversion is 'round up' or 'round down'
 *              rather than 'round to nearest', but is surprising otherwise.
 *      TCL_DD_E_FORMAT - This value is used to prepare numbers for %e format
 *              conversion (or for default floating->string if tcl_precision
 *              is not 0). It constructs a string of at most 'ndigits' digits,
 *              choosing the one that is closest to the given number (and
 *              resolving ties with 'round to even').  It is allowed to return
 *              fewer than 'ndigits' if the number converts exactly; if the
 *              TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG is supplied instead, it
 *              also returns fewer digits if the shorter string will still
 *              reconvert without loss to the given input number. In any case,
 *              strings of trailing zeroes are suppressed.
 *      TCL_DD_F_FORMAT - This value is used to prepare numbers for %f format
 *              conversion. It requests that conversion proceed until
 *              'ndigits' digits after the decimal point have been converted.
 *              It is possible for this format to result in a zero-length
 *              string if the number is sufficiently small. Again, it is
 *              permissible for TCL_DD_F_FORMAT to return fewer digits for a
 *              number that converts exactly, and changing the argument to
 *              TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow the routine
 *              also to return fewer digits if the shorter string will still
 *              reconvert without loss to the given input number. Strings of
 *              trailing zeroes are suppressed.
 *
 *      To any of these flags may be OR'ed TCL_DD_NO_QUICK; this flag requires
 *      all calculations to be done in exact arithmetic. Normally, E and F
 *      format with fewer than about 14 digits will be done with a quick
 *      floating point approximation and fall back on the exact arithmetic
 *      only if the input number is close enough to the midpoint between two
 *      decimal strings that more precision is needed to resolve which string
 *      is correct.
 *
 * The value stored in the 'decpt' argument on return may be negative
 * (indicating that the decimal point falls to the left of the string) or
 * greater than the length of the string. In addition, the value -9999 is used
 * as a sentinel to indicate that the string is one of the special values
 * "Infinity" and "NaN", and that no decimal point should be inserted.
 *
 *----------------------------------------------------------------------
 */

char *
TclDoubleDigits(
    double dv,                  /* Number to convert. */
    int ndigits,                /* Number of digits requested. */
    int flags,                  /* Conversion flags. */
    int *decpt,                 /* OUTPUT: Position of the decimal point. */
    int *sign,                  /* OUTPUT: 1 if the result is negative. */
    char **endPtr)              /* OUTPUT: If not NULL, receives a pointer to
                                 *         one character beyond the end of the
                                 *         returned string. */
{
    int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK);
                                /* Type of conversion being performed:
                                 * TCL_DD_SHORTEST0, TCL_DD_STEELE0,
                                 * TCL_DD_E_FORMAT, or TCL_DD_F_FORMAT. */
    Double d;                   /* Union for deconstructing doubles. */
    Tcl_WideUInt bw;            /* Integer significand. */
    int be;                     /* Power of 2 by which b must be multiplied */
    int bbits;                  /* Number of bits needed to represent b. */
    int denorm;                 /* Flag == 1 iff the input number was
                                 * denormalized. */
    int k;                      /* Estimate of floor(log10(d)). */
    int k_check;                /* Flag == 1 if d is near enough to a power of
                                 * ten that k must be checked. */
    int b2, b5, s2, s5;         /* Powers of 2 and 5 in the numerator and
                                 * denominator of intermediate results. */
    int ilim = -1, ilim1 = -1;  /* Number of digits to convert, and number to
                                 * convert if log10(d) has been
                                 * overestimated. */
    char *retval;               /* Return value from this function. */
    int i = -1;

    /*
     * Put the input number into a union for bit-whacking.
     */

    d.d = dv;

    /*
     * Handle the cases of negative numbers (by taking the absolute value:
     * this includes -Inf and -NaN!), infinity, Not a Number, and zero.
     */

    TakeAbsoluteValue(&d, sign);
    if ((d.w.word0 & EXP_MASK) == EXP_MASK) {
        return FormatInfAndNaN(&d, decpt, endPtr);
    }
    if (d.d == 0.0) {
        return FormatZero(decpt, endPtr);
    }

    /*
     * Unpack the floating point into a wide integer and an exponent.
     * Determine the number of bits that the big integer requires, and compute
     * a quick approximation (which may be one too high) of ceil(log10(d.d)).
     */

    denorm = ((d.w.word0 & EXP_MASK) == 0);
    DoubleToExpAndSig(d.d, &bw, &be, &bbits);
    k = ApproximateLog10(bw, be, bbits);
    k = BetterLog10(d.d, k, &k_check);

    /* At this point, we have:
     *    d is the number to convert.
     *    bw are significand and exponent: d == bw*2**be,
     *    bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits
     *    k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0 if we
     *      know that k is exactly ceil(log10(d)) and 1 if we need to check.
     *    We want a rational number
     *      r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5),
     *    with b2, b5, s2, s5 >= 0.  Note that the most significant decimal
     *    digit is floor(r) and that successive digits can be obtained by
     *    setting r <- 10*floor(r) (or b <= 10 * (b % S)).  Find appropriate
     *    b2, b5, s2, s5.
     */

    ComputeScale(be, k, &b2, &b5, &s2, &s5);

    /*
     * Correct an incorrect caller-supplied 'ndigits'.  Also determine:
     *  i = The maximum number of decimal digits that will be returned in the
     *      formatted string.  This is k + 1 + ndigits for F format, 18 for
     *      shortest and Steele, and ndigits for E format.
     *  ilim = The number of significant digits to convert if k has been
     *         guessed correctly. This is -1 for shortest and Steele (which
     *         stop when all significance has been lost), 'ndigits' for E
     *         format, and 'k + 1 + ndigits' for F format.
     *  ilim1 = The minimum number of significant digits to convert if k has
     *          been guessed 1 too high. This, too, is -1 for shortest and
     *          Steele, and 'ndigits' for E format, but it's 'ndigits-1' for F
     *          format.
     */

    SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1);

    /*
     * Try to do low-precision conversion in floating point rather than
     * resorting to expensive multiprecision arithmetic.
     */

    if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) {
        retval = QuickConversion(d.d, k, k_check, flags, i, ilim, ilim1,
                decpt, endPtr);
        if (retval != NULL) {
            return retval;
        }
    }

    /*
     * For shortening conversions, determine the upper and lower bounds for
     * the remainder at which we can stop.
     *   m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the high
     *        side, and
     *   m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the low
     *        side.
     * We may need to increase s2 to put m2plus, m2minus, b2 over a common
     * denominator.
     */

    if (flags & TCL_DD_SHORTEN_FLAG) {
        int m2minus = b2;
        int m2plus;
        int m5 = b5;
        int len = i;

        /*
         * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5) is 1/2 unit
         * in the least significant place of the floating point number.
         */

        if (denorm) {
            i = be + EXPONENT_BIAS + (FP_PRECISION-1);
        } else {
            i = 1 + FP_PRECISION - bbits;
        }
        b2 += i;
        s2 += i;

        /*
         * Reduce the fractions to lowest terms, since the above calculation
         * may have left excess powers of 2 in numerator and denominator.
         */

        CastOutPowersOf2(&b2, &m2minus, &s2);

        /*
         * In the special case where bw==1, the nearest floating point number
         * to it on the low side is 1/4 ulp below it. Adjust accordingly.
         */

        m2plus = m2minus;
        if (!denorm && bw == 1) {
            ++b2;
            ++s2;
            ++m2plus;
        }

        if (s5+1 < N_LOG2POW5 && s2+1 + log2pow5[s5+1] <= 64) {
            /*
             * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit word,
             * then all our intermediate calculations can be done using exact
             * 64-bit arithmetic with no need for expensive multiprecision
             * operations. (This will be true for all numbers in the range
             * [1.0e-3 .. 1.0e+24]).
             */

            return ShorteningInt64Conversion(&d, convType, bw, b2, b5, m2plus,
                    m2minus, m5, s2, s5, k, len, ilim, ilim1, decpt, endPtr);
        } else if (s5 == 0) {
            /*
             * The denominator is a power of 2, so we can replace division by
             * digit shifts. First we round up s2 to a multiple of MP_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 % MP_DIGIT_BIT != 0) {
                int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT);

                b2 += delta;
                m2plus += delta;
                m2minus += delta;
                s2 += delta;
            }
            return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5,
                    m2plus, m2minus, m5, s2/MP_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 MP_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 % MP_DIGIT_BIT != 0) {
                int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT);

                b2 += delta;
                s2 += delta;
            }
            return StrictBignumConversionPowD(&d, convType, bw, b2, b5,
                    s2/MP_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);
        }
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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.));
    log10_DIGIT_MAX = (int) floor(MP_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
}
\f
/*
 *----------------------------------------------------------------------
 *
 * TclFinalizeDoubleConversion --
 *
 *      Cleans up this file on exit.
 *
 * Results:
 *      None
 *
 * Side effects:
 *      Memory allocated by TclInitDoubleConversion is freed.
 *
 *----------------------------------------------------------------------
 */

void
TclFinalizeDoubleConversion(void)
{
    int i;

    ckfree(pow10_wide);
    for (i=0; i<9; ++i) {
        mp_clear(pow5 + i);
    }
    for (i=0; i < 5; ++i) {
        mp_clear(pow5_13 + i);
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * TclBignumToDouble --
 *
 *      Convert an arbitrary-precision integer to a native floating point
 *      number.
 *
 * Results:
 *      Returns the converted number. Sets errno to ERANGE if the number is
 *      too large to convert.
 *
 *----------------------------------------------------------------------
 */

double
TclBignumToDouble(
    const mp_int *a)                    /* Integer to convert. */
{
    mp_int b;
    int bits, shift, i, lsb;
    double r;


    /*
     * We need a 'mantBits'-bit significand.  Determine what shift will
     * give us that.
     */

    bits = mp_count_bits(a);
    if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
        errno = ERANGE;
        if (mp_isneg(a)) {
            return -HUGE_VAL;
        } else {
            return HUGE_VAL;
        }
    }
    shift = mantBits - bits;

    /*
     * If shift > 0, shift the significand left by the requisite number of
     * bits.  If shift == 0, the significand is already exactly 'mantBits'
     * in length.  If shift < 0, we will need to shift the significand right
     * by the requisite number of bits, and round it. If the '1-shift'
     * least significant bits are 0, but the 'shift'th bit is nonzero,
     * then the significand lies exactly between two values and must be
     * 'rounded to even'.
     */

    mp_init(&b);
    if (shift == 0) {
        mp_copy(a, &b);
    } else if (shift > 0) {
        mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
        lsb = mp_cnt_lsb(a);
        if (lsb == -1-shift) {

            /*
             * Round to even
             */

            mp_div_2d(a, -shift, &b, NULL);
            if (mp_isodd(&b)) {
                if (mp_isneg(&b)) {
                    mp_sub_d(&b, 1, &b);
                } else {
                    mp_add_d(&b, 1, &b);
                }
            }
        } else {

            /*
             * Ordinary rounding
             */

            mp_div_2d(a, -1-shift, &b, NULL);
            if (mp_isneg(&b)) {
                mp_sub_d(&b, 1, &b);
            } else {
                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, MP_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 (mp_isneg(a)) {
        return -r;
    } else {
        return r;
    }
}
\f
/*
 *----------------------------------------------------------------------
 *
 * TclCeil --
 *
 *      Computes the smallest floating point number that is at least the
 *      mp_int argument.
 *
 * Results:
 *      Returns the floating point number.
 *
 *----------------------------------------------------------------------
 */

double
TclCeil(
    const mp_int *a)                    /* Integer to convert. */
{
    double r = 0.0;
    mp_int b;

    mp_init(&b);
    if (mp_cmp_d(a, 0) == MP_LT) {
        mp_neg(a, &b);
        r = -TclFloor(&b);
    } else {
        int bits = mp_count_bits(a);

        if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
            r = HUGE_VAL;
        } else {
            int i, exact = 1, shift = mantBits - bits;

            if (shift > 0) {
                mp_mul_2d(a, shift, &b);
            } else if (shift < 0) {
                mp_int d;
                mp_init(&d);
                mp_div_2d(a, -shift, &b, &d);
                exact = mp_iszero(&d);
                mp_clear(&d);
            } else {
                mp_copy(a, &b);
            }
            if (!exact) {
                mp_add_d(&b, 1, &b);
            }
            for (i=b.used-1 ; i>=0 ; --i) {
                r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
            }
            r = ldexp(r, bits - mantBits);
        }
    }
    mp_clear(&b);
    return r;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * TclFloor --
 *
 *      Computes the largest floating point number less than or equal to the
 *      mp_int argument.
 *
 * Results:
 *      Returns the floating point value.
 *
 *----------------------------------------------------------------------
 */

double
TclFloor(
    const mp_int *a)                    /* Integer to convert. */
{
    double r = 0.0;
    mp_int b;

    mp_init(&b);
    if (mp_cmp_d(a, 0) == MP_LT) {
        mp_neg(a, &b);
        r = -TclCeil(&b);
    } else {
        int bits = mp_count_bits(a);

        if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
            r = DBL_MAX;
        } else {
            int i, shift = mantBits - bits;

            if (shift > 0) {
                mp_mul_2d(a, shift, &b);
            } else if (shift < 0) {
                mp_div_2d(a, -shift, &b, NULL);
            } else {
                mp_copy(a, &b);
            }
            for (i=b.used-1 ; i>=0 ; --i) {
                r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
            }
            r = ldexp(r, bits - mantBits);
        }
    }
    mp_clear(&b);
    return r;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * BignumToBiasedFrExp --
 *
 *      Convert an arbitrary-precision integer to a native floating point
 *      number in the range [0.5,1) times a power of two. NOTE: Intentionally
 *      converts to a number that's a few ulp too small, so that
 *      RefineApproximation will not overflow near the high end of the
 *      machine's arithmetic range.
 *
 * Results:
 *      Returns the converted number.
 *
 * Side effects:
 *      Stores the exponent of two in 'machexp'.
 *
 *----------------------------------------------------------------------
 */

static double
BignumToBiasedFrExp(
    const mp_int *a,            /* Integer to convert. */
    int *machexp)               /* Power of two. */
{
    mp_int b;
    int bits;
    int shift;
    int i;
    double r;

    /*
     * Determine how many bits we need, and extract that many from the input.
     * Round to nearest unit in the last place.
     */

    bits = mp_count_bits(a);
    shift = mantBits - 2 - bits;
    mp_init(&b);
    if (shift > 0) {
        mp_mul_2d(a, shift, &b);
    } else if (shift < 0) {
        mp_div_2d(a, -shift, &b, NULL);
    } else {
        mp_copy(a, &b);
    }

    /*
     * Accumulate the result, one mp_digit at a time.
     */

    r = 0.0;
    for (i=b.used-1; i>=0; --i) {
        r = ldexp(r, MP_DIGIT_BIT) + b.dp[i];
    }
    mp_clear(&b);

    /*
     * Return the result with the appropriate sign.
     */

    *machexp = bits - mantBits + 2;
    return (mp_isneg(a) ? -r : r);
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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;
}
\f
/*
 *----------------------------------------------------------------------
 *
 * 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 */
}
\f
/*
 *----------------------------------------------------------------------
 *
 * Nokia770Twiddle --
 *
 *      Transpose the two words of a number for Nokia 770 floating point
 *      handling.
 *
 *----------------------------------------------------------------------
 */
#ifdef IEEE_FLOATING_POINT
static Tcl_WideUInt
Nokia770Twiddle(
    Tcl_WideUInt w)             /* Number to transpose. */
{
    return (((w >> 32) & 0xFFFFFFFF) | (w << 32));
}
#endif
\f
/*
 *----------------------------------------------------------------------
 *
 * TclNokia770Doubles --
 *
 *      Transpose the two words of a number for Nokia 770 floating point
 *      handling.
 *
 *----------------------------------------------------------------------
 */

int
TclNokia770Doubles(void)
{
    return n770_fp;
}
\f
/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 4
 * fill-column: 78
 * End:
 */