/*
* 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 © 2005 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 "tclTomMath.h"
#include <float.h>
#include <math.h>
#ifdef _WIN32
#define copysign _copysign
#endif
#ifndef PRIx64
# define PRIx64 TCL_LL_MODIFIER "x"
#endif
/*
* 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_DECL \
fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING; \
fpu_control_t oldRoundingMode;
# define TCL_IEEE_DOUBLE_ROUNDING \
_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_DECL
# define TCL_IEEE_DOUBLE_ROUNDING \
ieee_flags("set","precision","double",NULL)
# define TCL_DEFAULT_DOUBLE_ROUNDING \
ieee_flags("clear","precision",NULL,NULL)
# endif
#endif
/*
* 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.
*/
#ifndef TCL_IEEE_DOUBLE_ROUNDING /* !__i386 || (!__GNUC__ && !__sun) */
# define TCL_IEEE_DOUBLE_ROUNDING_DECL
# 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[] = {
(Tcl_WideUInt) 1U, /* 5**0 */
(Tcl_WideUInt) 5U,
(Tcl_WideUInt) 25U,
(Tcl_WideUInt) 125U,
(Tcl_WideUInt) 625U,
(Tcl_WideUInt) 3125U, /* 5**5 */
(Tcl_WideUInt) 3125U*5U,
(Tcl_WideUInt) 3125U*25U,
(Tcl_WideUInt) 3125U*125U,
(Tcl_WideUInt) 3125U*625U,
(Tcl_WideUInt) 3125U*3125U, /* 5**10 */
(Tcl_WideUInt) 3125U*3125U*5U,
(Tcl_WideUInt) 3125U*3125U*25U,
(Tcl_WideUInt) 3125U*3125U*125U,
(Tcl_WideUInt) 3125U*3125U*625U,
(Tcl_WideUInt) 3125U*3125U*3125U, /* 5**15 */
(Tcl_WideUInt) 3125U*3125U*3125U*5U,
(Tcl_WideUInt) 3125U*3125U*3125U*25U,
(Tcl_WideUInt) 3125U*3125U*3125U*125U,
(Tcl_WideUInt) 3125U*3125U*3125U*625U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U, /* 5**20 */
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*5U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*25U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*125U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*625U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*3125U, /* 5**25 */
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*3125U*5U,
(Tcl_WideUInt) 3125U*3125U*3125U*3125U*3125U*25U /* 5**27 */
};
/*
* 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 mp_err MulPow5(mp_int *, unsigned, mp_int *) MP_WUR;
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 *, Tcl_WideUInt,
int, int, int, int, int, int, int, int, int,
int, int, int *, char **);
static char * StrictInt64Conversion(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, mp_int *);
static char * ShorteningBignumConversionPowD(Double *dPtr,
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(
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);
static char * ShorteningBignumConversion(Double *dPtr,
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(
Tcl_WideUInt bw, int b2,
int s2, int s5, int k, int len,
int ilim, int ilim1, int *decpt,
char **endPtr);
static double BignumToBiasedFrExp(const mp_int *big, int *machexp);
static double Pow10TimesFrExp(int exponent, double fraction,
int *machexp);
static double SafeLdExp(double fraction, int exponent);
#ifdef IEEE_FLOATING_POINT
static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
#endif
/*
*----------------------------------------------------------------------
*
* TclParseNumber --
*
* Scans bytes, interpreted as characters in Tcl's internal encoding, and
* parses the longest prefix that is the string representation of a
* number in a format recognized by Tcl.
*
* The arguments bytes, numBytes, and objPtr are the inputs which
* determine the string to be parsed. If bytes is non-NULL, it points to
* the first byte to be scanned. If bytes is NULL, then objPtr must be
* non-NULL, and the string representation of objPtr will be scanned
* (generated first, if necessary). The numBytes argument determines the
* number of bytes to be scanned. If numBytes is TCL_INDEX_NONE, the first NUL
* byte encountered will terminate the scan. Otherwise,
* 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. */
Tcl_Size 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, ZERO_D, BINARY,
HEXADECIMAL, OCTAL, DECIMAL,
LEADING_RADIX_POINT, FRACTION,
EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
#ifdef IEEE_FLOATING_POINT
, sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
#endif
} state = INITIAL;
enum State acceptState = INITIAL;
int signum = 0; /* Sign of the number being parsed. */
Tcl_WideUInt significandWide = 0;
/* Significand of the number being parsed (if
* no overflow). */
mp_int significandBig; /* Significand of the number being parsed (if
* it overflows significandWide). */
int significandOverflow = 0;/* Flag==1 iff significandBig is used. */
Tcl_WideUInt octalSignificandWide = 0;
/* Significand of an octal number; needed
* because we don't know whether a number with
* a leading zero is octal or decimal until
* we've scanned forward to a '.' or 'e'. */
mp_int octalSignificandBig; /* Significand of octal number once
* octalSignificandWide overflows. */
int octalSignificandOverflow = 0;
/* Flag==1 if octalSignificandBig is used. */
int numSigDigs = 0; /* Number of significant digits in the decimal
* significand. */
int numTrailZeros = 0; /* Number of trailing zeroes at the current
* point in the parse. */
int numDigitsAfterDp = 0; /* Number of digits scanned after the decimal
* point. */
int exponentSignum = 0; /* Signum of the exponent of a floating point
* number. */
long exponent = 0; /* Exponent of a floating point number. */
const char *p; /* Pointer to next character to scan. */
Tcl_Size len; /* Number of characters remaining after p. */
const char *acceptPoint; /* Pointer to position after last character in
* an acceptable number. */
Tcl_Size 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 */
mp_err err = MP_OKAY;
#define MOST_BITS (UWIDE_MAX >> 1)
/*
* Initialize bytes to start of the object's string rep if the caller
* didn't pass anything else.
*/
if (bytes == NULL) {
if (interp == NULL && endPtrPtr == NULL) {
if (TclHasInternalRep(objPtr, &tclDictType)) {
/* A dict can never be a (single) number */
return TCL_ERROR;
}
if (TclHasInternalRep(objPtr, &tclListType)) {
Tcl_Size length;
/* A list can only be a (single) number if its length == 1 */
TclListObjLength(NULL, objPtr, &length);
if (length != 1) {
return TCL_ERROR;
}
}
}
bytes = TclGetString(objPtr);
}
p = bytes;
len = numBytes;
acceptPoint = p;
acceptLen = len;
while (1) {
char c = len ? *p : '\0';
/*
* Filter out Numeric Whitespace. Expects:
*
* ::digit:: '_' ::digit::
*
* Verify current '_' is ok, then move on to next character,
* otherwise follow through on to error.
*/
if (c == '_' && !(flags & TCL_PARSE_NO_UNDERSCORE)) {
const char *before, *after;
if (p==bytes) {
/* Not allowed at beginning */
goto endgame;
}
/*
* span multiple numeric whitespace
* V
* example: 5___6
*/
for (before = (p - 1);
(before && *before == '_');
before = (before > p ? (before - 1) : NULL));
for (after = (p + 1);
(after && *after && *after == '_');
after = (*after && *after == '_') ? (after + 1) : NULL);
switch (state) {
case ZERO_B:
case BINARY:
if ((before && (*before != '0' && *before != '1')) ||
(after && (*after != '0' && *after != '1'))) {
/* Not a valid digit */
goto endgame;
}
break;
case ZERO_O:
case OCTAL:
if (((before && (*before < '0' || '7' < *before))) ||
((after && (*after < '0' || '7' < *after)))) {
goto endgame;
}
break;
case FRACTION:
case ZERO:
case ZERO_D:
case DECIMAL:
case LEADING_RADIX_POINT:
case EXPONENT_START:
case EXPONENT_SIGNUM:
case EXPONENT:
if ((!before || isdigit(UCHAR(*before))) &&
(!after || isdigit(UCHAR(*after)))) {
break;
}
if (after && *after=='(') {
/* could be function */
goto continue_num;
}
goto endgame;
case ZERO_X:
case HEXADECIMAL:
if ( (!before || isxdigit(UCHAR(*before))) &&
(!after || isxdigit(UCHAR(*after)))) {
break;
}
goto endgame;
default:
/*
* Not whitespace, but could be legal for other reasons.
* Continue number processing for current character.
*/
goto continue_num;
}
/* Valid whitespace found, move on to the next character */
goto next;
}
continue_num:
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') {
state = ZERO_O;
break;
}
if (c == 'd' || c == 'D') {
state = ZERO_D;
break;
}
goto decimal;
case OCTAL:
/*
* Scanned an optional + or -, followed by a string of octal
* digits. Acceptable inputs are more digits, period, or E. If 8
* or 9 is encountered, commit to floating point.
*/
acceptState = state;
acceptPoint = p;
acceptLen = len;
/* FALLTHROUGH */
case ZERO_O:
zeroo:
if (c == '0') {
numTrailZeros++;
state = OCTAL;
break;
} else if (c >= '1' && c <= '7') {
if (objPtr != NULL) {
shift = 3 * (numTrailZeros + 1);
significandOverflow = AccumulateDecimalDigit(
(unsigned)(c-'0'), numTrailZeros,
&significandWide, &significandBig,
significandOverflow);
if (!octalSignificandOverflow) {
/*
* Shifting by as many or more bits than are in the
* value being shifted is undefined behavior. Check
* for too large shifts first.
*/
if ((octalSignificandWide != 0)
&& (((size_t)shift >=
CHAR_BIT*sizeof(Tcl_WideUInt))
|| (octalSignificandWide >
(UWIDE_MAX >> shift)))) {
octalSignificandOverflow = 1;
err = mp_init_u64(&octalSignificandBig,
octalSignificandWide);
}
}
if (!octalSignificandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (octalSignificandWide != 0) {
octalSignificandWide <<= shift;
}
octalSignificandWide += c - '0';
} else {
if (err == MP_OKAY) {
err = mp_mul_2d(&octalSignificandBig, shift,
&octalSignificandBig);
}
if (err == MP_OKAY) {
err = mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
&octalSignificandBig);
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
}
if (numSigDigs != 0) {
numSigDigs += numTrailZeros+1;
} else {
numSigDigs = 1;
}
numTrailZeros = 0;
state = OCTAL;
break;
}
goto endgame;
/*
* Scanned 0x. If state is HEXADECIMAL, scanned at least one
* character following the 0x. The only acceptable inputs are
* hexadecimal digits.
*/
case HEXADECIMAL:
acceptState = state;
acceptPoint = p;
acceptLen = len;
/* FALLTHROUGH */
case ZERO_X:
zerox:
if (c == '0') {
numTrailZeros++;
state = HEXADECIMAL;
break;
} else if (isdigit(UCHAR(c))) {
d = (c-'0');
} else if (c >= 'A' && c <= 'F') {
d = (c-'A'+10);
} else if (c >= 'a' && c <= 'f') {
d = (c-'a'+10);
} else {
goto endgame;
}
if (objPtr != NULL) {
shift = 4 * (numTrailZeros + 1);
if (!significandOverflow) {
/*
* Shifting by as many or more bits than are in the
* value being shifted is undefined behavior. Check
* for too large shifts first.
*/
if (significandWide != 0 &&
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > (UWIDE_MAX >> shift))) {
significandOverflow = 1;
err = mp_init_u64(&significandBig,
significandWide);
}
}
if (!significandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (significandWide != 0) {
significandWide <<= shift;
}
significandWide += d;
} else if (err == MP_OKAY) {
err = mp_mul_2d(&significandBig, shift, &significandBig);
if (err == MP_OKAY) {
err = mp_add_d(&significandBig, (mp_digit) d, &significandBig);
}
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
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 as many or more bits than are in the
* value being shifted is undefined behavior. Check
* for too large shifts first.
*/
if (significandWide != 0 &&
((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
significandWide > (UWIDE_MAX >> shift))) {
significandOverflow = 1;
err = mp_init_u64(&significandBig,
significandWide);
}
}
if (!significandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (significandWide != 0) {
significandWide <<= shift;
}
significandWide += 1;
} else if (err == MP_OKAY) {
err = mp_mul_2d(&significandBig, shift, &significandBig);
if (err == MP_OKAY) {
err = mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
}
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
numTrailZeros = 0;
state = BINARY;
break;
case ZERO_D:
if (c == '0') {
numTrailZeros++;
} else if ( ! isdigit(UCHAR(c))) {
goto endgame;
}
state = DECIMAL;
flags |= TCL_PARSE_INTEGER_ONLY;
/* FALLTHROUGH */
case DECIMAL:
/*
* Scanned an optional + or - followed by a string of decimal
* digits.
*/
decimal:
acceptState = state;
acceptPoint = p;
acceptLen = len;
if (c == '0') {
numTrailZeros++;
state = DECIMAL;
break;
} else if (isdigit(UCHAR(c))) {
if (objPtr != NULL) {
significandOverflow = AccumulateDecimalDigit(
(unsigned)(c - '0'), numTrailZeros,
&significandWide, &significandBig,
significandOverflow);
}
numSigDigs += numTrailZeros+1;
numTrailZeros = 0;
state = DECIMAL;
break;
} else if (flags & TCL_PARSE_INTEGER_ONLY) {
goto endgame;
} else if (c == '.') {
state = FRACTION;
break;
} else if (c == 'E' || c == 'e') {
state = EXPONENT_START;
break;
}
goto endgame;
/*
* Found a decimal point. If no digits have yet been scanned, E is
* not allowed; otherwise, it introduces the exponent. If at least
* one digit has been found, we have a possible complete number.
*/
case FRACTION:
acceptState = state;
acceptPoint = p;
acceptLen = len;
if (c == 'E' || c=='e') {
state = EXPONENT_START;
break;
}
/* FALLTHROUGH */
case LEADING_RADIX_POINT:
if (c == '0') {
numDigitsAfterDp++;
numTrailZeros++;
state = FRACTION;
break;
} else if (isdigit(UCHAR(c))) {
numDigitsAfterDp++;
if (objPtr != NULL) {
significandOverflow = AccumulateDecimalDigit(
(unsigned)(c-'0'), numTrailZeros,
&significandWide, &significandBig,
significandOverflow);
}
if (numSigDigs != 0) {
numSigDigs += numTrailZeros+1;
} else {
numSigDigs = 1;
}
numTrailZeros = 0;
state = FRACTION;
break;
}
goto endgame;
case EXPONENT_START:
/*
* Scanned the E at the start of an exponent. Make sure a legal
* character follows before using the C library strtol routine,
* which allows whitespace.
*/
if (c == '+') {
state = EXPONENT_SIGNUM;
break;
} else if (c == '-') {
exponentSignum = 1;
state = EXPONENT_SIGNUM;
break;
}
/* FALLTHROUGH */
case EXPONENT_SIGNUM:
/*
* Found the E at the start of the exponent, followed by a sign
* character.
*/
if (isdigit(UCHAR(c))) {
exponent = c - '0';
state = EXPONENT;
break;
}
goto endgame;
case EXPONENT:
/*
* Found an exponent with at least one digit. Accumulate it,
* making sure to hard-pin it to LONG_MAX on overflow.
*/
acceptState = state;
acceptPoint = p;
acceptLen = len;
if (isdigit(UCHAR(c))) {
if (exponent < (LONG_MAX - 9) / 10) {
exponent = 10 * exponent + (c - '0');
} else {
exponent = LONG_MAX;
}
state = EXPONENT;
break;
}
goto endgame;
/*
* Parse out INFINITY by simply spelling it out. INF is accepted
* as an abbreviation; other prefices are not.
*/
case sI:
if (c == 'n' || c == 'N') {
state = sIN;
break;
}
goto endgame;
case sIN:
if (c == 'f' || c == 'F') {
state = sINF;
break;
}
goto endgame;
case sINF:
acceptState = state;
acceptPoint = p;
acceptLen = len;
if (c == 'i' || c == 'I') {
state = sINFI;
break;
}
goto endgame;
case sINFI:
if (c == 'n' || c == 'N') {
state = sINFIN;
break;
}
goto endgame;
case sINFIN:
if (c == 'i' || c == 'I') {
state = sINFINI;
break;
}
goto endgame;
case sINFINI:
if (c == 't' || c == 'T') {
state = sINFINIT;
break;
}
goto endgame;
case sINFINIT:
if (c == 'y' || c == 'Y') {
state = sINFINITY;
break;
}
goto endgame;
/*
* Parse NaN's.
*/
#ifdef IEEE_FLOATING_POINT
case sN:
if (c == 'a' || c == 'A') {
state = sNA;
break;
}
goto endgame;
case sNA:
if (c == 'n' || c == 'N') {
state = sNAN;
break;
}
goto endgame;
case sNAN:
acceptState = state;
acceptPoint = p;
acceptLen = len;
if (c == '(') {
state = sNANPAREN;
break;
}
goto endgame;
/*
* Parse NaN(hexdigits)
*/
case sNANHEX:
if (c == ')') {
state = sNANFINISH;
break;
}
/* FALLTHROUGH */
case sNANPAREN:
if (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;
}
next:
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.
* If the last char seen is the numeric whitespace character '_',
* backup to that.
*/
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 + 1 > 1) || (*p != '\0'))) {
status = TCL_ERROR;
}
} else {
*endPtrPtr = p;
}
}
/*
* Generate and store the appropriate internal rep.
*/
if (status == TCL_OK && objPtr != NULL) {
TclFreeInternalRep(objPtr);
switch (acceptState) {
case SIGNUM:
case ZERO_X:
case ZERO_O:
case ZERO_B:
case ZERO_D:
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;
err = mp_init_u64(&significandBig, significandWide);
}
if (shift) {
if (!significandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (significandWide != 0) {
significandWide <<= shift;
}
} else if (err == MP_OKAY) {
err = mp_mul_2d(&significandBig, shift, &significandBig);
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
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;
err = mp_init_u64(&significandBig, significandWide);
}
if (shift) {
if (!significandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (significandWide != 0) {
significandWide <<= shift;
}
} else if (err == MP_OKAY) {
err = mp_mul_2d(&significandBig, shift, &significandBig);
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
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;
err = mp_init_u64(&octalSignificandBig,
octalSignificandWide);
}
if (shift) {
if (!octalSignificandOverflow) {
/*
* When the significand is 0, it is possible for the
* amount to be shifted to equal or exceed the width
* of the significand. Do not shift when the
* significand is 0 to avoid undefined behavior.
*/
if (octalSignificandWide != 0) {
octalSignificandWide <<= shift;
}
} else if (err == MP_OKAY) {
err = mp_mul_2d(&octalSignificandBig, shift,
&octalSignificandBig);
}
}
if (!octalSignificandOverflow) {
if ((err == MP_OKAY) && (octalSignificandWide > (MOST_BITS + signum))) {
err = mp_init_u64(&octalSignificandBig,
octalSignificandWide);
octalSignificandOverflow = 1;
} else {
objPtr->typePtr = &tclIntType;
if (signum) {
objPtr->internalRep.wideValue =
(Tcl_WideInt)(-octalSignificandWide);
} else {
objPtr->internalRep.wideValue =
(Tcl_WideInt)octalSignificandWide;
}
}
}
if ((err == MP_OKAY) && octalSignificandOverflow) {
if (signum) {
err = mp_neg(&octalSignificandBig, &octalSignificandBig);
}
TclSetBignumInternalRep(objPtr, &octalSignificandBig);
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
break;
case ZERO:
case DECIMAL:
significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
&significandWide, &significandBig, significandOverflow);
if ((err == MP_OKAY) && !significandOverflow && (significandWide > MOST_BITS+signum)) {
significandOverflow = 1;
err = mp_init_u64(&significandBig, significandWide);
}
returnInteger:
if (!significandOverflow) {
if ((err == MP_OKAY) && (significandWide > MOST_BITS+signum)) {
err = mp_init_u64(&significandBig,
significandWide);
significandOverflow = 1;
} else {
objPtr->typePtr = &tclIntType;
if (signum) {
objPtr->internalRep.wideValue =
(Tcl_WideInt)(-significandWide);
} else {
objPtr->internalRep.wideValue =
(Tcl_WideInt)significandWide;
}
}
}
if ((err == MP_OKAY) && significandOverflow) {
if (signum) {
err = mp_neg(&significandBig, &significandBig);
}
TclSetBignumInternalRep(objPtr, &significandBig);
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
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_Size argc;
const char **argv;
if ((TclMaxListLength(bytes, TCL_INDEX_NONE, NULL) > 1)
&& Tcl_SplitList(NULL, bytes, &argc, &argv) == TCL_OK) {
Tcl_Free(argv);
Tcl_AppendToObj(msg, "a list", -1);
} else {
Tcl_AppendToObj(msg, "\"", -1);
Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, "");
Tcl_AppendToObj(msg, "\"", -1);
}
Tcl_SetObjResult(interp, msg);
Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", (char *)NULL);
}
}
/*
* Free memory.
*/
if (octalSignificandOverflow) {
mp_clear(&octalSignificandBig);
}
if (significandOverflow) {
mp_clear(&significandBig);
}
return status;
}
/*
*----------------------------------------------------------------------
*
* AccumulateDecimalDigit --
*
* Consume a decimal digit in a number being scanned.
*
* Results:
* Returns 1 if the number has overflowed to a bignum, 0 if it still fits
* in a wide integer.
*
* Side effects:
* Updates either the wide or bignum representation.
*
*----------------------------------------------------------------------
*/
static int
AccumulateDecimalDigit(
unsigned digit, /* Digit being scanned. */
int numZeros, /* Count of zero digits preceding the digit
* being scanned. */
Tcl_WideUInt *wideRepPtr, /* Representation of the partial number as a
* wide integer. */
mp_int *bignumRepPtr, /* Representation of the partial number as a
* bignum. */
int bignumFlag) /* Flag == 1 if the number overflowed previous
* to this digit. */
{
int i, n;
Tcl_WideUInt w;
/*
* Try wide multiplication first.
*/
if (!bignumFlag) {
w = *wideRepPtr;
if (w == 0) {
/*
* There's no need to multiply if the multiplicand is zero.
*/
*wideRepPtr = digit;
return 0;
} else if (numZeros >= maxpow10_wide
|| w > (UWIDE_MAX-digit)/pow10_wide[numZeros+1]) {
/*
* Wide multiplication will overflow. Expand the number to a
* bignum and fall through into the bignum case.
*/
if (mp_init_u64(bignumRepPtr, w) != MP_OKAY) {
return 0;
}
} 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.
*/
if ((mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
bignumRepPtr) != MP_OKAY)
|| (mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr) != MP_OKAY))
return 0;
} else {
mp_err err;
/*
* 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;
err = mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
for (i = 3; (err == MP_OKAY) && (i <= 7); ++i) {
if (n & (1 << i)) {
err = mp_mul(bignumRepPtr, pow5+i, bignumRepPtr);
}
}
while ((err == MP_OKAY) && (n >= 256)) {
err = mp_mul(bignumRepPtr, pow5+8, bignumRepPtr);
n -= 256;
}
if ((err != MP_OKAY)
|| (mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr) != MP_OKAY)
|| (mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr) != MP_OKAY)) {
return 0;
}
}
return 1;
}
/*
*----------------------------------------------------------------------
*
* MakeLowPrecisionDouble --
*
* Makes the double precision number, signum*significand*10**exponent.
*
* Results:
* Returns the constructed number.
*
* Common cases, where there are few enough digits that the number can be
* represented with at most roundoff, are handled specially here. If the
* number requires more than one rounded operation to compute, the code
* promotes the significand to a bignum and calls MakeHighPrecisionDouble
* to do it instead.
*
*----------------------------------------------------------------------
*/
static double
MakeLowPrecisionDouble(
int signum, /* 1 if the number is negative, 0 otherwise */
Tcl_WideUInt significand, /* Significand of the number */
int numSigDigs, /* Number of digits in the significand */
long exponent) /* Power of ten */
{
TCL_IEEE_DOUBLE_ROUNDING_DECL
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. We also make
* 'retval' volatile, so that it doesn't get promoted to a register.
*/
volatile double retval; /* Value of the number. */
/*
* Test for zero significand, which requires explicit construction
* of -0.0. (Unary minus returns a positive zero.)
*/
if (significand == 0) {
return copysign(0.0, -signum);
}
/*
* Set the FP control word for 53 bits, WARNING: It must be reset
* before returning.
*/
TCL_IEEE_DOUBLE_ROUNDING;
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 the significand to bignum and
* call MakeHighPrecisionDouble to do it the hard way.
*/
if (mp_init_u64(&significandBig, significand) != MP_OKAY) {
return 0.0;
}
retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
exponent);
mp_clear(&significandBig);
/*
* Come here to return the computed value.
*/
returnValue:
if (signum) {
retval = -retval;
}
/*
* On gcc on x86, restore the floating point mode word.
*/
TCL_DEFAULT_DOUBLE_ROUNDING;
return retval;
}
/*
*----------------------------------------------------------------------
*
* MakeHighPrecisionDouble --
*
* Makes the double precision number, signum*significand*10**exponent.
*
* Results:
* Returns the constructed number.
*
* MakeHighPrecisionDouble is used when arbitrary-precision arithmetic is
* needed to ensure correct rounding. It begins by calculating a
* low-precision approximation to the desired number, and then refines
* the answer in high precision.
*
*----------------------------------------------------------------------
*/
static double
MakeHighPrecisionDouble(
int signum, /* 1=negative, 0=nonnegative */
mp_int *significand, /* Exact significand of the number */
int numSigDigs, /* Number of significant digits */
long exponent) /* Power of 10 by which to multiply */
{
TCL_IEEE_DOUBLE_ROUNDING_DECL
int machexp = 0; /* 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. We also make
* 'retval' volatile to make sure that it doesn't get promoted to a
* register.
*/
volatile double retval;
/*
* A zero significand requires explicit construction of -0.0.
* (Unary minus returns positive zero.)
*/
if (mp_iszero(significand)) {
return copysign(0.0, -signum);
}
/*
* Set the 53-bit rounding mode. WARNING: It must be reset before
* returning.
*/
TCL_IEEE_DOUBLE_ROUNDING;
/*
* Make quick checks for over/underflow. Be careful to avoid
* integer overflow when calculating with 'exponent'.
*/
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;
}
/*
*----------------------------------------------------------------------
*
* MakeNaN --
*
* Makes a "Not a Number" given a set of bits to put in the tag bits
*
* Note that a signalling NaN is never returned.
*
*----------------------------------------------------------------------
*/
#ifdef IEEE_FLOATING_POINT
static double
MakeNaN(
int signum, /* Sign bit (1=negative, 0=nonnegative. */
Tcl_WideUInt tags) /* Tag bits to put in the NaN. */
{
union {
Tcl_WideUInt iv;
double dv;
} theNaN;
theNaN.iv = tags;
theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
if (signum) {
theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48;
} else {
theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
}
if (n770_fp) {
theNaN.iv = Nokia770Twiddle(theNaN.iv);
}
return theNaN.dv;
}
#endif
/*
*----------------------------------------------------------------------
*
* RefineApproximation --
*
* Given a poor approximation to a floating point number, returns a
* better one. (The better approximation is correct to within 1 ulp, and
* is entirely correct if the poor approximation is correct to 1 ulp.)
*
* Results:
* Returns the improved result.
*
*----------------------------------------------------------------------
*/
static double
RefineApproximation(
double approxResult, /* Approximate result of conversion. */
mp_int *exactSignificand, /* Integer significand. */
int exponent) /* Power of 10 to multiply by significand. */
{
int M2, M5; /* Powers of 2 and of 5 needed to put the
* decimal and binary numbers over a common
* denominator. */
double significand; /* Sigificand of the binary number. */
int binExponent; /* Exponent of the binary number. */
int msb; /* Most significant bit position of an
* intermediate result. */
int nDigits; /* Number of mp_digit's in an intermediate
* result. */
mp_int twoMv; /* Approx binary value expressed as an exact
* integer scaled by the multiplier 2M. */
mp_int twoMd; /* Exact decimal value expressed as an exact
* integer scaled by the multiplier 2M. */
int scale; /* Scale factor for M. */
int multiplier; /* Power of two to scale M. */
double num, den; /* Numerator and denominator of the correction
* term. */
double quot; /* Correction term. */
double minincr; /* Lower bound on the absolute value of the
* correction term. */
int 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;
mp_err err = MP_OKAY;
/*
* 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;
if (mp_init_size(&twoMv, nDigits) != MP_OKAY) {
return approxResult;
}
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) != MP_OKAY)) {
mp_clear(&twoMv);
return approxResult;
}
}
/*
* 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 course, a left shift.
*/
if (mp_init_copy(&twoMd, exactSignificand) != MP_OKAY) {
mp_clear(&twoMv);
return approxResult;
}
for (i = 0; (i <= 8); ++i) {
if ((M5 + exponent) & (1 << i)) {
err = mp_mul(&twoMd, pow5+i, &twoMd);
}
}
if (err == MP_OKAY) {
err = mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
}
/*
* Now let twoMd = twoMd - twoMv, the difference between the exact and
* approximate values.
*/
if (err == MP_OKAY) {
err = 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_u64(&twoMv, 1);
for (i = 0; (i <= 8) && (err == MP_OKAY); ++i) {
if (M5 & (1 << i)) {
err = mp_mul(&twoMv, pow5+i, &twoMv);
}
}
multiplier = M2 + scale + 1;
if (err != MP_OKAY) {
mp_clear(&twoMd);
mp_clear(&twoMv);
return approxResult;
} else if (multiplier > 0) {
err = mp_mul_2d(&twoMv, multiplier, &twoMv);
} else if (multiplier < 0) {
err = mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
}
if (err != MP_OKAY) {
mp_clear(&twoMd);
mp_clear(&twoMv);
return approxResult;
}
/*
* 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) {
err = mp_div_2d(&twoMv, shift, &twoMv, NULL);
if (err == MP_OKAY) {
err = mp_div_2d(&twoMd, shift, &twoMd, NULL);
}
}
if (err != MP_OKAY) {
mp_clear(&twoMd);
mp_clear(&twoMv);
return approxResult;
}
/*
* Convert the numerator and denominator of the corrector term accurately
* to floating point numbers.
*/
num = TclBignumToDouble(&twoMd);
den = TclBignumToDouble(&twoMv);
quot = SafeLdExp(num/den, scale);
minincr = SafeLdExp(1.0, binExponent-mantBits);
if (quot<0. && quot>-minincr) {
quot = -minincr;
} else if (quot>0. && quot<minincr) {
quot = minincr;
}
mp_clear(&twoMd);
mp_clear(&twoMv);
return approxResult + quot;
}
/*
*----------------------------------------------------------------------
*
* MultPow5 --
*
* Multiply a bignum by a power of 5.
*
* Side effects:
* Stores base*5**n in result.
*
*----------------------------------------------------------------------
*/
static inline mp_err
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;
mp_err err = MP_OKAY;
if (r != 0) {
err = mp_mul_d(p, dpow5[r], result);
p = result;
}
r = 0;
while ((err == MP_OKAY) && (n13 != 0)) {
if (n13 & 1) {
err = mp_mul(p, pow5_13+r, result);
p = result;
}
n13 >>= 1;
++r;
}
if ((err == MP_OKAY) && (p != result)) {
err = mp_copy(p, result);
}
return err;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
/*
*----------------------------------------------------------------------
*
* RequiredPrecision --
*
* Determines the number of bits needed to hold an integer.
*
* Results:
* Returns the position of the most significant bit (0 - 63). Returns 0
* if the number is zero.
*
*----------------------------------------------------------------------
*/
static int
RequiredPrecision(
Tcl_WideUInt w) /* Number to interrogate. */
{
int rv;
unsigned long wi;
if (w & ((Tcl_WideUInt) 0xFFFFFFFF << 32)) {
wi = (unsigned long) (w >> 32); rv = 32;
} else {
wi = (unsigned long) w; rv = 0;
}
if (wi & 0xFFFF0000) {
wi >>= 16; rv += 16;
}
if (wi & 0xFF00) {
wi >>= 8; rv += 8;
}
if (wi & 0xF0) {
wi >>= 4; rv += 4;
}
if (wi & 0xC) {
wi >>= 2; rv += 2;
}
if (wi & 0x2) {
wi >>= 1; ++rv;
}
if (wi & 0x1) {
++rv;
}
return rv;
}
/*
*----------------------------------------------------------------------
*
* DoubleToExpAndSig --
*
* Separates a 'double' into exponent and significand.
*
* Side effects:
* Stores the significand in '*significand' and the exponent in '*expon'
* so that dv == significand * 2.0**expon, and significand is odd. Also
* stores the position of the leftmost 1-bit in 'significand' in 'bits'.
*
*----------------------------------------------------------------------
*/
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;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
}
/*
*----------------------------------------------------------------------
*
* 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 'Tcl_Free'.
*
* 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 = (char *)Tcl_Alloc(9);
strcpy(retval, "Infinity");
if (endPtr) {
*endPtr = retval + 8;
}
} else {
retval = (char *)Tcl_Alloc(4);
strcpy(retval, "NaN");
if (endPtr) {
*endPtr = retval + 3;
}
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* FormatZero --
*
* Bailout to format a zero floating-point number.
*
* Results:
* Returns the constant string "0"
*
* Side effects:
* Stores 1 in '*decpt' and puts a pointer to the NUL byte terminating
* the string in '*endPtr' if 'endPtr' is not NULL.
*
*----------------------------------------------------------------------
*/
static inline char *
FormatZero(
int *decpt, /* Location of the decimal point. */
char **endPtr) /* Pointer to the end of the formatted data */
{
char *retval = (char *)Tcl_Alloc(2);
strcpy(retval, "0");
if (endPtr) {
*endPtr = retval+1;
}
*decpt = 0;
return retval;
}
/*
*----------------------------------------------------------------------
*
* ApproximateLog10 --
*
* Computes a two-term Taylor series approximation to the common log of a
* number, and computes the number's binary log.
*
* Results:
* Return an approximation to floor(log10(bw*2**be)) that is either exact
* or 1 too high.
*
*----------------------------------------------------------------------
*/
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;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
}
/*
*----------------------------------------------------------------------
*
* 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 flags, /* Type of conversion: TCL_DD_SHORTEST,
* 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 (flags & TCL_DD_CONVERSION_TYPE_MASK) {
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:
*iLimPtr = *iLim1Ptr = -1;
*iPtr = 18;
*ndigitsPtr = 0;
break;
}
}
/*
*----------------------------------------------------------------------
*
* 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;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
/*
*----------------------------------------------------------------------
*
* 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;
}
}
/*
*----------------------------------------------------------------------
*
* 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;
}
}
/*
*----------------------------------------------------------------------
*
* 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_SHORTEST */
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 = (char *)Tcl_Alloc(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 {
Tcl_Free(retval);
return NULL;
}
}
/*
* Format the digit string.
*/
if (flags & TCL_DD_SHORTEST) {
end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt);
} else {
end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt);
}
if (end == NULL) {
Tcl_Free(retval);
return NULL;
}
*end = '\0';
if (endPtr != NULL) {
*endPtr = end;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* CastOutPowersOf2 --
*
* Adjust the factors 'b2', 'm2', and 's2' to cast out common powers of 2
* from numerator and denominator in preparation for the 'bignum' method
* of floating point conversion.
*
*----------------------------------------------------------------------
*/
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;
}
}
/*
*----------------------------------------------------------------------
*
* 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. */
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 = (char *)Tcl_Alloc(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 roundoff of being exact?
*/
if (b < mplus || (b == mplus
&& (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
&& (dPtr->w.word1 & 1) == 0)) {
if (digit == 9) {
*s++ = '9';
s = BumpUp(s, retval, &k);
break;
}
++digit;
*s++ = '0' + digit;
break;
}
/*
* Have we converted all the requested digits?
*/
*s++ = '0' + digit;
if (i == ilim) {
if (2*b > S || (2*b == S && (digit & 1) != 0)) {
s = BumpUp(s, retval, &k);
}
break;
}
/*
* Advance to the next digit.
*/
b = 10 * b;
mplus = 10 * mplus;
mminus = 10 * mminus;
++i;
}
/*
* Endgame - store the location of the decimal point and the end of the
* string.
*/
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* StrictInt64Conversion --
*
* Converts a double-precision number to a fixed-length string of 'ilim'
* digits that reconverts exactly to the given number. ('ilim' should be
* replaced with 'ilim1' in the case where log10(d) has been
* overestimated). The numerator and denominator in David Gay's
* conversion algorithm are known to fit in Tcl_WideUInt, giving
* considerably faster arithmetic than mp_int's.
*
* Results:
* Returns the string of significant decimal digits, in newly allocated
* memory
*
* Side effects:
* Stores the location of the decimal point in '*decpt' and the location
* of the terminal null byte in '*endPtr'.
*
*----------------------------------------------------------------------
*/
static inline char *
StrictInt64Conversion(
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 = (char *)Tcl_Alloc(len + 1);
/* Output buffer. */
Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
/* Numerator of the fraction being
* converted. */
Tcl_WideUInt S = wuipow5[s5] << s2;
/* Denominator of the fraction being
* converted. */
int digit; /* Current output digit. */
char *s = retval; /* Cursor in the output buffer. */
int i; /* Current position in the output buffer. */
/*
* Adjust if the logarithm was guessed wrong.
*/
if (b < S) {
b = 10 * b;
ilim = ilim1;
--k;
}
/*
* Loop through the digits.
*/
i = 1;
for (;;) {
digit = (int)(b / S);
if (digit > 10) {
Tcl_Panic("wrong digit!");
}
b = b % S;
/*
* Have we converted all the requested digits?
*/
*s++ = '0' + digit;
if (i == ilim) {
if (2*b > S || (2*b == S && (digit & 1) != 0)) {
s = BumpUp(s, retval, &k);
} else {
while (*--s == '0') {
/* do nothing */
}
++s;
}
break;
}
/*
* Advance to the next digit.
*/
b = 10 * b;
++i;
}
/*
* Endgame - store the location of the decimal point and the end of the
* string.
*/
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* ShouldBankerRoundUpPowD --
*
* Test whether bankers' rounding should round a digit up. Assumption is
* made that the denominator of the fraction being tested is a power of
* 2**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;
}
/*
*----------------------------------------------------------------------
*
* 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 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)
*/
if ((mp_add(b, m, temp) != MP_OKAY) || (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;
}
}
return isodd;
}
/*
*----------------------------------------------------------------------
*
* ShorteningBignumConversionPowD --
*
* Converts a double-precision number to the shortest string of digits
* that reconverts exactly to the given number, or to 'ilim' digits if
* that will yield a shorter result. The denominator in David Gay's
* conversion algorithm is known to be a power of 2**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. */
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 = (char *)Tcl_Alloc(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;
mp_err err = MP_OKAY;
/*
* b = bw * 2**b2 * 5**b5
* mminus = 5**m5
*/
if ((retval == NULL) || (mp_init_u64(&b, bw) != MP_OKAY)) {
return NULL;
}
if (mp_init_set(&mminus, 1) != MP_OKAY) {
mp_clear(&b);
return NULL;
}
err = MulPow5(&b, b5, &b);
if (err == MP_OKAY) {
err = mp_mul_2d(&b, b2, &b);
}
/*
* Adjust if the logarithm was guessed wrong.
*/
if ((err == MP_OKAY) && (b.used <= sd)) {
err = mp_mul_d(&b, 10, &b);
++m2plus; ++m2minus; ++m5;
ilim = ilim1;
--k;
}
/*
* mminus = 5**m5 * 2**m2minus
* mplus = 5**m5 * 2**m2plus
*/
if (err == MP_OKAY) {
err = mp_mul_2d(&mminus, m2minus, &mminus);
}
if (err == MP_OKAY) {
err = MulPow5(&mminus, m5, &mminus);
}
if ((err == MP_OKAY) && (m2plus > m2minus)) {
err = mp_init_copy(&mplus, &mminus);
if (err == MP_OKAY) {
err = mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
}
}
if (err == MP_OKAY) {
err = 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 roundoff of being exact?
*/
r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
if (r1 == MP_LT || (r1 == MP_EQ
&& (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,
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.
*/
if (err == MP_OKAY) {
err = mp_mul_d(&b, 10, &b);
}
if (err == MP_OKAY) {
err = mp_mul_d(&mminus, 10, &mminus);
}
if ((err == MP_OKAY) && (m2plus > m2minus)) {
err = 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, (void *)NULL);
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return (err == MP_OKAY) ? retval : NULL;
}
/*
*----------------------------------------------------------------------
*
* 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(
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 = (char *)Tcl_Alloc(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_err err;
/*
* b = bw * 2**b2 * 5**b5
*/
if (mp_init_u64(&b, bw) != MP_OKAY) {
return NULL;
}
err = MulPow5(&b, b5, &b);
if (err == MP_OKAY) {
err = mp_mul_2d(&b, b2, &b);
}
/*
* Adjust if the logarithm was guessed wrong.
*/
if ((err == MP_OKAY) && (b.used <= sd)) {
err = mp_mul_d(&b, 10, &b);
ilim = ilim1;
--k;
}
/*
* Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
* by mp_digit extraction.
*/
i = 1;
while (err == MP_OKAY) {
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.
*/
err = mp_mul_d(&b, 10, &b);
++i;
}
/*
* Endgame - store the location of the decimal point and the end of the
* string.
*/
mp_clear(&b);
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* ShouldBankerRoundUp --
*
* Tests whether a digit should be rounded up or down when finishing
* bignum-based floating point conversion.
*
* Results:
* Returns 1 if the number needs to be rounded up, 0 otherwise.
*
*----------------------------------------------------------------------
*/
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_EQ:
return isodd;
case MP_GT:
return 1;
default:
return 0;
}
}
/*
*----------------------------------------------------------------------
*
* ShouldBankerRoundUpToNext --
*
* Tests whether the remainder is great enough to force rounding to the
* next higher digit.
*
* Results:
* Returns 1 if the number should be rounded up, 0 otherwise.
*
*----------------------------------------------------------------------
*/
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 isodd) /* 1 if the integer significand is odd. */
{
int r;
mp_int temp;
/*
* Compare b and S-m: this is the same as comparing B+m and S.
*/
if ((mp_init(&temp) != MP_OKAY) || (mp_add(b, m, &temp) != MP_OKAY)) {
return 0;
}
r = mp_cmp_mag(&temp, S);
mp_clear(&temp);
switch (r) {
case MP_EQ:
return isodd;
case MP_GT:
return 1;
default:
return 0;
}
}
/*
*----------------------------------------------------------------------
*
* ShorteningBignumConversion --
*
* Convert a floating point number to a variable-length digit string
* using the multiprecision method.
*
* Results:
* Returns the string of digits.
*
* Side effects:
* Stores the position of the decimal point in *decpt. Stores a pointer
* to the end of the number in *endPtr.
*
*----------------------------------------------------------------------
*/
static inline char *
ShorteningBignumConversion(
Double *dPtr, /* Original number being converted. */
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 = (char *)Tcl_Alloc(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. */
int minit = 1; /* Fudge factor for when we misguess k. */
int i;
int r1;
mp_err err;
/*
* b = bw * 2**b2 * 5**b5
* S = 2**s2 * 5*s5
*/
if ((retval == NULL) || (mp_init_u64(&b, bw) != MP_OKAY)) {
return NULL;
}
err = mp_mul_2d(&b, b2, &b);
if (err == MP_OKAY) {
err = mp_init_set(&S, 1);
}
if (err == MP_OKAY) {
err = MulPow5(&S, s5, &S);
}
if (err == MP_OKAY) {
err = mp_mul_2d(&S, s2, &S);
}
/*
* Handle the case where we guess the position of the decimal point wrong.
*/
if ((err == MP_OKAY) && (mp_cmp_mag(&b, &S) == MP_LT)) {
err = mp_mul_d(&b, 10, &b);
minit = 10;
ilim =ilim1;
--k;
}
/*
* mminus = 2**m2minus * 5**m5
*/
if (err == MP_OKAY) {
err = mp_init_set(&mminus, minit);
}
if (err == MP_OKAY) {
err = mp_mul_2d(&mminus, m2minus, &mminus);
}
if ((err == MP_OKAY) && (m2plus > m2minus)) {
err = mp_init_copy(&mplus, &mminus);
if (err == MP_OKAY) {
err = mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
}
}
/*
* Loop through the digits.
*/
if (err == MP_OKAY) {
err = mp_init(&dig);
}
i = 1;
while (err == MP_OKAY) {
err = 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 && (dPtr->w.word1 & 1) == 0)) {
err = 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,
dPtr->w.word1 & 1)) {
++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 ((err == MP_OKAY) && (i == ilim)) {
err = mp_mul_2d(&b, 1, &b);
if (ShouldBankerRoundUp(&b, &S, digit&1)) {
s = BumpUp(s, retval, &k);
}
break;
}
/*
* Advance to the next digit.
*/
if ((err == MP_OKAY) && (s5 > 0)) {
/*
* Can possibly shorten the denominator.
*/
err = mp_mul_2d(&b, 1, &b);
if (err == MP_OKAY) {
err = mp_mul_2d(&mminus, 1, &mminus);
}
if ((err == MP_OKAY) && (m2plus > m2minus)) {
err = mp_mul_2d(&mplus, 1, &mplus);
}
if (err == MP_OKAY) {
err = 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 if (err == MP_OKAY) {
err = mp_mul_d(&b, 10, &b);
if (err == MP_OKAY) {
err = mp_mul_d(&mminus, 10, &mminus);
}
if ((err == MP_OKAY) && (m2plus > m2minus)) {
err = 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, &dig, &S, (void *)NULL);
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* StrictBignumConversion --
*
* Convert a floating point number to a fixed-length digit string using
* the multiprecision method.
*
* Results:
* Returns the string of digits.
*
* Side effects:
* Stores the position of the decimal point in *decpt. Stores a pointer
* to the end of the number in *endPtr.
*
*----------------------------------------------------------------------
*/
static inline char *
StrictBignumConversion(
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 = (char *)Tcl_Alloc(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. */
int g; /* Size of the current digit ground. */
int i, j;
mp_err err;
/*
* b = bw * 2**b2 * 5**b5
* S = 2**s2 * 5*s5
*/
if (mp_init(&dig) != MP_OKAY) {
return NULL;
}
if (mp_init_u64(&b, bw) != MP_OKAY) {
mp_clear(&dig);
return NULL;
}
err = mp_mul_2d(&b, b2, &b);
if (err == MP_OKAY) {
err = mp_init_set(&S, 1);
}
if (err == MP_OKAY) {
err = MulPow5(&S, s5, &S);
if (err == MP_OKAY) {
err = 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) == MP_OKAY)) {
ilim =ilim1;
--k;
}
/*
* Convert the leading digit.
*/
i = 0;
err = 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) {
if ((mp_mul_2d(&b, 1, &b) == MP_OKAY) && ShouldBankerRoundUp(&b, &S, digit&1)) {
s = BumpUp(s, retval, &k);
}
} else {
while (err == MP_OKAY) {
/*
* Shift by a group of digits.
*/
g = ilim - i;
if (g > DIGIT_GROUP) {
g = DIGIT_GROUP;
}
if (s5 >= g) {
err = mp_div_d(&S, dpow5[g], &S, NULL);
s5 -= g;
} else if (s5 > 0) {
err = mp_div_d(&S, dpow5[s5], &S, NULL);
if (err == MP_OKAY) {
err = mp_mul_d(&b, dpow5[g - s5], &b);
}
s5 = 0;
} else {
err = mp_mul_d(&b, dpow5[g], &b);
}
if (err == MP_OKAY) {
err = 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.
*/
if ((err != MP_OKAY) || (mp_div(&b, &S, &dig, &b) != MP_OKAY) || (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) {
if ((mp_mul_2d(&b, 1, &b) == MP_OKAY) && 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, &dig, (void *)NULL);
*s = '\0';
*decpt = k;
if (endPtr) {
*endPtr = s;
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* TclDoubleDigits --
*
* Core of Tcl's conversion of double-precision floating point numbers to
* decimal.
*
* Results:
* Returns a newly-allocated string of digits.
*
* Side effects:
* Sets *decpt to the index of the character in the string before the
* place that the decimal point should go. If 'endPtr' is not NULL, sets
* endPtr to point to the terminating '\0' byte of the string. Sets *sign
* to 1 if a minus sign should be printed with the number, or 0 if a plus
* sign (or no sign) should appear.
*
* This function is a service routine that produces the string of digits for
* floating-point-to-decimal conversion. It can do a number of things
* according to the 'flags' argument. Valid values for 'flags' include:
* TCL_DD_SHORTEST - This is the default for floating point conversion.
* 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_E_FORMAT - This value is used to prepare numbers for %e format
* conversion. 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_SHORTEST 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_SHORTEST 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. */
{
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 ndigits for E format.
* ilim = The number of significant digits to convert if k has been
* guessed correctly. This is -1 for shortest (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 'ndigits' for E format, but it's 'ndigits-1' for F
* format.
*/
SetPrecisionLimits(flags, 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_SHORTEST) {
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, 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, 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, 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(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(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(bw, b2, s2, s5, k,
len, ilim, ilim1, decpt, endPtr);
}
}
}
/*
*----------------------------------------------------------------------
*
* TclInitDoubleConversion --
*
* Initializes constants that are needed for conversions to and from
* 'double'
*
* Results:
* None.
*
* Side effects:
* The log base 2 of the floating point radix, the number of bits in a
* double mantissa, and a table of the powers of five and ten are
* computed and stored.
*
*----------------------------------------------------------------------
*/
void
TclInitDoubleConversion(void)
{
int i;
int x;
Tcl_WideUInt u;
double d;
#ifdef IEEE_FLOATING_POINT
union {
double dv;
Tcl_WideUInt iv;
} bitwhack;
#endif
mp_err err = MP_OKAY;
#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 *)
Tcl_Alloc((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) {
err = err || mp_init(pow5 + i);
}
mp_set_u64(pow5, 5);
for (i=0; i<8; ++i) {
err = err || mp_sqr(pow5+i, pow5+i+1);
}
err = err || mp_init_u64(pow5_13, 1220703125);
for (i = 1; i < 5; ++i) {
err = err || mp_init(pow5_13 + i);
err = err || mp_sqr(pow5_13 + i - 1, pow5_13 + i);
}
if (err != MP_OKAY) {
Tcl_Panic("out of memory");
}
/*
* 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
}
/*
*----------------------------------------------------------------------
*
* TclFinalizeDoubleConversion --
*
* Cleans up this file on exit.
*
* Results:
* None
*
* Side effects:
* Memory allocated by TclInitDoubleConversion is freed.
*
*----------------------------------------------------------------------
*/
void
TclFinalizeDoubleConversion(void)
{
int i;
Tcl_Free(pow10_wide);
for (i=0; i<9; ++i) {
mp_clear(pow5 + i);
}
for (i=0; i < 5; ++i) {
mp_clear(pow5_13 + i);
}
}
/*
*----------------------------------------------------------------------
*
* Tcl_InitBignumFromDouble --
*
* Extracts the integer part of a double and converts it to an arbitrary
* precision integer.
*
* Results:
* None.
*
* Side effects:
* Initializes the bignum supplied, and stores the converted number in
* it.
*
*----------------------------------------------------------------------
*/
int
Tcl_InitBignumFromDouble(
Tcl_Interp *interp, /* For error message. */
double d, /* Number to convert. */
void *big) /* Place to store the result. */
{
double fract;
int expt;
mp_err err;
mp_int *b = (mp_int *)big;
/*
* Infinite values can't convert to bignum.
*/
if (isinf(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, (char *)NULL);
}
return TCL_ERROR;
}
fract = frexp(d, &expt);
if (expt <= 0) {
err = mp_init(b);
mp_zero(b);
} else {
Tcl_WideInt w = (Tcl_WideInt)ldexp(fract, mantBits);
int shift = expt - mantBits;
err = mp_init_i64(b, w);
if (err != MP_OKAY) {
/* just skip */
} else if (shift < 0) {
err = mp_div_2d(b, -shift, b, NULL);
} else if (shift > 0) {
err = mp_mul_2d(b, shift, b);
}
}
if (err != MP_OKAY) {
return TCL_ERROR;
}
return TCL_OK;
}
/*
*----------------------------------------------------------------------
*
* TclBignumToDouble --
*
* Convert an arbitrary-precision integer to a native floating point
* number.
*
* Results:
* Returns the converted number. Sets errno to ERANGE if the number is
* too large to convert.
*
*----------------------------------------------------------------------
*/
double
TclBignumToDouble(
const void *big) /* Integer to convert. */
{
mp_int b;
int bits, shift, i, lsb;
double r;
mp_err err;
const mp_int *a = (const mp_int *)big;
/*
* 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'.
*/
err = mp_init(&b);
if (err != MP_OKAY) {
/* just skip */
} else if (shift == 0) {
err = mp_copy(a, &b);
} else if (shift > 0) {
err = mp_mul_2d(a, shift, &b);
} else if (shift < 0) {
lsb = mp_cnt_lsb(a);
if (lsb == -1-shift) {
/*
* Round to even
*/
err = mp_div_2d(a, -shift, &b, NULL);
if ((err == MP_OKAY) && mp_isodd(&b)) {
if (mp_isneg(&b)) {
err = mp_sub_d(&b, 1, &b);
} else {
err = mp_add_d(&b, 1, &b);
}
}
} else {
/*
* Ordinary rounding
*/
err = mp_div_2d(a, -1-shift, &b, NULL);
if (err != MP_OKAY) {
/* just skip */
} else if (mp_isneg(&b)) {
err = mp_sub_d(&b, 1, &b);
} else {
err = mp_add_d(&b, 1, &b);
}
err = mp_div_2d(&b, 1, &b, NULL);
}
}
/*
* Accumulate the result, one mp_digit at a time.
*/
if (err != MP_OKAY) {
return 0.0;
}
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;
}
}
/*
*----------------------------------------------------------------------
*
* TclCeil --
*
* Computes the smallest floating point number that is at least the
* mp_int argument.
*
* Results:
* Returns the floating point number.
*
*----------------------------------------------------------------------
*/
double
TclCeil(
const void *big) /* Integer to convert. */
{
double r = 0.0;
mp_int b;
mp_err err;
const mp_int *a = (const mp_int *)big;
err = mp_init(&b);
if ((err == MP_OKAY) && mp_isneg(a)) {
err = 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 (err != MP_OKAY) {
/* just skip */
} else if (shift > 0) {
err = mp_mul_2d(a, shift, &b);
} else if (shift < 0) {
mp_int d;
err = mp_init(&d);
if (err == MP_OKAY) {
err = mp_div_2d(a, -shift, &b, &d);
}
exact = mp_iszero(&d);
mp_clear(&d);
} else {
err = mp_copy(a, &b);
}
if ((err == MP_OKAY) && !exact) {
err = mp_add_d(&b, 1, &b);
}
if (err != MP_OKAY) {
return 0.0;
}
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;
}
/*
*----------------------------------------------------------------------
*
* 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 void *big) /* Integer to convert. */
{
double r = 0.0;
mp_int b;
mp_err err;
const mp_int *a = (const mp_int *)big;
err = mp_init(&b);
if ((err == MP_OKAY) && mp_isneg(a)) {
err = 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) {
err = mp_mul_2d(a, shift, &b);
} else if (shift < 0) {
err = mp_div_2d(a, -shift, &b, NULL);
} else {
err = mp_copy(a, &b);
}
if (err != MP_OKAY) {
return 0.0;
}
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;
}
/*
*----------------------------------------------------------------------
*
* 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;
mp_err err = MP_OKAY;
/*
* 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;
if (mp_init(&b)) {
return 0.0;
}
if (shift > 0) {
err = mp_mul_2d(a, shift, &b);
} else if (shift < 0) {
err = mp_div_2d(a, -shift, &b, NULL);
} else {
err = mp_copy(a, &b);
}
/*
* Accumulate the result, one mp_digit at a time.
*/
r = 0.0;
if (err == MP_OKAY) {
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);
}
/*
*----------------------------------------------------------------------
*
* Pow10TimesFrExp --
*
* Multiply a power of ten by a number expressed as fraction and
* exponent.
*
* Results:
* Returns the significand of the result.
*
* Side effects:
* Overwrites the 'machexp' parameter with the exponent of the result.
*
* Assumes that 'exponent' is such that 10**exponent would be a double, even
* though 'fraction*10**(machexp+exponent)' might overflow.
*
*----------------------------------------------------------------------
*/
static double
Pow10TimesFrExp(
int exponent, /* Power of 10 to multiply by. */
double fraction, /* Significand of multiplicand. */
int *machexp) /* On input, exponent of multiplicand. On
* output, exponent of result. */
{
int i, j;
int expt = *machexp;
double retval = fraction;
if (exponent > 0) {
/*
* Multiply by 10**exponent.
*/
retval = frexp(retval * pow10vals[exponent & 0xF], &j);
expt += j;
for (i=4; i<9; ++i) {
if (exponent & (1<<i)) {
retval = frexp(retval * pow_10_2_n[i], &j);
expt += j;
}
}
} else if (exponent < 0) {
/*
* Divide by 10**-exponent.
*/
retval = frexp(retval / pow10vals[(-exponent) & 0xF], &j);
expt += j;
for (i=4; i<9; ++i) {
if ((-exponent) & (1<<i)) {
retval = frexp(retval / pow_10_2_n[i], &j);
expt += j;
}
}
}
*machexp = expt;
return retval;
}
/*
*----------------------------------------------------------------------
*
* SafeLdExp --
*
* Do an 'ldexp' operation, but handle denormals gracefully.
*
* Results:
* Returns the appropriately scaled value.
*
* On some platforms, 'ldexp' fails when presented with a number too
* small to represent as a normalized double. This routine does 'ldexp'
* in two steps for those numbers, to return correctly denormalized
* values.
*
*----------------------------------------------------------------------
*/
static double
SafeLdExp(
double fract,
int expt)
{
int minexpt = DBL_MIN_EXP * log2FLT_RADIX;
volatile double a, b, retval;
if (expt < minexpt) {
a = ldexp(fract, expt - mantBits - minexpt);
b = ldexp(1.0, mantBits + minexpt);
retval = a * b;
} else {
retval = ldexp(fract, expt);
}
return retval;
}
/*
*----------------------------------------------------------------------
*
* TclFormatNaN --
*
* Makes the string representation of a "Not a Number"
*
* Results:
* None.
*
* Side effects:
* Stores the string representation in the supplied buffer, which must be
* at least TCL_DOUBLE_SPACE characters.
*
*----------------------------------------------------------------------
*/
void
TclFormatNaN(
double value, /* The Not-a-Number to format. */
char *buffer) /* String representation. */
{
#ifndef IEEE_FLOATING_POINT
strcpy(buffer, "NaN");
return;
#else
union {
double dv;
uint64_t iv;
} bitwhack;
bitwhack.dv = value;
if (n770_fp) {
bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
}
if (bitwhack.iv & (UINT64_C(1) << 63)) {
bitwhack.iv &= ~ (UINT64_C(1) << 63);
*buffer++ = '-';
}
*buffer++ = 'N';
*buffer++ = 'a';
*buffer++ = 'N';
bitwhack.iv &= ((UINT64_C(1)) << 51) - 1;
if (bitwhack.iv != 0) {
snprintf(buffer, TCL_DOUBLE_SPACE, "(%" PRIx64 ")", bitwhack.iv);
} else {
*buffer = '\0';
}
#endif /* IEEE_FLOATING_POINT */
}
/*
*----------------------------------------------------------------------
*
* Nokia770Twiddle --
*
* Transpose the two words of a number for Nokia 770 floating point
* handling.
*
*----------------------------------------------------------------------
*/
#ifdef IEEE_FLOATING_POINT
static Tcl_WideUInt
Nokia770Twiddle(
Tcl_WideUInt w) /* Number to transpose. */
{
return (((w >> 32) & 0xFFFFFFFF) | (w << 32));
}
#endif
/*
*----------------------------------------------------------------------
*
* TclNokia770Doubles --
*
* Transpose the two words of a number for Nokia 770 floating point
* handling.
*
*----------------------------------------------------------------------
*/
int
TclNokia770Doubles(void)
{
return n770_fp;
}
/*
* Local Variables:
* mode: c
* c-basic-offset: 4
* fill-column: 78
* End:
*/
