diff options
Diffstat (limited to 'generic/tclStrToD.c')
| -rw-r--r-- | generic/tclStrToD.c | 3397 |
1 files changed, 1507 insertions, 1890 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c index 09fd1f3..cff9bdd 100644 --- a/generic/tclStrToD.c +++ b/generic/tclStrToD.c @@ -1,4 +1,6 @@ /* + *---------------------------------------------------------------------- + * * tclStrToD.c -- * * This file contains a collection of procedures for managing conversions @@ -7,25 +9,23 @@ * into strings of digits, and procedures for interconversion among * 'double' and 'mp_int' types. * - * Copyright © 2005 Kevin B. Kenny. All rights reserved. + * Copyright (c) 2005 by Kevin B. Kenny. All rights reserved. * * See the file "license.terms" for information on usage and redistribution of * this file, and for a DISCLAIMER OF ALL WARRANTIES. + *---------------------------------------------------------------------- */ #include "tclInt.h" -#include "tclTomMath.h" -#include <float.h> +#include "tommath.h" #include <math.h> -#ifdef _WIN32 -#define copysign _copysign -#endif - -#ifndef PRIx64 -# define PRIx64 TCL_LL_MODIFIER "x" -#endif +/* + * Define KILL_OCTAL to suppress interpretation of numbers with leading zero + * as octal. (Ceterum censeo: numeros octonarios delendos esse.) + */ +#undef KILL_OCTAL /* * This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754 @@ -38,76 +38,48 @@ #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 +#if defined(__GNUC__) && defined(__i386) +typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); +#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw)) +#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw)) +# define FPU_IEEE_ROUNDING 0x027f +# define ADJUST_FPU_CONTROL_WORD #endif -/* - * 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. + +/* Sun ProC needs sunmath for rounding control on x86 like gcc above. + * + * */ -#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) +#if defined(__sun) && defined(__i386) && !defined(__GNUC__) +#include <sunmath.h> #endif /* - * MIPS floating-point units need special settings in control registers to use - * gradual underflow as we expect. This fix is for the MIPSpro compiler. + * 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) +# define NAN_START 0x7ff4 +# define NAN_MASK (((Tcl_WideUInt) 1) << 50) #else -# define NAN_START 0x7FF8 -# define NAN_MASK (((Tcl_WideUInt) 1) << 51) +# define NAN_START 0x7ff8 +# define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif /* @@ -121,44 +93,45 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); #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 +/* 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) + * 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. - */ +#define FP_PRECISION 53 + /* Number of bits of significand plus the + * hidden bit */ +#define EXPONENT_BIAS 0x3ff + /* Bias of the exponent 0 */ + +/* Derived quantities */ + +#define TEN_PMAX 22 + /* floor(FP_PRECISION*log(2)/log(5)) */ +#define QUICK_MAX 14 + /* floor((FP_PRECISION-1)*log(2)/log(10)) - 1 */ +#define BLETCH 0x10 + /* Highest power of two that is greater than + * DBL_MAX_10_EXP, divided by 16 */ +#define DIGIT_GROUP 8 + /* floor(DIGIT_BIT*log(2)/log(10)) */ + +/* Union used to dismantle floating point numbers. */ typedef union Double { struct { @@ -189,7 +162,7 @@ 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 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 @@ -211,12 +184,10 @@ static int n770_fp; /* Flag is 1 on Nokia N770 floating point. * 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. */ + * 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. - */ +/* 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, @@ -225,10 +196,7 @@ static const mp_digit dpow5[13] = { 244140625 }; -/* - * Table of powers: pow5_13[n] = 5**(13*2**(n+1)) - */ - +/* 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[] = { @@ -261,35 +229,34 @@ static const int log2pow5[27] = { }; #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 const Tcl_WideUInt wuipow5[27] = { + (Tcl_WideUInt) 1, /* 5**0 */ + (Tcl_WideUInt) 5, + (Tcl_WideUInt) 25, + (Tcl_WideUInt) 125, + (Tcl_WideUInt) 625, + (Tcl_WideUInt) 3125, /* 5**5 */ + (Tcl_WideUInt) 3125*5, + (Tcl_WideUInt) 3125*25, + (Tcl_WideUInt) 3125*125, + (Tcl_WideUInt) 3125*625, + (Tcl_WideUInt) 3125*3125, /* 5**10 */ + (Tcl_WideUInt) 3125*3125*5, + (Tcl_WideUInt) 3125*3125*25, + (Tcl_WideUInt) 3125*3125*125, + (Tcl_WideUInt) 3125*3125*625, + (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */ + (Tcl_WideUInt) 3125*3125*3125*5, + (Tcl_WideUInt) 3125*3125*3125*25, + (Tcl_WideUInt) 3125*3125*3125*125, + (Tcl_WideUInt) 3125*3125*3125*625, + (Tcl_WideUInt) 3125*3125*3125*3125, /* 5**20 */ + (Tcl_WideUInt) 3125*3125*3125*3125*5, + (Tcl_WideUInt) 3125*3125*3125*3125*25, + (Tcl_WideUInt) 3125*3125*3125*3125*125, + (Tcl_WideUInt) 3125*3125*3125*3125*625, + (Tcl_WideUInt) 3125*3125*3125*3125*3125, /* 5**25 */ + (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */ }; /* @@ -299,78 +266,72 @@ static const Tcl_WideUInt wuipow5[] = { static int AccumulateDecimalDigit(unsigned, int, Tcl_WideUInt *, mp_int *, int); static double MakeHighPrecisionDouble(int signum, - mp_int *significand, int nSigDigs, long exponent); + mp_int *significand, int nSigDigs, int exponent); static double MakeLowPrecisionDouble(int signum, Tcl_WideUInt significand, int nSigDigs, - long exponent); -#ifdef IEEE_FLOATING_POINT + int exponent); 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 void MulPow5(mp_int*, unsigned, mp_int*); +static int NormalizeRightward(Tcl_WideUInt*); static int RequiredPrecision(Tcl_WideUInt); -static void DoubleToExpAndSig(double, Tcl_WideUInt *, int *, - int *); -static void TakeAbsoluteValue(Double *, int *); -static char * FormatInfAndNaN(Double *, int *, char **); -static char * FormatZero(int *, char **); +static 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, +static int BetterLog10(double, int, int*); +static void ComputeScale(int, int, int*, int*, int*, int*); +static void SetPrecisionLimits(int, int, int*, int*, int*, int*); +static char* BumpUp(char*, char*, int*); +static int AdjustRange(double*, int); +static char* ShorteningQuickFormat(double, int, int, double, + char*, int*); +static char* StrictQuickFormat(double, int, int, double, + char*, int*); +static char* QuickConversion(double, int, int, int, int, int, int, + int*, char**); +static void CastOutPowersOf2(int*, int*, int*); +static char* ShorteningInt64Conversion(Double*, int, Tcl_WideUInt, int, int, int, int, int, int, int, int, int, - int, int, int *, char **); -static char * StrictInt64Conversion(Tcl_WideUInt, + int, int, int*, char**); +static char* StrictInt64Conversion(Double*, int, Tcl_WideUInt, int, int, int, int, int, int, - int, int, int *, char **); -static int ShouldBankerRoundUpPowD(mp_int *, int, int); -static int ShouldBankerRoundUpToNextPowD(mp_int *, mp_int *, - int, int, mp_int *); -static char * ShorteningBignumConversionPowD(Double *dPtr, - Tcl_WideUInt bw, int b2, int b5, + int, int, int*, char**); +static int ShouldBankerRoundUpPowD(mp_int*, int, int); +static int ShouldBankerRoundUpToNextPowD(mp_int*, mp_int*, + int, int, int, mp_int*); +static char* ShorteningBignumConversionPowD(Double* dPtr, + int convType, Tcl_WideUInt bw, int b2, int b5, int m2plus, int m2minus, int m5, int sd, int k, int len, - int ilim, int ilim1, int *decpt, - char **endPtr); -static char * StrictBignumConversionPowD( + int ilim, int ilim1, int* decpt, + char** endPtr); +static char* StrictBignumConversionPowD(Double* dPtr, int convType, Tcl_WideUInt bw, int b2, int b5, int sd, int k, int len, - int ilim, int ilim1, int *decpt, - char **endPtr); -static int ShouldBankerRoundUp(mp_int *, mp_int *, int); -static int ShouldBankerRoundUpToNext(mp_int *, mp_int *, - mp_int *, int); -static char * ShorteningBignumConversion(Double *dPtr, + int ilim, int ilim1, int* decpt, + char** endPtr); +static int ShouldBankerRoundUp(mp_int*, mp_int*, int); +static int ShouldBankerRoundUpToNext(mp_int*, mp_int*, mp_int*, + int, int, mp_int*); +static char* ShorteningBignumConversion(Double* dPtr, int convType, Tcl_WideUInt bw, int b2, int m2plus, int m2minus, int s2, int s5, int k, int len, - int ilim, int ilim1, int *decpt, - char **endPtr); -static char * StrictBignumConversion( + int ilim, int ilim1, int* decpt, + char** endPtr); +static char* StrictBignumConversion(Double* dPtr, int convType, Tcl_WideUInt bw, int b2, int s2, int s5, int k, int len, - int ilim, int ilim1, int *decpt, - char **endPtr); -static double BignumToBiasedFrExp(const mp_int *big, int *machexp); + int ilim, int ilim1, int* decpt, + char** endPtr); +static double BignumToBiasedFrExp(mp_int *big, int *machexp); static double Pow10TimesFrExp(int exponent, double fraction, int *machexp); static double SafeLdExp(double fraction, int exponent); -#ifdef IEEE_FLOATING_POINT static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w); -#endif /* *---------------------------------------------------------------------- @@ -399,17 +360,14 @@ static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w); * - 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 + * - 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, + * - 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 + * - 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 @@ -492,7 +450,7 @@ TclParseNumber( { enum State { INITIAL, SIGNUM, ZERO, ZERO_X, - ZERO_O, ZERO_B, ZERO_D, BINARY, + ZERO_O, ZERO_B, BINARY, HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL, LEADING_RADIX_POINT, FRACTION, EXPONENT_START, EXPONENT_SIGNUM, EXPONENT, @@ -503,45 +461,45 @@ TclParseNumber( } state = INITIAL; enum State acceptState = INITIAL; - int signum = 0; /* Sign of the number being parsed. */ + int signum = 0; /* Sign of the number being parsed */ Tcl_WideUInt significandWide = 0; /* Significand of the number being parsed (if - * no overflow). */ + * no overflow) */ mp_int significandBig; /* Significand of the number being parsed (if - * it overflows significandWide). */ - int significandOverflow = 0;/* Flag==1 iff significandBig is used. */ + * 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'. */ + * we've scanned forward to a '.' or 'e' */ mp_int octalSignificandBig; /* Significand of octal number once - * octalSignificandWide overflows. */ + * octalSignificandWide overflows */ int octalSignificandOverflow = 0; - /* Flag==1 if octalSignificandBig is used. */ + /* Flag==1 if octalSignificandBig is used */ int numSigDigs = 0; /* Number of significant digits in the decimal - * significand. */ + * 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. */ + * point */ int exponentSignum = 0; /* Signum of the exponent of a floating point - * number. */ - long exponent = 0; /* Exponent of a floating point number. */ - const char *p; /* Pointer to next character to scan. */ - size_t len; /* Number of characters remaining after p. */ + * number */ + long exponent = 0; /* Exponent of a floating point number */ + const char *p; /* Pointer to next character to scan */ + size_t len; /* Number of characters remaining after p */ const char *acceptPoint; /* Pointer to position after last character in - * an acceptable number. */ + * an acceptable number */ size_t acceptLen; /* Number of characters following that * point. */ - int status = TCL_OK; /* Status to return to caller. */ + int status = TCL_OK; /* Status to return to caller */ char d = 0; /* Last hexadecimal digit scanned; initialized * to avoid a compiler warning. */ int shift = 0; /* Amount to shift when accumulating binary */ int explicitOctal = 0; - mp_err err = MP_OKAY; -#define MOST_BITS (UWIDE_MAX >> 1) +#define ALL_BITS (~(Tcl_WideUInt)0) +#define MOST_BITS (ALL_BITS >> 1) /* * Initialize bytes to start of the object's string rep if the caller @@ -549,20 +507,6 @@ TclParseNumber( */ 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)) { - int 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); } @@ -572,87 +516,6 @@ TclParseNumber( 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: @@ -661,7 +524,7 @@ TclParseNumber( * I, N, and whitespace. */ - if (TclIsSpaceProcM(c)) { + if (TclIsSpaceProc(c)) { if (flags & TCL_PARSE_NO_WHITESPACE) { goto endgame; } @@ -691,8 +554,6 @@ TclParseNumber( 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))) { @@ -719,18 +580,15 @@ TclParseNumber( 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'. + * 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; } @@ -741,25 +599,15 @@ TclParseNumber( goto zeroo; } if (c == 'b' || c == 'B') { - if ((flags & TCL_PARSE_OCTAL_ONLY)) { - goto endgame; - } state = ZERO_B; break; } - if (flags & TCL_PARSE_BINARY_ONLY) { - goto zerob; - } if (c == 'o' || c == 'O') { explicitOctal = 1; state = ZERO_O; break; } - if (c == 'd' || c == 'D') { - state = ZERO_D; - break; - } -#ifdef TCL_NO_DEPRECATED +#ifdef KILL_OCTAL goto decimal; #endif /* FALLTHROUGH */ @@ -791,45 +639,29 @@ TclParseNumber( 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. + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check for + * too large shifts first. */ if ((octalSignificandWide != 0) && (((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt)) || (octalSignificandWide > - (UWIDE_MAX >> shift)))) { + (~(Tcl_WideUInt)0 >> shift)))) { octalSignificandOverflow = 1; - err = mp_init_u64(&octalSignificandBig, + TclBNInitBignumFromWideUInt(&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'; + octalSignificandWide = + (octalSignificandWide << shift) + (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; + mp_mul_2d(&octalSignificandBig, shift, + &octalSignificandBig); + mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'), + &octalSignificandBig); } } if (numSigDigs != 0) { @@ -858,7 +690,7 @@ TclParseNumber( goto endgame; } -#ifndef TCL_NO_DEPRECATED +#ifndef KILL_OCTAL /* * Scanned a number with a leading zero that contains an 8, 9, @@ -926,41 +758,26 @@ TclParseNumber( 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. + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check for too + * large shifts first. */ if (significandWide != 0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || - significandWide > (UWIDE_MAX >> shift))) { + significandWide > (~(Tcl_WideUInt)0 >> shift))) { significandOverflow = 1; - err = mp_init_u64(&significandBig, + TclBNInitBignumFromWideUInt(&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); - } + significandWide = (significandWide << shift) + d; + } else { + mp_mul_2d(&significandBig, shift, &significandBig); + mp_add_d(&significandBig, (mp_digit) d, &significandBig); } } - if (err != MP_OKAY) { - return TCL_ERROR; - } numTrailZeros = 0; state = HEXADECIMAL; break; @@ -969,9 +786,7 @@ TclParseNumber( acceptState = state; acceptPoint = p; acceptLen = len; - /* FALLTHRU */ case ZERO_B: - zerob: if (c == '0') { numTrailZeros++; state = BINARY; @@ -983,62 +798,37 @@ TclParseNumber( 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. + * Shifting by more bits than are in the value being + * shifted is at least de facto nonportable. Check for too + * large shifts first. */ if (significandWide != 0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || - significandWide > (UWIDE_MAX >> shift))) { + significandWide > (~(Tcl_WideUInt)0 >> shift))) { significandOverflow = 1; - err = mp_init_u64(&significandBig, + TclBNInitBignumFromWideUInt(&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); - } + significandWide = (significandWide << shift) + 1; + } else { + mp_mul_2d(&significandBig, shift, &significandBig); + 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. */ -#ifdef TCL_NO_DEPRECATED +#ifdef KILL_OCTAL decimal: #endif acceptState = state; @@ -1248,7 +1038,7 @@ TclParseNumber( } /* FALLTHROUGH */ case sNANPAREN: - if (TclIsSpaceProcM(c)) { + if (TclIsSpaceProc(c)) { break; } if (numSigDigs < 13) { @@ -1276,7 +1066,6 @@ TclParseNumber( acceptLen = len; goto endgame; } - next: p++; len--; } @@ -1294,19 +1083,16 @@ TclParseNumber( } 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)) { + while (len != 0 && TclIsSpaceProc(*p)) { p++; len--; } @@ -1325,14 +1111,13 @@ TclParseNumber( */ if (status == TCL_OK && objPtr != NULL) { - TclFreeInternalRep(objPtr); + TclFreeIntRep(objPtr); switch (acceptState) { case SIGNUM: case BAD_OCTAL: case ZERO_X: case ZERO_O: case ZERO_B: - case ZERO_D: case LEADING_RADIX_POINT: case EXPONENT_START: case EXPONENT_SIGNUM: @@ -1347,35 +1132,24 @@ TclParseNumber( case sNA: case sNANPAREN: case sNANHEX: -#endif Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'", acceptState, bytes); +#endif case BINARY: shift = numTrailZeros; if (!significandOverflow && significandWide != 0 && ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || significandWide > (MOST_BITS + signum) >> shift)) { significandOverflow = 1; - err = mp_init_u64(&significandBig, significandWide); + TclBNInitBignumFromWideUInt(&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); + significandWide <<= shift; + } else { + mp_mul_2d(&significandBig, shift, &significandBig); } } - if (err != MP_OKAY) { - return TCL_ERROR; - } goto returnInteger; case HEXADECIMAL: @@ -1388,31 +1162,20 @@ TclParseNumber( ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || significandWide > (MOST_BITS + signum) >> shift)) { significandOverflow = 1; - err = mp_init_u64(&significandBig, significandWide); + TclBNInitBignumFromWideUInt(&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); + significandWide <<= shift; + } else { + mp_mul_2d(&significandBig, shift, &significandBig); } } - if (err != MP_OKAY) { - return TCL_ERROR; - } goto returnInteger; case OCTAL: /* - * Returning an octal integer. Final scaling step. + * Returning an octal integer. Final scaling step */ shift = 3 * numTrailZeros; @@ -1420,49 +1183,52 @@ TclParseNumber( ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) || octalSignificandWide > (MOST_BITS + signum) >> shift)) { octalSignificandOverflow = 1; - err = mp_init_u64(&octalSignificandBig, + TclBNInitBignumFromWideUInt(&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, + octalSignificandWide <<= shift; + } else { + mp_mul_2d(&octalSignificandBig, shift, &octalSignificandBig); } } if (!octalSignificandOverflow) { - if ((err == MP_OKAY) && (octalSignificandWide > (MOST_BITS + signum))) { - err = mp_init_u64(&octalSignificandBig, + if (octalSignificandWide > + (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { +#ifndef NO_WIDE_TYPE + if (octalSignificandWide <= (MOST_BITS + signum)) { + objPtr->typePtr = &tclWideIntType; + if (signum) { + objPtr->internalRep.wideValue = + - (Tcl_WideInt) octalSignificandWide; + } else { + objPtr->internalRep.wideValue = + (Tcl_WideInt) octalSignificandWide; + } + break; + } +#endif + TclBNInitBignumFromWideUInt(&octalSignificandBig, octalSignificandWide); octalSignificandOverflow = 1; } else { objPtr->typePtr = &tclIntType; if (signum) { - objPtr->internalRep.wideValue = - (Tcl_WideInt)(-octalSignificandWide); + objPtr->internalRep.longValue = + - (long) octalSignificandWide; } else { - objPtr->internalRep.wideValue = - (Tcl_WideInt)octalSignificandWide; + objPtr->internalRep.longValue = + (long) octalSignificandWide; } } } - if ((err == MP_OKAY) && octalSignificandOverflow) { + if (octalSignificandOverflow) { if (signum) { - err = mp_neg(&octalSignificandBig, &octalSignificandBig); + mp_neg(&octalSignificandBig, &octalSignificandBig); } - TclSetBignumInternalRep(objPtr, &octalSignificandBig); - } - if (err != MP_OKAY) { - return TCL_ERROR; + TclSetBignumIntRep(objPtr, &octalSignificandBig); } break; @@ -1470,35 +1236,46 @@ TclParseNumber( case DECIMAL: significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1, &significandWide, &significandBig, significandOverflow); - if ((err == MP_OKAY) && !significandOverflow && (significandWide > MOST_BITS+signum)) { + if (!significandOverflow && (significandWide > MOST_BITS+signum)) { significandOverflow = 1; - err = mp_init_u64(&significandBig, significandWide); + TclBNInitBignumFromWideUInt(&significandBig, significandWide); } returnInteger: if (!significandOverflow) { - if ((err == MP_OKAY) && (significandWide > MOST_BITS+signum)) { - err = mp_init_u64(&significandBig, + if (significandWide > + (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) { +#ifndef NO_WIDE_TYPE + if (significandWide <= MOST_BITS+signum) { + objPtr->typePtr = &tclWideIntType; + if (signum) { + objPtr->internalRep.wideValue = + - (Tcl_WideInt) significandWide; + } else { + objPtr->internalRep.wideValue = + (Tcl_WideInt) significandWide; + } + break; + } +#endif + TclBNInitBignumFromWideUInt(&significandBig, significandWide); significandOverflow = 1; } else { objPtr->typePtr = &tclIntType; if (signum) { - objPtr->internalRep.wideValue = - (Tcl_WideInt)(-significandWide); + objPtr->internalRep.longValue = + - (long) significandWide; } else { - objPtr->internalRep.wideValue = - (Tcl_WideInt)significandWide; + objPtr->internalRep.longValue = + (long) significandWide; } } } - if ((err == MP_OKAY) && significandOverflow) { + if (significandOverflow) { if (signum) { - err = mp_neg(&significandBig, &significandBig); + mp_neg(&significandBig, &significandBig); } - TclSetBignumInternalRep(objPtr, &significandBig); - } - if (err != MP_OKAY) { - return TCL_ERROR; + TclSetBignumIntRep(objPtr, &significandBig); } break; @@ -1515,45 +1292,16 @@ TclParseNumber( 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; - } + exponent = - exponent; } - - /* - * 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); + signum, significandWide, numSigDigs, + (numTrailZeros + exponent - numDigitsAfterDp)); } else { objPtr->internalRep.doubleValue = MakeHighPrecisionDouble( - signum, &significandBig, numSigDigs, exponent); + signum, &significandBig, numSigDigs, + (numTrailZeros + exponent - numDigitsAfterDp)); } break; @@ -1575,7 +1323,7 @@ TclParseNumber( break; #endif case INITIAL: - /* This case only to silence compiler warning. */ + /* This case only to silence compiler warning */ Tcl_Panic("TclParseNumber: state INITIAL can't happen here"); } } @@ -1586,16 +1334,18 @@ TclParseNumber( if (status != TCL_OK) { if (interp != NULL) { - Tcl_Obj *msg = Tcl_ObjPrintf("expected %s but got \"", - expected); + Tcl_Obj *msg; + TclNewLiteralStringObj(msg, "expected "); + Tcl_AppendToObj(msg, expected, -1); + Tcl_AppendToObj(msg, " but got \"", -1); Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, ""); Tcl_AppendToObj(msg, "\"", -1); if (state == BAD_OCTAL) { Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1); } Tcl_SetObjResult(interp, msg); - Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", (void *)NULL); + Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", NULL); } } @@ -1645,7 +1395,7 @@ AccumulateDecimalDigit( Tcl_WideUInt w; /* - * Try wide multiplication first. + * Try wide multiplication first */ if (!bignumFlag) { @@ -1658,20 +1408,17 @@ AccumulateDecimalDigit( *wideRepPtr = digit; return 0; } else if (numZeros >= maxpow10_wide - || w > (UWIDE_MAX-digit)/pow10_wide[numZeros+1]) { + || w > ((~(Tcl_WideUInt)0)-digit)/pow10_wide[numZeros+1]) { /* - * Wide multiplication will overflow. Expand the number to a - * bignum and fall through into the bignum case. + * 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; - } + TclBNInitBignumFromWideUInt(bignumRepPtr, w); } else { /* * Wide multiplication. */ - *wideRepPtr = w * pow10_wide[numZeros+1] + digit; return 0; } @@ -1686,37 +1433,32 @@ AccumulateDecimalDigit( * 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; + mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1], + bignumRepPtr); + mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr); } 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 + * seems implausible. We presume that DIGIT_BIT is at least 27. The * first multiplication, by up to 10**7, is done with a one-DIGIT - * multiply (this presumes that MP_DIGIT_BIT >= 24). + * multiply (this presumes that 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) { + mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr); + for (i=3; i<=7; ++i) { if (n & (1 << i)) { - err = mp_mul(bignumRepPtr, pow5+i, bignumRepPtr); + mp_mul(bignumRepPtr, pow5+i, bignumRepPtr); } } - while ((err == MP_OKAY) && (n >= 256)) { - err = mp_mul(bignumRepPtr, pow5+8, bignumRepPtr); + while (n >= 256) { + 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; - } + mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr); + mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr); } return 1; @@ -1746,35 +1488,32 @@ 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 */ + int exponent) /* Power of ten */ { - TCL_IEEE_DOUBLE_ROUNDING_DECL - - mp_int significandBig; /* Significand expressed as a bignum. */ + double retval; /* Value of the number */ + mp_int significandBig; /* Significand expressed as a bignum */ /* * With gcc on x86, the floating point rounding mode is double-extended. * This causes the result of double-precision calculations to be rounded * twice: once to the precision of double-extended and then again to the * precision of double. Double-rounding introduces gratuitous errors of 1 - * ulp, so we need to change rounding mode to 53-bits. We also make - * 'retval' volatile, so that it doesn't get promoted to a register. + * ulp, so we need to change rounding mode to 53-bits. */ - 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); - } +#if defined(__GNUC__) && defined(__i386) + fpu_control_t roundTo53Bits = 0x027f; + fpu_control_t oldRoundingMode; + _FPU_GETCW(oldRoundingMode); + _FPU_SETCW(roundTo53Bits); +#endif +#if defined(__sun) && defined(__i386) && !defined(__GNUC__) + ieee_flags("set","precision","double",NULL); +#endif /* - * Set the FP control word for 53 bits, WARNING: It must be reset - * before returning. + * Test for the easy cases. */ - TCL_IEEE_DOUBLE_ROUNDING; if (numSigDigs <= QUICK_MAX) { if (exponent >= 0) { @@ -1785,12 +1524,10 @@ MakeLowPrecisionDouble( * without special handling. */ - retval = (double) - ((Tcl_WideInt)significand * pow10vals[exponent]); + 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 @@ -1799,8 +1536,8 @@ MakeLowPrecisionDouble( * with only one roundoff. */ - volatile double factor = (double) - ((Tcl_WideInt)significand * pow10vals[diff]); + volatile double factor = + (double)(Tcl_WideInt)significand * pow10vals[diff]; retval = factor * pow10vals[exponent-diff]; goto returnValue; } @@ -1813,21 +1550,18 @@ MakeLowPrecisionDouble( * only one rounding. */ - retval = (double) - ((Tcl_WideInt)significand / pow10vals[-exponent]); + retval = (double)(Tcl_WideInt)significand / pow10vals[-exponent]; goto returnValue; } } } /* - * All the easy cases have failed. Promote the significand to bignum and + * All the easy cases have failed. Promote ths significand to bignum and * call MakeHighPrecisionDouble to do it the hard way. */ - if (mp_init_u64(&significandBig, significand) != MP_OKAY) { - return 0.0; - } + TclBNInitBignumFromWideUInt(&significandBig, significand); retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs, exponent); mp_clear(&significandBig); @@ -1845,7 +1579,12 @@ MakeLowPrecisionDouble( * On gcc on x86, restore the floating point mode word. */ - TCL_DEFAULT_DOUBLE_ROUNDING; +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW(oldRoundingMode); +#endif +#if defined(__sun) && defined(__i386) && !defined(__GNUC__) + ieee_flags("clear","precision",NULL,NULL); +#endif return retval; } @@ -1873,46 +1612,39 @@ 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 */ + int exponent) /* Power of 10 by which to multiply */ { - TCL_IEEE_DOUBLE_ROUNDING_DECL - - int machexp = 0; /* Machine exponent of a power of 10. */ + double retval; + int machexp; /* Machine exponent of a power of 10 */ /* * With gcc on x86, the floating point rounding mode is double-extended. * This causes the result of double-precision calculations to be rounded * twice: once to the precision of double-extended and then again to the * precision of double. Double-rounding introduces gratuitous errors of 1 - * ulp, so we need to change rounding mode to 53-bits. We also make - * 'retval' volatile to make sure that it doesn't get promoted to a - * register. + * ulp, so we need to change rounding mode to 53-bits. */ - 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); - } +#if defined(__GNUC__) && defined(__i386) + fpu_control_t roundTo53Bits = 0x027f; + fpu_control_t oldRoundingMode; + _FPU_GETCW(oldRoundingMode); + _FPU_SETCW(roundTo53Bits); +#endif +#if defined(__sun) && defined(__i386) && !defined(__GNUC__) + ieee_flags("set","precision","double",NULL); +#endif /* - * Set the 53-bit rounding mode. WARNING: It must be reset before - * returning. + * Quick checks for over/underflow. */ - 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) { + if (numSigDigs+exponent-1 > maxDigits) { retval = HUGE_VAL; goto returnValue; - } else if (exponent < 0 && numSigDigs+exponent < minDigits+1) { - retval = 0.0; + } + if (numSigDigs+exponent-1 < minDigits) { + retval = 0; goto returnValue; } @@ -1932,9 +1664,9 @@ MakeHighPrecisionDouble( goto returnValue; } retval = SafeLdExp(retval, machexp); - if (tiny == 0.0) { - tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits); - } + if (tiny == 0.0) { + tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits); + } if (retval < tiny) { retval = tiny; } @@ -1960,8 +1692,12 @@ MakeHighPrecisionDouble( * On gcc on x86, restore the floating point mode word. */ - TCL_DEFAULT_DOUBLE_ROUNDING; - +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW(oldRoundingMode); +#endif +#if defined(__sun) && defined(__i386) && !defined(__GNUC__) + ieee_flags("clear","precision",NULL,NULL); +#endif return retval; } @@ -1980,8 +1716,8 @@ MakeHighPrecisionDouble( #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. */ + int signum, /* Sign bit (1=negative, 0=nonnegative */ + Tcl_WideUInt tags) /* Tag bits to put in the NaN */ { union { Tcl_WideUInt iv; @@ -2019,41 +1755,37 @@ MakeNaN( static double RefineApproximation( - double approxResult, /* Approximate result of conversion. */ - mp_int *exactSignificand, /* Integer significand. */ - int exponent) /* Power of 10 to multiply by significand. */ + 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. */ + double significand; /* Sigificand of the binary number */ + int binExponent; /* Exponent of the binary number */ int msb; /* Most significant bit position of an - * intermediate result. */ + * intermediate result */ int nDigits; /* Number of mp_digit's in an intermediate - * result. */ + * result */ mp_int twoMv; /* Approx binary value expressed as an exact - * integer scaled by the multiplier 2M. */ + * 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. */ + * 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. */ + * 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, @@ -2063,22 +1795,13 @@ RefineApproximation( 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. + * Find a common denominator for the decimal and binary fractions. The + * common denominator will be 2**M2 + 5**M5. */ + significand = frexp(approxResult, &binExponent); i = mantBits - binExponent; if (i < 0) { M2 = 0; @@ -2089,153 +1812,102 @@ RefineApproximation( M5 = 0; } else { M5 = -exponent; - if (M5 - 1 > M2) { - M2 = M5 - 1; + 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. + * The floating point number is significand*2**binExponent. Compute the + * large integer significand*2**(binExponent+M2+1). The 2**-1 bit of the + * significand (the most significant) corresponds to the + * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold + * that quantity, then convert the significand to a large integer, scaled + * appropriately. Then multiply by the appropriate power of 5. */ msb = binExponent + M2; /* 1008 */ - nDigits = msb / MP_DIGIT_BIT + 1; - if (mp_init_size(&twoMv, nDigits) != MP_OKAY) { - return approxResult; - } - i = (msb % MP_DIGIT_BIT + 1); + nDigits = msb / DIGIT_BIT + 1; + mp_init_size(&twoMv, nDigits); + i = (msb % DIGIT_BIT + 1); twoMv.used = nDigits; significand *= SafeLdExp(1.0, i); while (--nDigits >= 0) { twoMv.dp[nDigits] = (mp_digit) significand; significand -= (mp_digit) significand; - significand = SafeLdExp(significand, MP_DIGIT_BIT); + significand = SafeLdExp(significand, 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; + if (M5 & (1 << i)) { + mp_mul(&twoMv, pow5+i, &twoMv); } } /* - * Compute twoMd as 2*M*d, where d is the exact value. - * This is done by multiplying by 5**(M5+exponent) and then multiplying - * by 2**(M5+exponent+1), which is, of course, a left shift. + * Collect the decimal significand as a high precision integer. The least + * significant bit corresponds to bit M2+exponent+1 so it will need to be + * shifted left by that many bits after being multiplied by + * 5**(M5+exponent). */ - 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); + mp_init_copy(&twoMd, exactSignificand); + for (i=0; i<=8; ++i) { + if ((M5+exponent) & (1 << i)) { + 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); - } + mp_mul_2d(&twoMd, M2+exponent+1, &twoMd); + mp_sub(&twoMd, &twoMv, &twoMd); /* * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction * term. Because 2M may well overflow a double, we need to scale the - * denominator by a factor of 2**binExponent-mantBits. Place that factor - * times 1/2 ULP into twoMd. + * denominator by a factor of 2**binExponent-mantBits */ scale = binExponent - mantBits - 1; - mp_set_u64(&twoMv, 1); - for (i = 0; (i <= 8) && (err == MP_OKAY); ++i) { + + mp_set(&twoMv, 1); + for (i=0; i<=8; ++i) { if (M5 & (1 << i)) { - err = mp_mul(&twoMv, pow5+i, &twoMv); + 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); + if (multiplier > 0) { + 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; + mp_div_2d(&twoMv, -multiplier, &twoMv, NULL); } - /* - * Will the eventual correction term be less than, equal to, or - * greater than 1/2 ULP? - */ - switch (mp_cmp_mag(&twoMd, &twoMv)) { case MP_LT: /* - * If the error is less than 1/2 ULP, there's no correction to make. + * If the result is less than unity, the error is less than 1/2 unit in + * the last place, so there's no correction to make. */ - mp_clear(&twoMd); - mp_clear(&twoMv); + 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. + * If the result is exactly unity, 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); + 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. */ @@ -2259,55 +1931,51 @@ RefineApproximation( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * MultPow5 -- * * Multiply a bignum by a power of 5. * * Side effects: - * Stores base*5**n in result. + * 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. */ +inline static void +MulPow5(mp_int* base, /* Number to multiply */ + unsigned n, /* Power of 5 to multiply by */ + mp_int* result) /* Place to store the result */ { - mp_int *p = base; + 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); + mp_mul_d(p, dpow5[r], result); p = result; } r = 0; - while ((err == MP_OKAY) && (n13 != 0)) { + while (n13 != 0) { if (n13 & 1) { - err = mp_mul(p, pow5_13+r, result); + mp_mul(p, pow5_13+r, result); p = result; } n13 >>= 1; ++r; } - if ((err == MP_OKAY) && (p != result)) { - err = mp_copy(p, result); + if (p != result) { + mp_copy(p, result); } - return err; } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * NormalizeRightward -- * - * Shifts a number rightward until it is odd (that is, until the least - * significant bit is nonzero. + * 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. @@ -2315,26 +1983,25 @@ MulPow5( * 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. */ +inline static int +NormalizeRightward(Tcl_WideUInt* wPtr) + /* INOUT: Number to shift */ { int rv = 0; Tcl_WideUInt w = *wPtr; - - if (!(w & (Tcl_WideUInt) 0xFFFFFFFF)) { + if (!(w & (Tcl_WideUInt) 0xffffffff)) { w >>= 32; rv += 32; } - if (!(w & (Tcl_WideUInt) 0xFFFF)) { + if (!(w & (Tcl_WideUInt) 0xffff)) { w >>= 16; rv += 16; } - if (!(w & (Tcl_WideUInt) 0xFF)) { + if (!(w & (Tcl_WideUInt) 0xff)) { w >>= 8; rv += 8; } - if (!(w & (Tcl_WideUInt) 0xF)) { + if (!(w & (Tcl_WideUInt) 0xf)) { w >>= 4; rv += 4; } if (!(w & 0x3)) { @@ -2346,43 +2013,42 @@ NormalizeRightward( *wPtr = w; return rv; } - + /* - *---------------------------------------------------------------------- + *-----------------------------------------------------------------------------0 * * RequiredPrecision -- * - * Determines the number of bits needed to hold an integer. + * Determines the number of bits needed to hold an intger. * * Results: - * Returns the position of the most significant bit (0 - 63). Returns 0 - * if the number is zero. + * 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. */ +RequiredPrecision(Tcl_WideUInt w) + /* Number to interrogate */ { int rv; unsigned long wi; - - if (w & ((Tcl_WideUInt) 0xFFFFFFFF << 32)) { + if (w & ((Tcl_WideUInt) 0xffffffff << 32)) { wi = (unsigned long) (w >> 32); rv = 32; } else { wi = (unsigned long) w; rv = 0; } - if (wi & 0xFFFF0000) { + if (wi & 0xffff0000) { wi >>= 16; rv += 16; } - if (wi & 0xFF00) { + if (wi & 0xff00) { wi >>= 8; rv += 8; } - if (wi & 0xF0) { + if (wi & 0xf0) { wi >>= 4; rv += 4; } - if (wi & 0xC) { + if (wi & 0xc) { wi >>= 2; rv += 2; } if (wi & 0x2) { @@ -2393,40 +2059,38 @@ RequiredPrecision( } 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'. + * 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. */ +inline static void +DoubleToExpAndSig(double dv, /* Number to convert */ + Tcl_WideUInt* significand, + /* OUTPUT: Significand of the number */ + int* expon, /* OUTPUT: Exponent to multiply the number by */ + int* bits) /* OUTPUT: Number of significant bits */ { - Double d; /* Number being converted. */ - Tcl_WideUInt z; /* Significand under construction. */ - int de; /* Exponent of the number. */ - int k; /* Bit count. */ + 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. - */ + /* Extract exponent and significand */ de = (d.w.word0 & EXP_MASK) >> EXP_SHIFT; z = d.q & SIG_MASK; @@ -2442,25 +2106,24 @@ DoubleToExpAndSig( } *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. + * 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. */ +inline static void +TakeAbsoluteValue(Double* d, /* Number to replace with absolute value */ + int* sign) /* Place to put the signum */ { if (d->w.word0 & SIGN_BIT) { *sign = 1; @@ -2469,42 +2132,41 @@ TakeAbsoluteValue( *sign = 0; } } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * FormatInfAndNaN -- * * Bailout for formatting infinities and Not-A-Number. * * Results: - * Returns one of the strings 'Infinity' and 'NaN'. The string returned - * must be freed by the caller using 'ckfree'. + * Returns one of the strings 'Infinity' and 'NaN'. * * Side effects: - * Stores 9999 in *decpt, and sets '*endPtr' to designate the terminating - * NUL byte of the string if 'endPtr' is not NULL. + * Stores 9999 in *decpt, and sets '*endPtr' to designate the + * terminating NUL byte of the string if 'endPtr' is not NULL. * - *---------------------------------------------------------------------- + * The string returned must be freed by the caller using 'ckfree'. + * + *----------------------------------------------------------------------------- */ -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 */ +inline static char* +FormatInfAndNaN(Double* d, /* Exceptional number to format */ + int* decpt, /* Decimal point to set to a bogus value */ + char** endPtr) /* Pointer to the end of the formatted data */ { - char *retval; - + char* retval; *decpt = 9999; if (!(d->w.word1) && !(d->w.word0 & HI_ORDER_SIG_MASK)) { - retval = (char *)ckalloc(9); + retval = ckalloc(9); strcpy(retval, "Infinity"); if (endPtr) { *endPtr = retval + 8; } } else { - retval = (char *)ckalloc(4); + retval = ckalloc(4); strcpy(retval, "NaN"); if (endPtr) { *endPtr = retval + 3; @@ -2512,9 +2174,9 @@ FormatInfAndNaN( } return retval; } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * FormatZero -- * @@ -2527,16 +2189,14 @@ FormatInfAndNaN( * 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 */ +inline static char* +FormatZero(int* decpt, /* Location of the decimal point */ + char** endPtr) /* Pointer to the end of the formatted data */ { - char *retval = (char *)ckalloc(2); - + char* retval = ckalloc(2); strcpy(retval, "0"); if (endPtr) { *endPtr = retval+1; @@ -2544,31 +2204,31 @@ FormatZero( *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. + * 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. + * 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. */ +inline static int +ApproximateLog10(Tcl_WideUInt bw, + /* Integer significand of the number */ + int be, /* Power of two to scale bw */ + int bbits) /* Number of bits of precision in bw */ { - int i; /* Log base 2 of the number. */ + int i; /* Log base 2 of the number */ int k; /* Floor(Log base 10 of the number) */ - double ds; /* Mantissa of the number. */ + double ds; /* Mantissa of the number */ Double d2; /* @@ -2582,16 +2242,17 @@ ApproximateLog10( 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; + + LOG10_3HALVES_PLUS_FUDGE + + LOG10_2 * i; k = (int) ds; if (k > ds) { --k; } return k; } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * BetterLog10 -- * @@ -2599,27 +2260,24 @@ ApproximateLog10( * 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. + * 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). + * 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. */ +inline static int +BetterLog10(double d, /* Original number to format */ + int k, /* Characteristic(Log base 10) of the number */ + int* k_check) /* Flag == 1 if k is inexact */ { /* - * Performance hack. If k is in the range 0..TEN_PMAX, then we can use a - * powers-of-ten table to check it. + * 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--; @@ -2630,41 +2288,40 @@ BetterLog10( } 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. + * 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. + * 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. */ +inline static void +ComputeScale(int be, /* Exponent part of number: d = bw * 2**be */ + int k, /* Characteristic of log10(number) */ + int* b2, /* OUTPUT: Power of 2 in the numerator */ + int* b5, /* OUTPUT: Power of 5 in the numerator */ + int* s2, /* OUTPUT: Power of 2 in the denominator */ + int* s5) /* OUTPUT: Power of 5 in the denominator */ { + /* - * Scale numerator and denominator powers of 2 so that the input binary - * number is the ratio of integers. + * 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; @@ -2674,10 +2331,9 @@ ComputeScale( } /* - * Scale numerator and denominator so that the output decimal number is - * the ratio of integers. + * Scale numerator and denominator so that the output decimal number + * is the ratio of integers */ - if (k >= 0) { *b5 = 0; *s5 = k; @@ -2690,44 +2346,55 @@ ComputeScale( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * SetPrecisionLimits -- * - * Determines how many digits of significance should be computed (and, - * hence, how much memory need be allocated) for formatting a floating - * point number. + * 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 '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. + * 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. */ +inline static void +SetPrecisionLimits(int convType, + /* Type of conversion: + * TCL_DD_SHORTEST + * TCL_DD_STEELE0 + * TCL_DD_E_FMT + * TCL_DD_F_FMT */ + int k, /* Floor(log10(number to convert)) */ + int* ndigitsPtr, + /* IN/OUT: Number of digits requested + * (Will be adjusted if needed) */ + int* iPtr, /* OUT: Maximum number of digits + * to return */ + int *iLimPtr,/* OUT: Number of digits of significance + * if the bignum method is used.*/ + int *iLim1Ptr) + /* OUT: Number of digits of significance + * if the quick method is used. */ { - switch (flags & TCL_DD_CONVERSION_TYPE_MASK) { + switch(convType) { + case TCL_DD_SHORTEST0: + case TCL_DD_STEELE0: + *iLimPtr = *iLim1Ptr = -1; + *iPtr = 18; + *ndigitsPtr = 0; + break; case TCL_DD_E_FORMAT: if (*ndigitsPtr <= 0) { *ndigitsPtr = 1; @@ -2743,37 +2410,37 @@ SetPrecisionLimits( } break; default: - *iLimPtr = *iLim1Ptr = -1; - *iPtr = 18; - *ndigitsPtr = 0; - break; + *iPtr = -1; + *iLimPtr = -1; + *iLim1Ptr = -1; + Tcl_Panic("impossible conversion type in TclDoubleDigits"); } } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * BumpUp -- * - * Increases a string of digits ending in a series of nines to designate - * the next higher number. xxxxb9999... -> xxxx(b+1)0000... + * 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. + * 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. */ + +inline static char* +BumpUp(char* s, /* Cursor pointing one past the end of the + * string */ + char* retval, /* Start of the string of digits */ + int* kPtr) /* Position of the decimal point */ { while (*--s == '9') { if (s == retval) { @@ -2788,28 +2455,27 @@ BumpUp( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * AdjustRange -- * - * Rescales a 'double' in preparation for formatting it using the 'quick' - * double-to-string method. + * 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. + * 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)) */ +inline static int +AdjustRange(double* dPtr, /* INOUT: Number to adjust */ + int k) /* IN: floor(log10(d)) */ { int ieps; /* Number of roundoff errors that have - * accumulated. */ - double d = *dPtr; /* Number to adjust. */ + * accumulated */ + double d = *dPtr; /* Number to adjust */ double ds; int i, j, j1; @@ -2819,8 +2485,7 @@ AdjustRange( /* * The number must be reduced to bring it into range. */ - - ds = tens[k & 0xF]; + ds = tens[k & 0xf]; j = k >> 4; if (j & BLETCH) { j &= (BLETCH-1); @@ -2838,10 +2503,9 @@ AdjustRange( d /= ds; } else if ((j1 = -k) != 0) { /* - * The number must be increased to bring it into range. + * The number must be increased to bring it into range */ - - d *= tens[j1 & 0xF]; + d *= tens[j1 & 0xf]; i = 0; for (j = j1>>4; j; j>>=1) { if (j & 1) { @@ -2857,52 +2521,52 @@ AdjustRange( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * 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. + * 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. + * 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. */ +inline static char* +ShorteningQuickFormat(double d, /* Number to convert */ + int k, /* floor(log10(d)) */ + int ilim, /* Number of significant digits to return */ + double eps, + /* Estimated roundoff error */ + char* retval, + /* Buffer to receive the digit string */ + int* kPtr) + /* Pointer to stash the position of + * the decimal point */ { - char *s = retval; /* Cursor in the return value. */ - int digit; /* Current digit. */ + 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. - */ + /* 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. + * Truncate the conversion if the string of digits is within + * 1/2 ulp of the actual value. */ if (d < eps) { @@ -2916,7 +2580,7 @@ ShorteningQuickFormat( /* * Bail out if the conversion fails to converge to a sufficiently - * precise value. + * precise value */ if (++i >= ilim) { @@ -2933,44 +2597,40 @@ ShorteningQuickFormat( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * StrictQuickFormat -- * - * Convert a double precision number of a string of a precise number of - * digits, using the 'quick' double precision method. + * 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. + * 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. */ +inline static char* +StrictQuickFormat(double d, /* Number to convert */ + int k, /* floor(log10(d)) */ + int ilim, /* Number of significant digits to return */ + double eps, /* Estimated roundoff error */ + char* retval, /* Start of the digit string */ + int* kPtr) /* Pointer to stash the position of + * the decimal point */ { - char *s = retval; /* Cursor in the return value. */ - int digit; /* Current digit of the answer. */ + 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. - */ - + /* Extract a digit */ digit = (int) d; d -= digit; if (d == 0.0) { @@ -2979,10 +2639,9 @@ StrictQuickFormat( *s++ = '0' + digit; /* - * When the given digit count is reached, handle trailing strings of 0 - * and 9. + * When the given digit count is reached, handle trailing strings + * of 0 and 9. */ - if (i == ilim) { if (d > 0.5 + eps) { *kPtr = k; @@ -2999,17 +2658,14 @@ StrictQuickFormat( } } - /* - * Advance to the next digit. - */ - + /* Advance to the next digit */ ++i; d *= 10.0; } } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * QuickConversion -- * @@ -3018,72 +2674,67 @@ StrictQuickFormat( * therefore be used for the intermediate results. * * Results: - * Returns the converted string, or NULL if the bignum method must be - * used. + * 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. */ +inline static char* +QuickConversion(double e, /* Number to format */ + int k, /* floor(log10(d)), approximately */ + int k_check, /* 0 if k is exact, 1 if it may be too high */ + int flags, /* Flags passed to dtoa: + * TCL_DD_SHORTEN_FLAG */ + int len, /* Length of the return value */ + int ilim, /* Number of digits to store */ + int ilim1, /* Number of digits to store if we + * musguessed k */ + int* decpt, /* OUTPUT: Location of the decimal point */ + char** endPtr) /* OUTPUT: Pointer to the terminal null byte */ { int ieps; /* Number of 1-ulp roundoff errors that have - * accumulated in the calculation. */ - Double eps; /* Estimated roundoff error. */ - char *retval; /* Returned string. */ - char *end; /* Pointer to the terminal null byte in the - * returned string. */ + * 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). + * 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 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; + d *= 10.0; ++ieps; } /* - * Compute estimated roundoff error. + * 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. + * Handle the peculiar case where the result has no significant + * digits. */ - - retval = (char *)ckalloc(len + 1); + retval = ckalloc(len + 1); if (ilim == 0) { - d = d - 5.; + d -= 5.; if (d > eps.d) { *retval = '1'; *decpt = k; @@ -3097,11 +2748,9 @@ QuickConversion( } } - /* - * Format the digit string. - */ + /* Format the digit string */ - if (flags & TCL_DD_SHORTEST) { + if (flags & TCL_DD_SHORTEN_FLAG) { end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt); } else { end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt); @@ -3118,97 +2767,106 @@ QuickConversion( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * 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. + * 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. */ +inline static void +CastOutPowersOf2(int* b2, /* Power of 2 to multiply the significand */ + int* m2, /* Power of 2 to multiply 1/2 ulp */ + int* s2) /* Power of 2 to multiply the common + * denominator */ { int i; - if (*m2 > 0 && *s2 > 0) { /* Find the smallest power of 2 in the - * numerator. */ - if (*m2 < *s2) { /* Find the lowest common denominator. */ + * numerator */ + if (*m2 < *s2) { /* Find the lowest common denominatorr */ i = *m2; } else { i = *s2; } - *b2 -= i; /* Reduce to lowest terms. */ + *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. + * 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 + * 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'. + * 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. */ +inline static char* +ShorteningInt64Conversion(Double* dPtr, + /* Original number to convert */ + int convType, + /* Type of conversion (shortest, Steele, + E format, F format) */ + Tcl_WideUInt bw, + /* Integer significand */ + int b2, int b5, + /* Scale factor for the significand + * in the numerator */ + int m2plus, int m2minus, int m5, + /* Scale factors for 1/2 ulp in + * the numerator (will be different if + * bw == 1 */ + int s2, int s5, + /* Scale factors for the denominator */ + int k, + /* Number of output digits before the decimal + * point */ + int len, + /* Number of digits to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Position of the terminal '\0' + * at the end of the returned string */ { - char *retval = (char *)ckalloc(len + 1); - /* Output buffer. */ + + char* retval = ckalloc(len + 1); + /* Output buffer */ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2; - /* Numerator of the fraction being - * converted. */ + /* 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. */ + * 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. - */ + /* Adjust if the logarithm was guessed wrong */ if (b < S) { b = 10 * b; @@ -3217,16 +2875,12 @@ ShorteningInt64Conversion( --k; } - /* - * Compute roundoff ranges. - */ + /* Compute roundoff ranges */ mplus = wuipow5[m5] << m2plus; mminus = wuipow5[m5] << m2minus; - /* - * Loop through the digits. - */ + /* Loop through the digits */ i = 1; for (;;) { @@ -3240,15 +2894,17 @@ ShorteningInt64Conversion( * Does the current digit put us on the low side of the exact value * but within within roundoff of being exact? */ - - if (b < mplus || (b == mplus + if (b < mplus + || (b == mplus + && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { /* - * Make sure we shouldn't be rounding *up* instead, in case the - * next number above is closer. + * 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)) { + if (2 * b > S + || (2 * b == S + && (digit & 1) != 0)) { ++digit; if (digit == 10) { *s++ = '9'; @@ -3257,9 +2913,7 @@ ShorteningInt64Conversion( } } - /* - * Stash the current digit. - */ + /* Stash the current digit */ *s++ = '0' + digit; break; @@ -3269,8 +2923,9 @@ ShorteningInt64Conversion( * Does one plus the current digit put us within roundoff of the * number? */ - - if (b > S - mminus || (b == S - mminus + if (b > S - mminus + || (b == S - mminus + && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { if (digit == 9) { *s++ = '9'; @@ -3285,18 +2940,16 @@ ShorteningInt64Conversion( /* * Have we converted all the requested digits? */ - *s++ = '0' + digit; if (i == ilim) { - if (2*b > S || (2*b == S && (digit & 1) != 0)) { + if (2*b > S + || (2*b == S && (digit & 1) != 0)) { s = BumpUp(s, retval, &k); } break; } - /* - * Advance to the next digit. - */ + /* Advance to the next digit */ b = 10 * b; mplus = 10 * mplus; @@ -3308,7 +2961,6 @@ ShorteningInt64Conversion( * Endgame - store the location of the decimal point and the end of the * string. */ - *s = '\0'; *decpt = k; if (endPtr) { @@ -3318,58 +2970,69 @@ ShorteningInt64Conversion( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * 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. + * 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 + * 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'. + * 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. */ +inline static char* +StrictInt64Conversion(Double* dPtr, + /* Original number to convert */ + int convType, + /* Type of conversion (shortest, Steele, + E format, F format) */ + Tcl_WideUInt bw, + /* Integer significand */ + int b2, int b5, + /* Scale factor for the significand + * in the numerator */ + int s2, int s5, + /* Scale factors for the denominator */ + int k, + /* Number of output digits before the decimal + * point */ + int len, + /* Number of digits to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Position of the terminal '\0' + * at the end of the returned string */ { - char *retval = (char *)ckalloc(len + 1); - /* Output buffer. */ + + char* retval = ckalloc(len + 1); + /* Output buffer */ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2; - /* Numerator of the fraction being - * converted. */ + /* 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. */ + * 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. - */ + /* Adjust if the logarithm was guessed wrong */ if (b < S) { b = 10 * b; @@ -3377,9 +3040,7 @@ StrictInt64Conversion( --k; } - /* - * Loop through the digits. - */ + /* Loop through the digits */ i = 1; for (;;) { @@ -3392,10 +3053,10 @@ StrictInt64Conversion( /* * Have we converted all the requested digits? */ - *s++ = '0' + digit; if (i == ilim) { - if (2*b > S || (2*b == S && (digit & 1) != 0)) { + if (2*b > S + || (2*b == S && (digit & 1) != 0)) { s = BumpUp(s, retval, &k); } else { while (*--s == '0') { @@ -3406,9 +3067,7 @@ StrictInt64Conversion( break; } - /* - * Advance to the next digit. - */ + /* Advance to the next digit */ b = 10 * b; ++i; @@ -3418,7 +3077,6 @@ StrictInt64Conversion( * Endgame - store the location of the decimal point and the end of the * string. */ - *s = '\0'; *decpt = k; if (endPtr) { @@ -3428,30 +3086,30 @@ StrictInt64Conversion( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * 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. + * Test whether bankers' rounding should round a digit up. Assumption + * is made that the denominator of the fraction being tested is + * a power of 2**DIGIT_BIT. * * Results: - * Returns 1 iff the fraction is more than 1/2, or if the fraction is - * exactly 1/2 and the digit is odd. + * 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. */ +inline static int +ShouldBankerRoundUpPowD(mp_int* b, + /* Numerator of the fraction */ + int sd, /* Denominator is 2**(sd*DIGIT_BIT) */ + int isodd) + /* 1 if the digit is odd, 0 if even */ { int i; - static const mp_digit topbit = ((mp_digit)1) << (MP_DIGIT_BIT - 1); - + static const mp_digit topbit = (1<<(DIGIT_BIT-1)); if (b->used < sd || (b->dp[sd-1] & topbit) == 0) { return 0; } @@ -3467,37 +3125,45 @@ ShouldBankerRoundUpPowD( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * ShouldBankerRoundUpToNextPowD -- * - * Tests whether bankers' rounding will round down in the "denominator is - * a power of 2**MP_DIGIT" case. + * 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. */ +inline static int +ShouldBankerRoundUpToNextPowD(mp_int* b, + /* Numerator of the fraction */ + mp_int* m, + /* Numerator of the rounding tolerance */ + int sd, + /* Common denominator is 2**(sd*DIGIT_BIT) */ + int convType, + /* Conversion type: STEELE defeats + * round-to-even (Not sure why one wants to + * do this; I copied it from Gay) FIXME */ + int isodd, + /* 1 if the integer significand is odd */ + mp_int* temp) + /* Work area for the calculation */ { int i; /* - * Compare B and S-m - which is the same as comparing B+m and S - which we - * do by computing b+m and doing a bitwhack compare against - * 2**(MP_DIGIT_BIT*sd) + * Compare B and S-m -- which is the same as comparing B+m and S -- + * which we do by computing b+m and doing a bitwhack compare against + * 2**(DIGIT_BIT*sd) */ - - if ((mp_add(b, m, temp) != MP_OKAY) || (temp->used <= sd)) { /* Too few digits to be > s */ + mp_add(b, m, temp); + if (temp->used <= sd) { /* too few digits to be > S */ return 0; } if (temp->used > sd+1 || temp->dp[sd] > 1) { @@ -3505,91 +3171,98 @@ ShouldBankerRoundUpToNextPowD( return 1; } for (i = sd-1; i >= 0; --i) { - /* Check for ==s */ + /* check for ==s */ if (temp->dp[i] != 0) { /* > s */ return 1; } } + if (convType == TCL_DD_STEELE0) { + /* biased rounding */ + return 0; + } return isodd; } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * ShorteningBignumConversionPowD -- * - * Converts a double-precision number to the shortest string of digits - * that reconverts exactly to the given number, or to 'ilim' digits if - * that will yield a shorter result. The denominator in David Gay's - * conversion algorithm is known to be a power of 2**MP_DIGIT_BIT, and hence - * the division in the main loop may be replaced by a digit shift and - * mask. + * Converts a double-precision number to the shortest string of + * digits that reconverts exactly to the given number, or to + * 'ilim' digits if that will yield a shorter result. The denominator + * in David Gay's conversion algorithm is known to be a power of + * 2**DIGIT_BIT, and hence the division in the main loop may be replaced + * by a digit shift and mask. * * Results: - * Returns the string of significant decimal digits, in newly allocated - * memory + * 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'. + * 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. */ +inline static char* +ShorteningBignumConversionPowD(Double* dPtr, + /* Original number to convert */ + int convType, + /* Type of conversion (shortest, Steele, + E format, F format) */ + Tcl_WideUInt bw, + /* Integer significand */ + int b2, int b5, + /* Scale factor for the significand + * in the numerator */ + int m2plus, int m2minus, int m5, + /* Scale factors for 1/2 ulp in + * the numerator (will be different if + * bw == 1 */ + int sd, + /* Scale factor for the denominator */ + int k, + /* Number of output digits before the decimal + * point */ + int len, + /* Number of digits to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Position of the terminal '\0' + * at the end of the returned string */ { - char *retval = (char *)ckalloc(len + 1); - /* Output buffer. */ - mp_int b; /* Numerator of the fraction being - * converted. */ - mp_int mplus, mminus; /* Bounds for roundoff. */ - mp_digit digit; /* Current output digit. */ - char *s = retval; /* Cursor in the output buffer. */ - int i; /* Index in the output buffer. */ + + char* retval = ckalloc(len + 1); + /* Output buffer */ + mp_int b; /* Numerator of the fraction being converted */ + mp_int mplus, mminus; /* Bounds for roundoff */ + mp_digit digit; /* Current output digit */ + char* s = retval; /* Cursor in the output buffer */ + int i; /* Index in the output buffer */ mp_int temp; int r1; - 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); - } + TclBNInitBignumFromWideUInt(&b, bw); + mp_init_set_int(&mminus, 1); + MulPow5(&b, b5, &b); + mp_mul_2d(&b, b2, &b); - /* - * Adjust if the logarithm was guessed wrong. - */ + /* Adjust if the logarithm was guessed wrong */ - if ((err == MP_OKAY) && (b.used <= sd)) { - err = mp_mul_d(&b, 10, &b); + if (b.used <= sd) { + mp_mul_d(&b, 10, &b); ++m2plus; ++m2minus; ++m5; ilim = ilim1; --k; @@ -3600,26 +3273,16 @@ ShorteningBignumConversionPowD( * 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); + mp_mul_2d(&mminus, m2minus, &mminus); + MulPow5(&mminus, m5, &mminus); + if (m2plus > m2minus) { + mp_init_copy(&mplus, &mminus); + mp_mul_2d(&mplus, m2plus-m2minus, &mplus); } + mp_init(&temp); - /* - * Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT) - * by mp_digit extraction. - */ + /* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT) + * by mp_digit extraction */ i = 0; for (;;) { @@ -3639,13 +3302,14 @@ ShorteningBignumConversionPowD( */ r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus); - if (r1 == MP_LT || (r1 == MP_EQ + if (r1 == MP_LT + || (r1 == MP_EQ + && convType != TCL_DD_STEELE0 && (dPtr->w.word1 & 1) == 0)) { /* - * Make sure we shouldn't be rounding *up* instead, in case the - * next number above is closer. + * 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) { @@ -3655,9 +3319,7 @@ ShorteningBignumConversionPowD( } } - /* - * Stash the last digit. - */ + /* Stash the last digit */ *s++ = '0' + digit; break; @@ -3669,7 +3331,8 @@ ShorteningBignumConversionPowD( */ if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd, - dPtr->w.word1 & 1, &temp)) { + convType, dPtr->w.word1 & 1, + &temp)) { if (digit == 9) { *s++ = '9'; s = BumpUp(s, retval, &k); @@ -3683,7 +3346,6 @@ ShorteningBignumConversionPowD( /* * Have we converted all the requested digits? */ - *s++ = '0' + digit; if (i == ilim) { if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { @@ -3692,18 +3354,12 @@ ShorteningBignumConversionPowD( break; } - /* - * Advance to the next digit. - */ + /* 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); + mp_mul_d(&b, 10, &b); + mp_mul_d(&mminus, 10, &mminus); + if (m2plus > m2minus) { + mp_mul_2d(&mminus, m2plus-m2minus, &mplus); } ++i; } @@ -3712,94 +3368,101 @@ ShorteningBignumConversionPowD( * 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); + mp_clear_multi(&b, &mminus, &temp, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } - return (err == MP_OKAY) ? retval : NULL; + return retval; } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * StrictBignumConversionPowD -- * - * Converts a double-precision number to a fixed-lengt string of 'ilim' - * digits (or 'ilim1' if log10(d) has been overestimated). The - * denominator in David Gay's conversion algorithm is known to be a power - * of 2**MP_DIGIT_BIT, and hence the division in the main loop may be - * replaced by a digit shift and mask. + * Converts a double-precision number to a fixed-lengt string of + * 'ilim' digits (or 'ilim1' if log10(d) has been overestimated.) + * The denominator in David Gay's conversion algorithm is known to + * be a power of 2**DIGIT_BIT, and hence the division in the main + * loop may be replaced by a digit shift and mask. * * Results: - * Returns the string of significant decimal digits, in newly allocated - * memory. + * 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'. + * 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. */ +inline static char* +StrictBignumConversionPowD(Double* dPtr, + /* Original number to convert */ + int convType, + /* Type of conversion (shortest, Steele, + E format, F format) */ + Tcl_WideUInt bw, + /* Integer significand */ + int b2, int b5, + /* Scale factor for the significand + * in the numerator */ + int sd, + /* Scale factor for the denominator */ + int k, + /* Number of output digits before the decimal + * point */ + int len, + /* Number of digits to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Position of the terminal '\0' + * at the end of the returned string */ { - char *retval = (char *)ckalloc(len + 1); - /* Output buffer. */ - mp_int b; /* Numerator of the fraction being - * converted. */ - mp_digit digit; /* Current output digit. */ - char *s = retval; /* Cursor in the output buffer. */ - int i; /* Index in the output buffer. */ - mp_err err; + + char* retval = ckalloc(len + 1); + /* Output buffer */ + mp_int b; /* Numerator of the fraction being converted */ + mp_digit digit; /* Current output digit */ + char* s = retval; /* Cursor in the output buffer */ + int i; /* Index in the output buffer */ + mp_int temp; /* * b = bw * 2**b2 * 5**b5 */ - 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); - } + TclBNInitBignumFromWideUInt(&b, bw); + MulPow5(&b, b5, &b); + mp_mul_2d(&b, b2, &b); - /* - * Adjust if the logarithm was guessed wrong. - */ + /* Adjust if the logarithm was guessed wrong */ - if ((err == MP_OKAY) && (b.used <= sd)) { - err = mp_mul_d(&b, 10, &b); + if (b.used <= sd) { + mp_mul_d(&b, 10, &b); ilim = ilim1; --k; } + mp_init(&temp); /* - * Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT) - * by mp_digit extraction. + * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT) + * by mp_digit extraction */ i = 1; - while (err == MP_OKAY) { + for (;;) { if (b.used <= sd) { digit = 0; } else { @@ -3807,31 +3470,28 @@ StrictBignumConversionPowD( if (b.used > sd+1 || digit >= 10) { Tcl_Panic("wrong digit!"); } - --b.used; - mp_clamp(&b); + --b.used; mp_clamp(&b); } /* * Have we converted all the requested digits? */ - *s++ = '0' + digit; if (i == ilim) { if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { s = BumpUp(s, retval, &k); + } else { + while (*--s == '0') { + /* do nothing */ + } + ++s; } - while (*--s == '0') { - /* do nothing */ - } - ++s; break; } - /* - * Advance to the next digit. - */ + /* Advance to the next digit */ - err = mp_mul_d(&b, 10, &b); + mp_mul_d(&b, 10, &b); ++i; } @@ -3839,8 +3499,7 @@ StrictBignumConversionPowD( * Endgame - store the location of the decimal point and the end of the * string. */ - - mp_clear(&b); + mp_clear_multi(&b, &temp, NULL); *s = '\0'; *decpt = k; if (endPtr) { @@ -3850,7 +3509,7 @@ StrictBignumConversionPowD( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * ShouldBankerRoundUp -- * @@ -3860,74 +3519,82 @@ StrictBignumConversionPowD( * 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. */ +inline static int +ShouldBankerRoundUp(mp_int* twor, + /* 2x the remainder from thd division that + * produced the last digit */ + mp_int* S, /* Denominator */ + int isodd) /* Flag == 1 if the last digit is odd */ { int r = mp_cmp_mag(twor, S); - switch (r) { + case MP_LT: + return 0; case MP_EQ: return isodd; case MP_GT: return 1; - default: - return 0; } + Tcl_Panic("in ShouldBankerRoundUp, trichotomy fails!"); + return 0; } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * ShouldBankerRoundUpToNext -- * - * Tests whether the remainder is great enough to force rounding to the - * next higher digit. + * 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 +inline static int +ShouldBankerRoundUpToNext(mp_int* b, + /* Remainder from the division that produced * the last digit. */ - mp_int *m, /* Numerator of the rounding tolerance. */ - mp_int *S, /* Denominator. */ - int isodd) /* 1 if the integer significand is odd. */ + mp_int* m, + /* Numerator of the rounding tolerance */ + mp_int* S, + /* Denominator */ + int convType, + /* Conversion type: STEELE0 defeats + * round-to-even. (Not sure why one would + * want this; I coped it from Gay. FIXME */ + int isodd, + /* 1 if the integer significand is odd */ + mp_int* temp) + /* Work area needed for the calculation */ { int r; - 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); + /* Compare b and S-m: this is the same as comparing B+m and S. */ + mp_add(b, m, temp); + r = mp_cmp_mag(temp, S); switch(r) { + case MP_LT: + return 0; case MP_EQ: - return isodd; + if (convType == TCL_DD_STEELE0) { + return 0; + } else { + return isodd; + } case MP_GT: return 1; - default: - return 0; } + Tcl_Panic("in ShouldBankerRoundUpToNext, trichotomy fails!"); + return 0; } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * ShorteningBignumConversion -- * @@ -3938,97 +3605,90 @@ ShouldBankerRoundUpToNext( * 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. + * 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 */ +inline static char* +ShorteningBignumConversion(Double* dPtr, + /* Original number being converted */ + int convType, + /* Conversion type */ + Tcl_WideUInt bw, + /* Integer significand and exponent */ + int b2, + /* Scale factor for the significand */ + int m2plus, int m2minus, + /* Scale factors for 1/2 ulp in numerator */ + int s2, int s5, + /* Scale factors for denominator */ + int k, + /* Guessed position of the decimal point */ + int len, + /* Size of the digit buffer to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Pointer to the end of the number */ { - char *retval = (char *)ckalloc(len+1); - /* Buffer of digits to return. */ - char *s = retval; /* Cursor in the return value. */ - mp_int b; /* Numerator of the result. */ - mp_int mminus; /* 1/2 ulp below the result. */ - mp_int mplus; /* 1/2 ulp above the result. */ - mp_int S; /* Denominator of the result. */ - mp_int dig; /* Current digit of the result. */ - int digit; /* Current digit of the result. */ - int minit = 1; /* Fudge factor for when we misguess k. */ + char* retval = ckalloc(len+1); + /* Buffer of digits to return */ + char* s = retval; /* Cursor in the return value */ + mp_int b; /* Numerator of the result */ + mp_int mminus; /* 1/2 ulp below the result */ + mp_int mplus; /* 1/2 ulp above the result */ + mp_int S; /* Denominator of the result */ + mp_int dig; /* Current digit of the result */ + int digit; /* Current digit of the result */ + mp_int temp; /* Work area */ + int minit = 1; /* Fudge factor for when we misguess k */ int i; int r1; - 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); - } + TclBNInitBignumFromWideUInt(&b, bw); + mp_mul_2d(&b, b2, &b); + mp_init_set_int(&S, 1); + MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S); /* - * Handle the case where we guess the position of the decimal point wrong. + * 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); + if (mp_cmp_mag(&b, &S) == MP_LT) { + mp_mul_d(&b, 10, &b); minit = 10; ilim =ilim1; --k; } - /* - * mminus = 2**m2minus * 5**m5 - */ + /* 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); - } + mp_init_set_int(&mminus, minit); + mp_mul_2d(&mminus, m2minus, &mminus); + if (m2plus > m2minus) { + mp_init_copy(&mplus, &mminus); + mp_mul_2d(&mplus, m2plus-m2minus, &mplus); } + mp_init(&temp); - /* - * Loop through the digits. - */ + /* Loop through the digits */ - if (err == MP_OKAY) { - err = mp_init(&dig); - } + mp_init(&dig); i = 1; - while (err == MP_OKAY) { - err = mp_div(&b, &S, &dig, &b); + for (;;) { + mp_div(&b, &S, &dig, &b); if (dig.used > 1 || dig.dp[0] >= 10) { Tcl_Panic("wrong digit!"); } @@ -4040,8 +3700,11 @@ ShorteningBignumConversion( */ 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 (r1 == MP_LT + || (r1 == MP_EQ + && convType != TCL_DD_STEELE0 + && (dPtr->w.word1 & 1) == 0)) { + mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { ++digit; if (digit == 10) { @@ -4055,12 +3718,12 @@ ShorteningBignumConversion( } /* - * Does the current digit leave us with a remainder large enough to - * commit to rounding up to the next higher digit? + * 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)) { + if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType, + dPtr->w.word1 & 1, &temp)) { ++digit; if (digit == 10) { *s++ = '9'; @@ -4071,47 +3734,36 @@ ShorteningBignumConversion( break; } - /* - * Have we converted all the requested digits? - */ + /* Have we converted all the requested digits? */ *s++ = '0' + digit; - if ((err == MP_OKAY) && (i == ilim)) { - err = mp_mul_2d(&b, 1, &b); + if (i == ilim) { + mp_mul_2d(&b, 1, &b); if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); } break; } - /* - * Advance to the next digit. - */ + /* Advance to the next digit */ - if ((err == MP_OKAY) && (s5 > 0)) { - /* - * Can possibly shorten the denominator. - */ + if (s5 > 0) { - 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); + /* Can possibly shorten the denominator */ + mp_mul_2d(&b, 1, &b); + mp_mul_2d(&mminus, 1, &mminus); + if (m2plus > m2minus) { + mp_mul_2d(&mplus, 1, &mplus); } + mp_div_d(&S, 5, &S, NULL); --s5; - /* - * IDEA: It might possibly be a win to fall back to int64_t - * arithmetic here if S < 2**64/10. But it's a win only for - * a fairly narrow range of magnitudes so perhaps not worth - * bothering. We already know that we shorten the - * denominator by at least 1 mp_digit, perhaps 2, as we do - * the conversion for 17 digits of significance. + * IDEA: It might possibly be a win to fall back to + * int64 arithmetic here if S < 2**64/10. But it's + * a win only for a fairly narrow range of magnitudes + * so perhaps not worth bothering. We already know that + * we shorten the denominator by at least 1 mp_digit, perhaps + * 2. as we do the conversion for 17 digits of significance. * Possible savings: * 10**26 1 trip through loop before fallback possible * 10**27 1 trip @@ -4130,310 +3782,311 @@ ShorteningBignumConversion( * 10**40 14 trips * 10**41 15 trips * 10**42 16 trips - * thereafter no gain. + * 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); + } else { + mp_mul_d(&b, 10, &b); + mp_mul_d(&mminus, 10, &mminus); + if (m2plus > m2minus) { + mp_mul_2d(&mplus, 10, &mplus); } } ++i; } + /* * Endgame - store the location of the decimal point and the end of the * string. */ - if (m2plus > m2minus) { mp_clear(&mplus); } - mp_clear_multi(&b, &mminus, &dig, &S, (void *)NULL); + mp_clear_multi(&b, &mminus, &temp, &dig, &S, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; + } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * StrictBignumConversion -- * - * Convert a floating point number to a fixed-length digit string using - * the multiprecision method. + * 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. + * 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 */ +inline static char* +StrictBignumConversion(Double* dPtr, + /* Original number being converted */ + int convType, + /* Conversion type */ + Tcl_WideUInt bw, + /* Integer significand and exponent */ + int b2, /* Scale factor for the significand */ + int s2, int s5, + /* Scale factors for denominator */ + int k, /* Guessed position of the decimal point */ + int len, /* Size of the digit buffer to allocate */ + int ilim, + /* Number of digits to convert if b >= s */ + int ilim1, + /* Number of digits to convert if b < s */ + int* decpt, + /* OUTPUT: Position of the decimal point */ + char** endPtr) + /* OUTPUT: Pointer to the end of the number */ { - char *retval = (char *)ckalloc(len+1); - /* Buffer of digits to return. */ - char *s = retval; /* Cursor in the return value. */ - mp_int b; /* Numerator of the result. */ - mp_int S; /* Denominator of the result. */ - mp_int dig; /* Current digit of the result. */ - int digit; /* Current digit of the result. */ - int g; /* Size of the current digit ground. */ + char* retval = ckalloc(len+1); + /* Buffer of digits to return */ + char* s = retval; /* Cursor in the return value */ + mp_int b; /* Numerator of the result */ + mp_int S; /* Denominator of the result */ + mp_int dig; /* Current digit of the result */ + int digit; /* Current digit of the result */ + mp_int temp; /* Work area */ + int g; /* Size of the current digit groun */ int i, j; - 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); - } - } + mp_init_multi(&temp, &dig, NULL); + TclBNInitBignumFromWideUInt(&b, bw); + mp_mul_2d(&b, b2, &b); + mp_init_set_int(&S, 1); + MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S); /* - * Handle the case where we guess the position of the decimal point wrong. + * 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)) { + if (mp_cmp_mag(&b, &S) == MP_LT) { + mp_mul_d(&b, 10, &b); ilim =ilim1; --k; } - /* - * Convert the leading digit. - */ + /* Convert the leading digit */ i = 0; - err = mp_div(&b, &S, &dig, &b); + 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? - */ + /* 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)) { + mp_mul_2d(&b, 1, &b); + if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); } } else { - while (err == MP_OKAY) { - /* - * Shift by a group of digits. - */ + + for (;;) { + + /* 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); + 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); - } + mp_div_d(&S, dpow5[s5], &S, NULL); + 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); + mp_mul_d(&b, dpow5[g], &b); } + mp_mul_2d(&b, g, &b); /* - * As with the shortening bignum conversion, it's possible at this - * point that we will have reduced the denominator to less than - * 2**64/10, at which point it would be possible to fall back to - * to int64_t arithmetic. But the potential payoff is tremendously - * less - unless we're working in F format - because we know that - * three groups of digits will always suffice for %#.17e, the - * longest format that doesn't introduce empty precision. - * - * Extract the next group of digits. + * As with the shortening bignum conversion, it's possible at + * this point that we will have reduced the denominator to + * less than 2**64/10, at which point it would be possible to + * fall back to to int64 arithmetic. But the potential payoff + * is tremendously less - unless we're working in F format - + * because we know that three groups of digits will always + * suffice for %#.17e, the longest format that doesn't introduce + * empty precision. */ + /* Extract the next group of digits */ - if ((err != MP_OKAY) || (mp_div(&b, &S, &dig, &b) != MP_OKAY) || (dig.used > 1)) { + mp_div(&b, &S, &dig, &b); + if (dig.used > 1) { Tcl_Panic("wrong digit!"); } digit = dig.dp[0]; for (j = g-1; j >= 0; --j) { int t = itens[j]; - *s++ = digit / t + '0'; digit %= t; } i += g; - /* - * Have we converted all the requested digits? - */ + /* Have we converted all the requested digits? */ if (i == ilim) { - if ((mp_mul_2d(&b, 1, &b) == MP_OKAY) && ShouldBankerRoundUp(&b, &S, digit&1)) { + mp_mul_2d(&b, 1, &b); + if (ShouldBankerRoundUp(&b, &S, digit&1)) { s = BumpUp(s, retval, &k); + } else { + while (*--s == '0') { + /* do nothing */ + } + ++s; } - break; + 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); + mp_clear_multi(&b, &S, &temp, &dig, NULL); *s = '\0'; *decpt = k; if (endPtr) { *endPtr = s; } return retval; + } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * TclDoubleDigits -- * - * Core of Tcl's conversion of double-precision floating point numbers to - * decimal. + * 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. + * 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 + * This function is a service routine that produces the string of digits + * for floating-point-to-decimal conversion. It can do a number of things * according to the 'flags' argument. Valid values for 'flags' include: - * TCL_DD_SHORTEST - This is the default for floating point conversion if - * ::tcl_precision is 0. It constructs the shortest string of - * digits that will reconvert to the given number when scanned. + * TCL_DD_SHORTEST - This is the default for floating point conversion + * if ::tcl_precision is 0. It constructs the shortest string + * of digits that will reconvert to the given number when scanned. * For floating point numbers that are exactly between two * decimal numbers, it resolves using the 'round to even' rule. * With this value, the 'ndigits' parameter is ignored. - * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e format - * conversion (or for default floating->string if tcl_precision - * is not 0). It constructs a string of at most 'ndigits' digits, - * choosing the one that is closest to the given number (and - * resolving ties with 'round to even'). It is allowed to return - * fewer than 'ndigits' if the number converts exactly; if the - * TCL_DD_E_FORMAT|TCL_DD_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 + * TCL_DD_STEELE - This value is not recommended and may be removed + * in the future. It follows the conversion algorithm outlined + * in "How to Print Floating-Point Numbers Accurately" by + * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, + * pp. 112-126]. This rule has the effect of rendering 1e23 + * as 9.9999999999999999e22 - which is a 'better' approximation + * in the sense that it will reconvert correctly even if + * a subsequent input conversion is 'round up' or 'round down' + * rather than 'round to nearest', but is surprising otherwise. + * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e + * format conversion (or for default floating->string if + * tcl_precision is not 0). It constructs a string of at most + * 'ndigits' digits, choosing the one that is closest to the + * given number (and resolving ties with 'round to even'). + * It is allowed to return fewer than 'ndigits' if the number + * converts exactly; if the TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG + * is supplied instead, it also returns fewer digits if the + * shorter string will still reconvert to the given input number. + * In any case, strings of trailing zeroes are suppressed. + * TCL_DD_F_FORMAT - This value is used to prepare numbers for %f + * format conversion. It requests that conversion proceed until * 'ndigits' digits after the decimal point have been converted. * It is possible for this format to result in a zero-length - * string if the number is sufficiently small. Again, it is - * permissible for TCL_DD_F_FORMAT to return fewer digits for a - * number that converts exactly, and changing the argument to - * TCL_DD_F_FORMAT|TCL_DD_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. + * string if the number is sufficiently small. Again, it + * is permissible for TCL_DD_F_FORMAT to return fewer digits + * for a number that converts exactly, and changing the + * argument to TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow + * the routine also to return fewer digits if the shorter string + * will still reconvert without loss to the given input number. + * Strings of trailing zeroes are suppressed. + * + * To any of these flags may be OR'ed TCL_DD_NO_QUICK; this flag + * requires all calculations to be done in exact arithmetic. Normally, + * E and F format with fewer than about 14 digits will be done with + * a quick floating point approximation and fall back on the exact + * arithmetic only if the input number is close enough to the + * midpoint between two decimal strings that more precision is needed + * to resolve which string is correct. * * The value stored in the 'decpt' argument on return may be negative - * (indicating that the decimal point falls to the left of the string) or - * greater than the length of the string. In addition, the value -9999 is used - * as a sentinel to indicate that the string is one of the special values - * "Infinity" and "NaN", and that no decimal point should be inserted. + * (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. */ +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 convType = (flags & TCL_DD_CONVERSION_TYPE_MASK); + /* Type of conversion being performed + * TCL_DD_SHORTEST0 + * TCL_DD_STEELE0 + * TCL_DD_E_FORMAT + * TCL_DD_F_FORMAT */ + Double d; /* Union for deconstructing doubles */ + Tcl_WideUInt bw; /* Integer significand */ int be; /* Power of 2 by which b must be multiplied */ - int bbits; /* Number of bits needed to represent b. */ + int 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. */ + * 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. */ + * 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. - */ + /* Put the input number into a union for bit-whacking */ d.d = dv; @@ -4452,10 +4105,10 @@ TclDoubleDigits( /* * 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)). + * 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); @@ -4465,70 +4118,71 @@ TclDoubleDigits( * 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. + * 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. + * 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: + * 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. + * shortest and Steele, and ndigits for E format. + * ilim = The number of significant digits to convert if + * k has been guessed correctly. This is -1 for shortest and Steele + * (which stop when all significance has been lost), 'ndigits' + * for E format, and 'k + 1 + ndigits' for F format. + * ilim1 = The minimum number of significant digits to convert if + * k has been guessed 1 too high. This, too, is -1 for shortest + * and Steele, and 'ndigits' for E format, but it's 'ndigits-1' + * for F format. */ - SetPrecisionLimits(flags, k, &ndigits, &i, &ilim, &ilim1); + SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1); /* - * Try to do low-precision conversion in floating point rather than - * resorting to expensive multiprecision arithmetic. + * 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) { + if ((retval = QuickConversion(d.d, k, k_check, flags, + i, ilim, ilim1, + decpt, endPtr)) != NULL) { return retval; } } /* - * For shortening conversions, determine the upper and lower bounds for - * the remainder at which we can stop. - * m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the high - * side, and - * m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the low - * side. - * We may need to increase s2 to put m2plus, m2minus, b2 over a common - * denominator. + * 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) { + if (flags & TCL_DD_SHORTEN_FLAG) { int m2minus = b2; int m2plus; int m5 = b5; int len = i; /* - * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5) is 1/2 unit - * in the least significant place of the floating point number. + * 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 { @@ -4539,16 +4193,14 @@ TclDoubleDigits( /* * Reduce the fractions to lowest terms, since the above calculation - * may have left excess powers of 2 in numerator and denominator. + * 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; @@ -4556,56 +4208,60 @@ TclDoubleDigits( ++m2plus; } - if (s5+1 < N_LOG2POW5 && s2+1 + log2pow5[s5+1] < 64) { + 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]). + * 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); + return ShorteningInt64Conversion(&d, convType, bw, b2, b5, + m2plus, m2minus, m5, + s2, s5, k, len, ilim, ilim1, + decpt, endPtr); } else if (s5 == 0) { /* - * The denominator is a power of 2, so we can replace division by - * digit shifts. First we round up s2 to a multiple of MP_DIGIT_BIT, - * and adjust m2 and b2 accordingly. Then we launch into a version - * of the comparison that's specialized for the 'power of mp_digit - * in the denominator' case. + * The denominator is a power of 2, so we can replace division + * by digit shifts. First we round up s2 to a multiple of + * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch + * into a version of the comparison that's specialized for + * the 'power of mp_digit in the denominator' case. */ - - if (s2 % MP_DIGIT_BIT != 0) { - int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT); - + if (s2 % DIGIT_BIT != 0) { + int delta = DIGIT_BIT - (s2 % 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); + return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5, + m2plus, m2minus, m5, + s2/DIGIT_BIT, k, len, + ilim, ilim1, decpt, endPtr); } else { + /* - * Alas, there's no helpful special case; use full-up bignum - * arithmetic for the conversion. + * 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); + return ShorteningBignumConversion(&d, convType, bw, + b2, m2plus, m2minus, + s2, s5, k, len, + ilim, ilim1, decpt, endPtr); + } + } else { - /* - * Non-shortening conversion. - */ + + /* Non-shortening conversion */ int len = i; - /* - * Reduce numerator and denominator to lowest terms. - */ + /* Reduce numerator and denominator to lowest terms */ if (b2 >= s2 && s2 > 0) { b2 -= s2; s2 = 0; @@ -4613,46 +4269,48 @@ TclDoubleDigits( s2 -= b2; b2 = 0; } - if (s5+1 < N_LOG2POW5 && s2+1 + log2pow5[s5+1] < 64) { + 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. + * 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); + return StrictInt64Conversion(&d, convType, bw, b2, b5, + s2, s5, k, len, ilim, ilim1, + decpt, endPtr); + } else if (s5 == 0) { /* - * The denominator is a power of 2, so we can replace division by - * digit shifts. First we round up s2 to a multiple of MP_DIGIT_BIT, - * and adjust m2 and b2 accordingly. Then we launch into a version - * of the comparison that's specialized for the 'power of mp_digit - * in the denominator' case. + * The denominator is a power of 2, so we can replace division + * by digit shifts. First we round up s2 to a multiple of + * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch + * into a version of the comparison that's specialized for + * the 'power of mp_digit in the denominator' case. */ - - if (s2 % MP_DIGIT_BIT != 0) { - int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT); - + if (s2 % DIGIT_BIT != 0) { + int delta = DIGIT_BIT - (s2 % DIGIT_BIT); b2 += delta; s2 += delta; } - return StrictBignumConversionPowD(bw, b2, b5, - s2/MP_DIGIT_BIT, k, len, ilim, ilim1, decpt, endPtr); + return StrictBignumConversionPowD(&d, convType, bw, b2, b5, + s2/DIGIT_BIT, k, len, + ilim, ilim1, decpt, endPtr); } else { /* - * There are no helpful special cases, but at least we know in - * advance how many digits we will convert. We can run the - * conversion in steps of DIGIT_GROUP digits, so as to have many - * fewer mp_int divisions. + * 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); + return StrictBignumConversion(&d, convType, bw, b2, s2, s5, + k, len, ilim, ilim1, decpt, endPtr); } } } + /* *---------------------------------------------------------------------- @@ -4680,13 +4338,14 @@ TclInitDoubleConversion(void) 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; @@ -4711,8 +4370,8 @@ TclInitDoubleConversion(void) pow10_wide[i] = u; /* - * Determine how many bits of precision a double has, and how many decimal - * digits that represents. + * 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) { @@ -4723,8 +4382,8 @@ TclInitDoubleConversion(void) d = 1.0; /* - * Initialize a table of powers of ten that can be exactly represented in - * a double. + * 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)); @@ -4743,19 +4402,16 @@ TclInitDoubleConversion(void) */ for (i=0; i<9; ++i) { - err = err || mp_init(pow5 + i); + mp_init(pow5 + i); } - mp_set_u64(pow5, 5); + mp_set(pow5, 5); for (i=0; i<8; ++i) { - err = err || mp_sqr(pow5+i, pow5+i+1); + mp_sqr(pow5+i, pow5+i+1); } - err = err || mp_init_u64(pow5_13, 1220703125); + mp_init_set_int(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"); + mp_init(pow5_13 + i); + mp_sqr(pow5_13 + i - 1, pow5_13 + i); } /* @@ -4769,7 +4425,7 @@ TclInitDoubleConversion(void) + 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.)); + log10_DIGIT_MAX = (int) floor(DIGIT_BIT * log(2.) / log(10.)); /* * Nokia 770's software-emulated floating point is "middle endian": the @@ -4781,9 +4437,9 @@ TclInitDoubleConversion(void) #ifdef IEEE_FLOATING_POINT bitwhack.dv = 1.000000238418579; /* 3ff0 0000 4000 0000 */ - if ((bitwhack.iv >> 32) == 0x3FF00000) { + if ((bitwhack.iv >> 32) == 0x3ff00000) { n770_fp = 0; - } else if ((bitwhack.iv & 0xFFFFFFFF) == 0x3FF00000) { + } else if ((bitwhack.iv & 0xffffffff) == 0x3ff00000) { n770_fp = 1; } else { Tcl_Panic("unknown floating point word order on this machine"); @@ -4812,13 +4468,10 @@ TclFinalizeDoubleConversion(void) { int i; - ckfree(pow10_wide); + ckfree((char *) pow10_wide); for (i=0; i<9; ++i) { mp_clear(pow5 + i); } - for (i=0; i < 5; ++i) { - mp_clear(pow5_13 + i); - } } /* @@ -4841,49 +4494,42 @@ TclFinalizeDoubleConversion(void) int Tcl_InitBignumFromDouble( - Tcl_Interp *interp, /* For error message. */ - double d, /* Number to convert. */ - void *big) /* Place to store the result. */ + Tcl_Interp *interp, /* For error message */ + double d, /* Number to convert */ + mp_int *b) /* 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 (TclIsInfinite(d)) { if (interp != NULL) { const char *s = "integer value too large to represent"; Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1)); - Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, (void *)NULL); + Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL); } return TCL_ERROR; } - fract = frexp(d, &expt); + fract = frexp(d,&expt); if (expt <= 0) { - err = mp_init(b); + mp_init(b); mp_zero(b); } else { - Tcl_WideInt w = (Tcl_WideInt)ldexp(fract, mantBits); + 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); + TclBNInitBignumFromWideInt(b, w); + if (shift < 0) { + mp_div_2d(b, -shift, b, NULL); } else if (shift > 0) { - err = mp_mul_2d(b, shift, b); + mp_mul_2d(b, shift, b); } } - if (err != MP_OKAY) { - return TCL_ERROR; - } return TCL_OK; } @@ -4904,13 +4550,11 @@ Tcl_InitBignumFromDouble( double TclBignumToDouble( - const void *big) /* Integer to convert. */ + mp_int *a) /* Integer to convert. */ { mp_int b; int bits, shift, i, lsb; double r; - mp_err err; - const mp_int *a = (const mp_int *)big; /* @@ -4921,10 +4565,10 @@ TclBignumToDouble( bits = mp_count_bits(a); if (bits > DBL_MAX_EXP*log2FLT_RADIX) { errno = ERANGE; - if (mp_isneg(a)) { - return -HUGE_VAL; - } else { + if (a->sign == MP_ZPOS) { return HUGE_VAL; + } else { + return -HUGE_VAL; } } shift = mantBits - bits; @@ -4939,13 +4583,11 @@ TclBignumToDouble( * 'rounded to even'. */ - err = mp_init(&b); - if (err != MP_OKAY) { - /* just skip */ - } else if (shift == 0) { - err = mp_copy(a, &b); + mp_init(&b); + if (shift == 0) { + mp_copy(a, &b); } else if (shift > 0) { - err = mp_mul_2d(a, shift, &b); + mp_mul_2d(a, shift, &b); } else if (shift < 0) { lsb = mp_cnt_lsb(a); if (lsb == -1-shift) { @@ -4954,12 +4596,12 @@ TclBignumToDouble( * 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); + mp_div_2d(a, -shift, &b, NULL); + if (mp_isodd(&b)) { + if (b.sign == MP_ZPOS) { + mp_add_d(&b, 1, &b); } else { - err = mp_add_d(&b, 1, &b); + mp_sub_d(&b, 1, &b); } } } else { @@ -4968,15 +4610,13 @@ TclBignumToDouble( * 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); + mp_div_2d(a, -1-shift, &b, NULL); + if (b.sign == MP_ZPOS) { + mp_add_d(&b, 1, &b); } else { - err = mp_add_d(&b, 1, &b); + mp_sub_d(&b, 1, &b); } - err = mp_div_2d(&b, 1, &b, NULL); + mp_div_2d(&b, 1, &b, NULL); } } @@ -4984,12 +4624,9 @@ TclBignumToDouble( * 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]; + for (i=b.used-1 ; i>=0 ; --i) { + r = ldexp(r, DIGIT_BIT) + b.dp[i]; } mp_clear(&b); @@ -5003,15 +4640,15 @@ TclBignumToDouble( * Return the result with the appropriate sign. */ - if (mp_isneg(a)) { - return -r; - } else { + if (a->sign == MP_ZPOS) { return r; + } else { + return -r; } } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * TclCeil -- * @@ -5021,21 +4658,19 @@ TclBignumToDouble( * Results: * Returns the floating point number. * - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- */ double TclCeil( - const void *big) /* Integer to convert. */ + mp_int *a) /* 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); + mp_init(&b); + if (mp_cmp_d(a, 0) == MP_LT) { + mp_neg(a, &b); r = -TclFloor(&b); } else { int bits = mp_count_bits(a); @@ -5045,29 +4680,22 @@ TclCeil( } 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); + if (shift > 0) { + 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); - } + mp_init(&d); + 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); + mp_copy(a, &b); } - if (err != MP_OKAY) { - return 0.0; + if (!exact) { + mp_add_d(&b, 1, &b); } for (i=b.used-1 ; i>=0 ; --i) { - r = ldexp(r, MP_DIGIT_BIT) + b.dp[i]; + r = ldexp(r, DIGIT_BIT) + b.dp[i]; } r = ldexp(r, bits - mantBits); } @@ -5075,33 +4703,31 @@ TclCeil( mp_clear(&b); return r; } - + /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * * TclFloor -- * - * Computes the largest floating point number less than or equal to the - * mp_int argument. + * 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. */ + mp_int *a) /* 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); + mp_init(&b); + if (mp_cmp_d(a, 0) == MP_LT) { + mp_neg(a, &b); r = -TclCeil(&b); } else { int bits = mp_count_bits(a); @@ -5112,17 +4738,14 @@ TclFloor( int i, shift = mantBits - bits; if (shift > 0) { - err = mp_mul_2d(a, shift, &b); + mp_mul_2d(a, shift, &b); } else if (shift < 0) { - err = mp_div_2d(a, -shift, &b, NULL); + mp_div_2d(a, -shift, &b, NULL); } else { - err = mp_copy(a, &b); - } - if (err != MP_OKAY) { - return 0.0; + mp_copy(a, &b); } for (i=b.used-1 ; i>=0 ; --i) { - r = ldexp(r, MP_DIGIT_BIT) + b.dp[i]; + r = ldexp(r, DIGIT_BIT) + b.dp[i]; } r = ldexp(r, bits - mantBits); } @@ -5153,15 +4776,14 @@ TclFloor( static double BignumToBiasedFrExp( - const mp_int *a, /* Integer to convert. */ - int *machexp) /* Power of two. */ + 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. @@ -5170,15 +4792,13 @@ BignumToBiasedFrExp( bits = mp_count_bits(a); shift = mantBits - 2 - bits; - if (mp_init(&b)) { - return 0.0; - } + mp_init(&b); if (shift > 0) { - err = mp_mul_2d(a, shift, &b); + mp_mul_2d(a, shift, &b); } else if (shift < 0) { - err = mp_div_2d(a, -shift, &b, NULL); + mp_div_2d(a, -shift, &b, NULL); } else { - err = mp_copy(a, &b); + mp_copy(a, &b); } /* @@ -5186,10 +4806,8 @@ BignumToBiasedFrExp( */ r = 0.0; - if (err == MP_OKAY) { - for (i=b.used-1; i>=0; --i) { - r = ldexp(r, MP_DIGIT_BIT) + b.dp[i]; - } + for (i=b.used-1; i>=0; --i) { + r = ldexp(r, DIGIT_BIT) + b.dp[i]; } mp_clear(&b); @@ -5198,7 +4816,7 @@ BignumToBiasedFrExp( */ *machexp = bits - mantBits + 2; - return (mp_isneg(a) ? -r : r); + return ((a->sign == MP_ZPOS) ? r : -r); } /* @@ -5223,8 +4841,8 @@ BignumToBiasedFrExp( static double Pow10TimesFrExp( - int exponent, /* Power of 10 to multiply by. */ - double fraction, /* Significand of multiplicand. */ + 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. */ { @@ -5234,10 +4852,10 @@ Pow10TimesFrExp( if (exponent > 0) { /* - * Multiply by 10**exponent. + * Multiply by 10**exponent */ - retval = frexp(retval * pow10vals[exponent & 0xF], &j); + retval = frexp(retval * pow10vals[exponent&0xf], &j); expt += j; for (i=4; i<9; ++i) { if (exponent & (1<<i)) { @@ -5247,10 +4865,10 @@ Pow10TimesFrExp( } } else if (exponent < 0) { /* - * Divide by 10**-exponent. + * Divide by 10**-exponent */ - retval = frexp(retval / pow10vals[(-exponent) & 0xF], &j); + retval = frexp(retval / pow10vals[(-exponent) & 0xf], &j); expt += j; for (i=4; i<9; ++i) { if ((-exponent) & (1<<i)) { @@ -5328,23 +4946,23 @@ TclFormatNaN( #else union { double dv; - uint64_t iv; + Tcl_WideUInt 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); + if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { + bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63); *buffer++ = '-'; } *buffer++ = 'N'; *buffer++ = 'a'; *buffer++ = 'N'; - bitwhack.iv &= ((UINT64_C(1)) << 51) - 1; + bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1; if (bitwhack.iv != 0) { - snprintf(buffer, TCL_DOUBLE_SPACE, "(%" PRIx64 ")", bitwhack.iv); + sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv); } else { *buffer = '\0'; } @@ -5356,27 +4974,26 @@ TclFormatNaN( * * Nokia770Twiddle -- * - * Transpose the two words of a number for Nokia 770 floating point - * handling. + * 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. */ + Tcl_WideUInt w) /* Number to transpose */ { - return (((w >> 32) & 0xFFFFFFFF) | (w << 32)); + return (((w >> 32) & 0xffffffff) | (w << 32)); } -#endif /* *---------------------------------------------------------------------- * * TclNokia770Doubles -- * - * Transpose the two words of a number for Nokia 770 floating point - * handling. + * Transpose the two words of a number for Nokia 770 floating + * point handling. * *---------------------------------------------------------------------- */ |