summaryrefslogtreecommitdiffstats
path: root/generic/tclStrToD.c
diff options
context:
space:
mode:
Diffstat (limited to 'generic/tclStrToD.c')
-rwxr-xr-xgeneric/tclStrToD.c2772
1 files changed, 2474 insertions, 298 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c
index d0a5345..8171df0 100755
--- a/generic/tclStrToD.c
+++ b/generic/tclStrToD.c
@@ -14,7 +14,7 @@
* See the file "license.terms" for information on usage and redistribution of
* this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclStrToD.c,v 1.46 2010/05/21 12:43:29 nijtmans Exp $
+ * RCS: @(#) $Id: tclStrToD.c,v 1.47 2010/11/28 23:20:11 kennykb Exp $
*
*----------------------------------------------------------------------
*/
@@ -90,6 +90,66 @@ typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
* runtime).
*/
+/* Magic constants */
+
+#define LOG10_2 0.3010299956639812
+#define TWO_OVER_3LOG10 0.28952965460216784
+#define LOG10_3HALVES_PLUS_FUDGE 0.1760912590558
+
+/* Definitions of the parts of an IEEE754-format floating point number */
+
+#define SIGN_BIT 0x80000000
+ /* Mask for the sign bit in the first
+ * word of a double */
+#define EXP_MASK 0x7ff00000
+ /* Mask for the exponent field in the
+ * first word of a double */
+#define EXP_SHIFT 20
+ /* Shift count to make the exponent an
+ * integer */
+#define HIDDEN_BIT (((Tcl_WideUInt) 0x00100000) << 32)
+ /* Hidden 1 bit for the significand */
+#define HI_ORDER_SIG_MASK 0x000fffff
+ /* Mask for the high-order part of the
+ * significand in the first word of a
+ * double */
+#define SIG_MASK (((Tcl_WideUInt) HI_ORDER_SIG_MASK << 32) \
+ | 0xffffffff)
+ /* Mask for the 52-bit significand. */
+#define FP_PRECISION 53
+ /* Number of bits of significand plus the
+ * hidden bit */
+#define EXPONENT_BIAS 0x3ff
+ /* Bias of the exponent 0 */
+
+/* Derived quantities */
+
+#define TEN_PMAX 22
+ /* floor(FP_PRECISION*log(2)/log(5)) */
+#define QUICK_MAX 14
+ /* floor((FP_PRECISION-1)*log(2)/log(10)) - 1 */
+#define BLETCH 0x10
+ /* Highest power of two that is greater than
+ * DBL_MAX_10_EXP, divided by 16 */
+#define DIGIT_GROUP 8
+ /* floor(DIGIT_BIT*log(2)/log(10)) */
+
+/* Union used to dismantle floating point numbers. */
+
+typedef union Double {
+ struct {
+#ifdef WORDS_BIGENDIAN
+ int word0;
+ int word1;
+#else
+ int word1;
+ int word0;
+#endif
+ } w;
+ double d;
+ Tcl_WideUInt q;
+} Double;
+
static int maxpow10_wide; /* The powers of ten that can be represented
* exactly as wide integers. */
static Tcl_WideUInt *pow10_wide;
@@ -123,6 +183,7 @@ static const double pow_10_2_n[] = { /* Inexact higher powers of ten. */
1.0e+128,
1.0e+256
};
+
static int n770_fp; /* Flag is 1 on Nokia N770 floating point.
* Nokia's floating point has the words
* reversed: if big-endian is 7654 3210,
@@ -131,27 +192,161 @@ static int n770_fp; /* Flag is 1 on Nokia N770 floating point.
* little-endian within the 32-bit words
* but big-endian between them. */
+/* Table of powers of 5 that are small enough to fit in an mp_digit. */
+
+static const mp_digit dpow5[13] = {
+ 1, 5, 25, 125,
+ 625, 3125, 15625, 78125,
+ 390625, 1953125, 9765625, 48828125,
+ 244140625
+};
+
+/* Table of powers: pow5_13[n] = 5**(13*2**(n+1)) */
+static mp_int pow5_13[5]; /* Table of powers: 5**13, 5**26, 5**52,
+ * 5**104, 5**208 */
+static const double tens[] = {
+ 1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09,
+ 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
+ 1e20, 1e21, 1e22
+};
+
+static const int itens [] = {
+ 1,
+ 10,
+ 100,
+ 1000,
+ 10000,
+ 100000,
+ 1000000,
+ 10000000,
+ 100000000
+};
+
+static const Tcl_WideUInt wtens[] = {
+ 1, 10, 100, 1000, 10000, 100000, 1000000,
+ (Tcl_WideUInt) 1000000*10, (Tcl_WideUInt) 1000000*100,
+ (Tcl_WideUInt) 1000000*1000, (Tcl_WideUInt) 1000000*10000,
+ (Tcl_WideUInt) 1000000*100000, (Tcl_WideUInt) 1000000*1000000,
+ (Tcl_WideUInt) 1000000*1000000*10, (Tcl_WideUInt) 1000000*1000000*100,
+ (Tcl_WideUInt) 1000000*1000000*1000,(Tcl_WideUInt) 1000000*1000000*10000
+
+};
+
+static const double bigtens[] = {
+ 1e016, 1e032, 1e064, 1e128, 1e256
+};
+#define N_BIGTENS 5
+
+static const int log2pow5[27] = {
+ 01, 3, 5, 7, 10, 12, 14, 17, 19, 21,
+ 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
+ 47, 49, 52, 54, 56, 59, 61
+};
+#define N_LOG2POW5 27
+
+static const Tcl_WideUInt wuipow5[27] = {
+ (Tcl_WideUInt) 1, /* 5**0 */
+ (Tcl_WideUInt) 5,
+ (Tcl_WideUInt) 25,
+ (Tcl_WideUInt) 125,
+ (Tcl_WideUInt) 625,
+ (Tcl_WideUInt) 3125, /* 5**5 */
+ (Tcl_WideUInt) 3125*5,
+ (Tcl_WideUInt) 3125*25,
+ (Tcl_WideUInt) 3125*125,
+ (Tcl_WideUInt) 3125*625,
+ (Tcl_WideUInt) 3125*3125, /* 5**10 */
+ (Tcl_WideUInt) 3125*3125*5,
+ (Tcl_WideUInt) 3125*3125*25,
+ (Tcl_WideUInt) 3125*3125*125,
+ (Tcl_WideUInt) 3125*3125*625,
+ (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */
+ (Tcl_WideUInt) 3125*3125*3125*5,
+ (Tcl_WideUInt) 3125*3125*3125*25,
+ (Tcl_WideUInt) 3125*3125*3125*125,
+ (Tcl_WideUInt) 3125*3125*3125*625,
+ (Tcl_WideUInt) 3125*3125*3125*3125, /* 5**20 */
+ (Tcl_WideUInt) 3125*3125*3125*3125*5,
+ (Tcl_WideUInt) 3125*3125*3125*3125*25,
+ (Tcl_WideUInt) 3125*3125*3125*3125*125,
+ (Tcl_WideUInt) 3125*3125*3125*3125*625,
+ (Tcl_WideUInt) 3125*3125*3125*3125*3125, /* 5**25 */
+ (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */
+};
+
/*
* Static functions defined in this file.
*/
-static double AbsoluteValue(double v, int *signum);
static int AccumulateDecimalDigit(unsigned, int,
Tcl_WideUInt *, mp_int *, int);
-static double BignumToBiasedFrExp(const mp_int *big, int *machexp);
-static int GetIntegerTimesPower(double v, mp_int *r, int *e);
static double MakeHighPrecisionDouble(int signum,
mp_int *significand, int nSigDigs, int exponent);
static double MakeLowPrecisionDouble(int signum,
Tcl_WideUInt significand, int nSigDigs,
int exponent);
static double MakeNaN(int signum, Tcl_WideUInt tag);
-static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
-static double Pow10TimesFrExp(int exponent, double fraction,
- int *machexp);
static double RefineApproximation(double approx,
mp_int *exactSignificand, int exponent);
+static void MulPow5(mp_int*, unsigned, mp_int*);
+static int NormalizeRightward(Tcl_WideUInt*);
+static int RequiredPrecision(Tcl_WideUInt);
+static void DoubleToExpAndSig(double, Tcl_WideUInt*, int*, int*);
+static void TakeAbsoluteValue(Double*, int*);
+static char* FormatInfAndNaN(Double*, int*, char**);
+static char* FormatZero(int*, char**);
+static int ApproximateLog10(Tcl_WideUInt, int, int);
+static int BetterLog10(double, int, int*);
+static void ComputeScale(int, int, int*, int*, int*, int*);
+static void SetPrecisionLimits(int, int, int*, int*, int*, int*);
+static char* BumpUp(char*, char*, int*);
+static int AdjustRange(double*, int);
+static char* ShorteningQuickFormat(double, int, int, double,
+ char*, int*);
+static char* StrictQuickFormat(double, int, int, double,
+ char*, int*);
+static char* QuickConversion(double, int, int, int, int, int, int,
+ int*, char**);
+static void CastOutPowersOf2(int*, int*, int*);
+static char* ShorteningInt64Conversion(Double*, int, Tcl_WideUInt,
+ int, int, int, int, int, int, int, int, int,
+ int, int, int*, char**);
+static char* StrictInt64Conversion(Double*, int, Tcl_WideUInt,
+ int, int, int, int, int, int,
+ int, int, int*, char**);
+static int ShouldBankerRoundUpPowD(mp_int*, int, int);
+static int ShouldBankerRoundUpToNextPowD(mp_int*, mp_int*,
+ int, int, int, mp_int*);
+static char* ShorteningBignumConversionPowD(Double* dPtr,
+ int convType, Tcl_WideUInt bw, int b2, int b5,
+ int m2plus, int m2minus, int m5,
+ int sd, int k, int len,
+ int ilim, int ilim1, int* decpt,
+ char** endPtr);
+static char* StrictBignumConversionPowD(Double* dPtr, int convType,
+ Tcl_WideUInt bw, int b2, int b5,
+ int sd, int k, int len,
+ int ilim, int ilim1, int* decpt,
+ char** endPtr);
+static int ShouldBankerRoundUp(mp_int*, mp_int*, int);
+static int ShouldBankerRoundUpToNext(mp_int*, mp_int*, mp_int*,
+ int, int, mp_int*);
+static char* ShorteningBignumConversion(Double* dPtr, int convType,
+ Tcl_WideUInt bw, int b2,
+ int m2plus, int m2minus,
+ int s2, int s5, int k, int len,
+ int ilim, int ilim1, int* decpt,
+ char** endPtr);
+static char* StrictBignumConversion(Double* dPtr, int convType,
+ Tcl_WideUInt bw, int b2,
+ int s2, int s5, int k, int len,
+ int ilim, int ilim1, int* decpt,
+ char** endPtr);
+static double BignumToBiasedFrExp(const mp_int *big, int *machexp);
+static double Pow10TimesFrExp(int exponent, double fraction,
+ int *machexp);
static double SafeLdExp(double fraction, int exponent);
+static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
/*
*----------------------------------------------------------------------
@@ -1732,414 +1927,2362 @@ RefineApproximation(
}
/*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
*
- * TclDoubleDigits --
+ * MultPow5 --
+ *
+ * Multiply a bignum by a power of 5.
+ *
+ * Side effects:
+ * Stores base*5**n in result
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+MulPow5(mp_int* base, /* Number to multiply */
+ unsigned n, /* Power of 5 to multiply by */
+ mp_int* result) /* Place to store the result */
+{
+ mp_int* p = base;
+ int n13 = n / 13;
+ int r = n % 13;
+ if (r != 0) {
+ mp_mul_d(p, dpow5[r], result);
+ p = result;
+ }
+ r = 0;
+ while (n13 != 0) {
+ if (n13 & 1) {
+ mp_mul(p, pow5_13+r, result);
+ p = result;
+ }
+ n13 >>= 1;
+ ++r;
+ }
+ if (p != result) {
+ mp_copy(p, result);
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * NormalizeRightward --
*
- * Converts a double to a string of digits.
+ * Shifts a number rightward until it is odd (that is, until the
+ * least significant bit is nonzero.
*
* Results:
- * Returns the position of the character in the string after which the
- * decimal point should appear. Since the string contains only
- * significant digits, the position may be less than zero or greater than
- * the length of the string.
+ * Returns the number of bit positions by which the number was shifted.
*
* Side effects:
- * Stores the digits in the given buffer and sets 'signum' according to
- * the sign of the number.
+ * Shifts the number in place; *wPtr is replaced by the shifted number.
*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
+ */
+inline static int
+NormalizeRightward(Tcl_WideUInt* wPtr)
+ /* INOUT: Number to shift */
+{
+ int rv = 0;
+ Tcl_WideUInt w = *wPtr;
+ if (!(w & (Tcl_WideUInt) 0xffffffff)) {
+ w >>= 32; rv += 32;
+ }
+ if (!(w & (Tcl_WideUInt) 0xffff)) {
+ w >>= 16; rv += 16;
+ }
+ if (!(w & (Tcl_WideUInt) 0xff)) {
+ w >>= 8; rv += 8;
+ }
+ if (!(w & (Tcl_WideUInt) 0xf)) {
+ w >>= 4; rv += 4;
+ }
+ if (!(w & 0x3)) {
+ w >>= 2; rv += 2;
+ }
+ if (!(w & 0x1)) {
+ w >>= 1; ++rv;
+ }
+ *wPtr = w;
+ return rv;
+}
+
+/*
+ *-----------------------------------------------------------------------------0
+ *
+ * RequiredPrecision --
+ *
+ * Determines the number of bits needed to hold an intger.
+ *
+ * Results:
+ * Returns the position of the most significant bit (0 - 63).
+ * Returns 0 if the number is zero.
+ *
+ *----------------------------------------------------------------------------
*/
-int
-TclDoubleDigits(
- char *buffer, /* Buffer in which to store the result, must
- * have at least 18 chars */
- double v, /* Number to convert. Must be finite, and not
- * NaN */
- int *signum) /* Output: 1 if the number is negative.
- * Should handle -0 correctly on the IEEE
- * architecture. */
+static int
+RequiredPrecision(Tcl_WideUInt w)
+ /* Number to interrogate */
{
- int e; /* Power of FLT_RADIX that satisfies
- * v = f * FLT_RADIX**e */
- int lowOK, highOK;
- mp_int r; /* Scaled significand. */
- mp_int s; /* Divisor such that v = r / s */
- int smallestSig; /* Flag == 1 iff v's significand is the
- * smallest that can be represented. */
- mp_int mplus; /* Scaled epsilon: (r + 2* mplus) == v(+)
- * where v(+) is the floating point successor
- * of v. */
- mp_int mminus; /* Scaled epsilon: (r - 2*mminus) == v(-)
- * where v(-) is the floating point
- * predecessor of v. */
- mp_int temp;
- int rfac2 = 0; /* Powers of 2 and 5 by which large */
- int rfac5 = 0; /* integers should be scaled. */
- int sfac2 = 0;
- int sfac5 = 0;
- int mplusfac2 = 0;
- int mminusfac2 = 0;
- char c;
- int i, k, n;
+ int rv;
+ unsigned long wi;
+ if (w & ((Tcl_WideUInt) 0xffffffff << 32)) {
+ wi = w >> 32; rv = 32;
+ } else {
+ wi = w; rv = 0;
+ }
+ if (wi & 0xffff0000) {
+ wi >>= 16; rv += 16;
+ }
+ if (wi & 0xff00) {
+ wi >>= 8; rv += 8;
+ }
+ if (wi & 0xf0) {
+ wi >>= 4; rv += 4;
+ }
+ if (wi & 0xc) {
+ wi >>= 2; rv += 2;
+ }
+ if (wi & 0x2) {
+ wi >>= 1; ++rv;
+ }
+ if (wi & 0x1) {
+ ++rv;
+ }
+ return rv;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * DoubleToExpAndSig --
+ *
+ * Separates a 'double' into exponent and significand.
+ *
+ * Side effects:
+ * Stores the significand in '*significand' and the exponent in
+ * '*expon' so that dv == significand * 2.0**expon, and significand
+ * is odd. Also stores the position of the leftmost 1-bit in 'significand'
+ * in 'bits'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+DoubleToExpAndSig(double dv, /* Number to convert */
+ Tcl_WideUInt* significand,
+ /* OUTPUT: Significand of the number */
+ int* expon, /* OUTPUT: Exponent to multiply the number by */
+ int* bits) /* OUTPUT: Number of significant bits */
+{
+ Double d; /* Number being converted */
+ Tcl_WideUInt z; /* Significand under construction */
+ int de; /* Exponent of the number */
+ int k; /* Bit count */
+
+ d.d = dv;
+
+ /* Extract exponent and significand */
+
+ de = (d.w.word0 & EXP_MASK) >> EXP_SHIFT;
+ z = d.q & SIG_MASK;
+ if (de != 0) {
+ z |= HIDDEN_BIT;
+ k = NormalizeRightward(&z);
+ *bits = FP_PRECISION - k;
+ *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1);
+ } else {
+ k = NormalizeRightward(&z);
+ *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1) + 1;
+ *bits = RequiredPrecision(z);
+ }
+ *significand = z;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TakeAbsoluteValue --
+ *
+ * Takes the absolute value of a 'double' including 0, Inf and NaN
+ *
+ * Side effects:
+ * The 'double' in *d is replaced with its absolute value. The
+ * signum is stored in 'sign': 1 for negative, 0 for nonnegative.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+TakeAbsoluteValue(Double* d, /* Number to replace with absolute value */
+ int* sign) /* Place to put the signum */
+{
+ if (d->w.word0 & SIGN_BIT) {
+ *sign = 1;
+ d->w.word0 &= ~SIGN_BIT;
+ } else {
+ *sign = 0;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * FormatInfAndNaN --
+ *
+ * Bailout for formatting infinities and Not-A-Number.
+ *
+ * Results:
+ * Returns one of the strings 'Infinity' and 'NaN'.
+ *
+ * Side effects:
+ * Stores 9999 in *decpt, and sets '*endPtr' to designate the
+ * terminating NUL byte of the string if 'endPtr' is not NULL.
+ *
+ * The string returned must be freed by the caller using 'ckfree'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+FormatInfAndNaN(Double* d, /* Exceptional number to format */
+ int* decpt, /* Decimal point to set to a bogus value */
+ char** endPtr) /* Pointer to the end of the formatted data */
+{
+ char* retval;
+ *decpt = 9999;
+ if (!(d->w.word1) && !(d->w.word0 & HI_ORDER_SIG_MASK)) {
+ retval = ckalloc(9);
+ strcpy(retval, "Infinity");
+ if (endPtr) {
+ *endPtr = retval + 8;
+ }
+ } else {
+ retval = ckalloc(4);
+ strcpy(retval, "NaN");
+ if (endPtr) {
+ *endPtr = retval + 3;
+ }
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * FormatZero --
+ *
+ * Bailout to format a zero floating-point number.
+ *
+ * Results:
+ * Returns the constant string "0"
+ *
+ * Side effects:
+ * Stores 1 in '*decpt' and puts a pointer to the NUL byte terminating
+ * the string in '*endPtr' if 'endPtr' is not NULL.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+FormatZero(int* decpt, /* Location of the decimal point */
+ char** endPtr) /* Pointer to the end of the formatted data */
+{
+ char* retval = ckalloc(2);
+ strcpy(retval, "0");
+ if (endPtr) {
+ *endPtr = retval+1;
+ }
+ *decpt = 0;
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ApproximateLog10 --
+ *
+ * Computes a two-term Taylor series approximation to the common
+ * log of a number, and computes the number's binary log.
+ *
+ * Results:
+ * Return an approximation to floor(log10(bw*2**be)) that is either
+ * exact or 1 too high.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ApproximateLog10(Tcl_WideUInt bw,
+ /* Integer significand of the number */
+ int be, /* Power of two to scale bw */
+ int bbits) /* Number of bits of precision in bw */
+{
+ int i; /* Log base 2 of the number */
+ int k; /* Floor(Log base 10 of the number) */
+ double ds; /* Mantissa of the number */
+ Double d2;
/*
- * Split the number into absolute value and signum.
+ * Compute i and d2 such that d = d2*2**i, and 1 < d2 < 2.
+ * Compute an approximation to log10(d),
+ * log10(d) ~ log10(2) * i + log10(1.5)
+ * + (significand-1.5)/(1.5 * log(10))
*/
- v = AbsoluteValue(v, signum);
+ d2.q = bw << (FP_PRECISION - bbits) & SIG_MASK;
+ d2.w.word0 |= (EXPONENT_BIAS) << EXP_SHIFT;
+ i = be + bbits - 1;
+ ds = (d2.d - 1.5) * TWO_OVER_3LOG10
+ + LOG10_3HALVES_PLUS_FUDGE
+ + LOG10_2 * i;
+ k = (int) ds;
+ if (k > ds) {
+ --k;
+ }
+ return k;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * BetterLog10 --
+ *
+ * Improves the result of ApproximateLog10 for numbers in the range
+ * 1 .. 10**(TEN_PMAX)-1
+ *
+ * Side effects:
+ * Sets k_check to 0 if the new result is known to be exact, and to
+ * 1 if it may still be one too high.
+ *
+ * Results:
+ * Returns the improved approximation to log10(d)
+ *
+ *-----------------------------------------------------------------------------
+ */
- /*
- * Handle zero specially.
+inline static int
+BetterLog10(double d, /* Original number to format */
+ int k, /* Characteristic(Log base 10) of the number */
+ int* k_check) /* Flag == 1 if k is inexact */
+{
+ /*
+ * Performance hack. If k is in the range 0..TEN_PMAX, then we can
+ * use a powers-of-ten table to check it.
*/
+ if (k >= 0 && k <= TEN_PMAX) {
+ if (d < tens[k]) {
+ k--;
+ }
+ *k_check = 0;
+ } else {
+ *k_check = 1;
+ }
+ return k;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ComputeScale --
+ *
+ * Prepares to format a floating-point number as decimal.
+ *
+ * Parameters:
+ * floor(log10*x) is k (or possibly k-1). floor(log2(x) is i.
+ * The significand of x requires bbits bits to represent.
+ *
+ * Results:
+ * Determines integers b2, b5, s2, s5 so that sig*2**b2*5**b5/2**s2*2**s5
+ * exactly represents the value of the x/10**k. This value will lie
+ * in the range [1 .. 10), and allows for computing successive digits
+ * by multiplying sig%10 by 10.
+ *
+ *-----------------------------------------------------------------------------
+ */
- if (v == 0.0) {
- *buffer++ = '0';
- *buffer++ = '\0';
- return 1;
+inline static void
+ComputeScale(int be, /* Exponent part of number: d = bw * 2**be */
+ int k, /* Characteristic of log10(number) */
+ int* b2, /* OUTPUT: Power of 2 in the numerator */
+ int* b5, /* OUTPUT: Power of 5 in the numerator */
+ int* s2, /* OUTPUT: Power of 2 in the denominator */
+ int* s5) /* OUTPUT: Power of 5 in the denominator */
+{
+
+ /*
+ * Scale numerator and denominator powers of 2 so that the
+ * input binary number is the ratio of integers
+ */
+ if (be <= 0) {
+ *b2 = 0;
+ *s2 = -be;
+ } else {
+ *b2 = be;
+ *s2 = 0;
}
+ /*
+ * Scale numerator and denominator so that the output decimal number
+ * is the ratio of integers
+ */
+ if (k >= 0) {
+ *b5 = 0;
+ *s5 = k;
+ *s2 += k;
+ } else {
+ *b2 -= k;
+ *b5 = -k;
+ *s5 = 0;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * SetPrecisionLimits --
+ *
+ * Determines how many digits of significance should be computed
+ * (and, hence, how much memory need be allocated) for formatting a
+ * floating point number.
+ *
+ * Given that 'k' is floor(log10(x)):
+ * if 'shortest' format is used, there will be at most 18 digits in the result.
+ * if 'F' format is used, there will be at most 'ndigits' + k + 1 digits
+ * if 'E' format is used, there will be exactly 'ndigits' digits.
+ *
+ * Side effects:
+ * Adjusts '*ndigitsPtr' to have a valid value.
+ * Stores the maximum memory allocation needed in *iPtr.
+ * Sets '*iLimPtr' to the limiting number of digits to convert if k
+ * has been guessed correctly, and '*iLim1Ptr' to the limiting number
+ * of digits to convert if k has been guessed to be one too high.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static void
+SetPrecisionLimits(int convType,
+ /* Type of conversion:
+ * TCL_DD_SHORTEST
+ * TCL_DD_STEELE0
+ * TCL_DD_E_FMT
+ * TCL_DD_F_FMT */
+ int k, /* Floor(log10(number to convert)) */
+ int* ndigitsPtr,
+ /* IN/OUT: Number of digits requested
+ * (Will be adjusted if needed) */
+ int* iPtr, /* OUT: Maximum number of digits
+ * to return */
+ int *iLimPtr,/* OUT: Number of digits of significance
+ * if the bignum method is used.*/
+ int *iLim1Ptr)
+ /* OUT: Number of digits of significance
+ * if the quick method is used. */
+{
+ switch(convType) {
+ case TCL_DD_SHORTEST0:
+ case TCL_DD_STEELE0:
+ *iLimPtr = *iLim1Ptr = -1;
+ *iPtr = 18;
+ *ndigitsPtr = 0;
+ break;
+ case TCL_DD_E_FORMAT:
+ if (*ndigitsPtr <= 0) {
+ *ndigitsPtr = 1;
+ }
+ *iLimPtr = *iLim1Ptr = *iPtr = *ndigitsPtr;
+ break;
+ case TCL_DD_F_FORMAT:
+ *iPtr = *ndigitsPtr + k + 1;
+ *iLimPtr = *iPtr;
+ *iLim1Ptr = *iPtr - 1;
+ if (*iPtr <= 0) {
+ *iPtr = 1;
+ }
+ break;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * BumpUp --
+ *
+ * Increases a string of digits ending in a series of nines to
+ * designate the next higher number. xxxxb9999... -> xxxx(b+1)0000...
+ *
+ * Results:
+ * Returns a pointer to the end of the adjusted string.
+ *
+ * Side effects:
+ * In the case that the string consists solely of '999999', sets it
+ * to "1" and moves the decimal point (*kPtr) one place to the right.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+
+inline static char*
+BumpUp(char* s, /* Cursor pointing one past the end of the
+ * string */
+ char* retval, /* Start of the string of digits */
+ int* kPtr) /* Position of the decimal point */
+{
+ while (*--s == '9') {
+ if (s == retval) {
+ ++(*kPtr);
+ *s = '1';
+ return s+1;
+ }
+ }
+ ++*s;
+ ++s;
+ return s;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * AdjustRange --
+ *
+ * Rescales a 'double' in preparation for formatting it using the
+ * 'quick' double-to-string method.
+ *
+ * Results:
+ * Returns the precision that has been lost in the prescaling as
+ * a count of units in the least significant place.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+AdjustRange(double* dPtr, /* INOUT: Number to adjust */
+ int k) /* IN: floor(log10(d)) */
+{
+ int ieps; /* Number of roundoff errors that have
+ * accumulated */
+ double d = *dPtr; /* Number to adjust */
+ double ds;
+ int i, j, j1;
+
+ ieps = 2;
+
+ if (k > 0) {
+ /*
+ * The number must be reduced to bring it into range.
+ */
+ ds = tens[k & 0xf];
+ j = k >> 4;
+ if (j & BLETCH) {
+ j &= (BLETCH-1);
+ d /= bigtens[N_BIGTENS - 1];
+ ieps++;
+ }
+ i = 0;
+ for (; j != 0; j>>=1) {
+ if (j & 1) {
+ ds *= bigtens[i];
+ ++ieps;
+ }
+ ++i;
+ }
+ d /= ds;
+ } else if ((j1 = -k) != 0) {
+ /*
+ * The number must be increased to bring it into range
+ */
+ d *= tens[j1 & 0xf];
+ i = 0;
+ for (j = j1>>4; j; j>>=1) {
+ if (j & 1) {
+ ieps++;
+ d *= bigtens[i];
+ }
+ ++i;
+ }
+ }
+
+ *dPtr = d;
+ return ieps;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningQuickFormat --
+ *
+ * Returns a 'quick' format of a double precision number to a string
+ * of digits, preferring a shorter string of digits if the shorter
+ * string is still within 1/2 ulp of the number.
+ *
+ * Results:
+ * Returns the string of digits. Returns NULL if the 'quick' method
+ * fails and the bignum method must be used.
+ *
+ * Side effects:
+ * Stores the position of the decimal point at '*kPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningQuickFormat(double d, /* Number to convert */
+ int k, /* floor(log10(d)) */
+ int ilim, /* Number of significant digits to return */
+ double eps,
+ /* Estimated roundoff error */
+ char* retval,
+ /* Buffer to receive the digit string */
+ int* kPtr)
+ /* Pointer to stash the position of
+ * the decimal point */
+{
+ char* s = retval; /* Cursor in the return value */
+ int digit; /* Current digit */
+ int i;
+
+ eps = 0.5 / tens[ilim-1] - eps;
+ i = 0;
+ for (;;) {
+ /* Convert a digit */
+
+ digit = d;
+ d -= digit;
+ *s++ = '0' + digit;
+
+ /*
+ * Truncate the conversion if the string of digits is within
+ * 1/2 ulp of the actual value.
+ */
+
+ if (d < eps) {
+ *kPtr = k;
+ return s;
+ }
+ if ((1. - d) < eps) {
+ *kPtr = k;
+ return BumpUp(s, retval, kPtr);
+ }
+
+ /*
+ * Bail out if the conversion fails to converge to a sufficiently
+ * precise value
+ */
+
+ if (++i >= ilim) {
+ return NULL;
+ }
+
+ /*
+ * Bring the next digit to the integer part.
+ */
+
+ eps *= 10;
+ d *= 10.0;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictQuickFormat --
+ *
+ * Convert a double precision number of a string of a precise number
+ * of digits, using the 'quick' double precision method.
+ *
+ * Results:
+ * Returns the digit string, or NULL if the bignum method must be
+ * used to do the formatting.
+ *
+ * Side effects:
+ * Stores the position of the decimal point in '*kPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictQuickFormat(double d, /* Number to convert */
+ int k, /* floor(log10(d)) */
+ int ilim, /* Number of significant digits to return */
+ double eps, /* Estimated roundoff error */
+ char* retval, /* Start of the digit string */
+ int* kPtr) /* Pointer to stash the position of
+ * the decimal point */
+{
+ char* s = retval; /* Cursor in the return value */
+ int digit; /* Current digit of the answer */
+ int i;
+
+ eps *= tens[ilim-1];
+ i = 1;
+ for (;;) {
+ /* Extract a digit */
+ digit = d;
+ d -= digit;
+ if (d == 0.0) {
+ ilim = i;
+ }
+ *s++ = '0' + digit;
+
+ /*
+ * When the given digit count is reached, handle trailing strings
+ * of 0 and 9.
+ */
+ if (i == ilim) {
+ if (d > 0.5 + eps) {
+ *kPtr = k;
+ return BumpUp(s, retval, kPtr);
+ } else if (d < 0.5 - eps) {
+ while (*--s == '0') {
+ /* do nothing */
+ }
+ s++;
+ *kPtr = k;
+ return s;
+ } else {
+ return NULL;
+ }
+ }
+
+ /* Advance to the next digit */
+ ++i;
+ d *= 10.0;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * QuickConversion --
+ *
+ * Converts a floating point number the 'quick' way, when only a limited
+ * number of digits is required and floating point arithmetic can
+ * therefore be used for the intermediate results.
+ *
+ * Results:
+ * Returns the converted string, or NULL if the bignum method must
+ * be used.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+QuickConversion(double d, /* Number to format */
+ int k, /* floor(log10(d)), approximately */
+ int k_check, /* 0 if k is exact, 1 if it may be too high */
+ int flags, /* Flags passed to dtoa:
+ * TCL_DD_SHORTEN_FLAG */
+ int len, /* Length of the return value */
+ int ilim, /* Number of digits to store */
+ int ilim1, /* Number of digits to store if we
+ * musguessed k */
+ int* decpt, /* OUTPUT: Location of the decimal point */
+ char** endPtr) /* OUTPUT: Pointer to the terminal null byte */
+{
+ int ieps; /* Number of 1-ulp roundoff errors that have
+ * accumulated in the calculation*/
+ Double eps; /* Estimated roundoff error */
+ char* retval; /* Returned string */
+ char* end; /* Pointer to the terminal null byte in the
+ * returned string */
+
/*
- * Find a large integer r, and integer e, such that
- * v = r * FLT_RADIX**e
- * and r is as small as possible. Also determine whether the significand
- * is the smallest possible.
+ * Bring d into the range [1 .. 10)
*/
+ ieps = AdjustRange(&d, k);
- smallestSig = GetIntegerTimesPower(v, &r, &e);
+ /*
+ * If the guessed value of k didn't get d into range, adjust it
+ * by one. If that leaves us outside the range in which quick format
+ * is accurate, bail out.
+ */
+ if (k_check && d < 1. && ilim > 0) {
+ if (ilim1 < 0) {
+ return NULL;
+ }
+ ilim = ilim1;
+ --k;
+ d *= 10.0;
+ ++ieps;
+ }
- lowOK = highOK = (mp_iseven(&r));
+ /*
+ * Compute estimated roundoff error
+ */
+ eps.d = ieps * d + 7.;
+ eps.w.word0 -= (FP_PRECISION-1) << EXP_SHIFT;
/*
- * We are going to want to develop integers r, s, mplus, and mminus such
- * that v = r / s, v(+)-v / 2 = mplus / s; v-v(-) / 2 = mminus / s and
- * then scale either s or r, mplus, mminus by an appropriate power of ten.
- *
- * We actually do this by keeping track of the powers of 2 and 5 by which
- * f is multiplied to yield v and by which 1 is multiplied to yield s,
- * mplus, and mminus.
+ * Handle the peculiar case where the result has no significant
+ * digits.
*/
+ retval = ckalloc(len + 1);
+ if (ilim == 0) {
+ d -= 5.;
+ if (d > eps.d) {
+ *retval = '1';
+ *decpt = k;
+ return retval;
+ } else if (d < -eps.d) {
+ *decpt = k;
+ return retval;
+ } else {
+ ckfree(retval);
+ return NULL;
+ }
+ }
- if (e >= 0) {
- int bits = e * log2FLT_RADIX;
+ /* Format the digit string */
- if (!smallestSig) {
- /*
- * Normal case, m+ and m- are both FLT_RADIX**e
- */
+ if (flags & TCL_DD_SHORTEN_FLAG) {
+ end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt);
+ } else {
+ end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt);
+ }
+ if (end == NULL) {
+ ckfree(retval);
+ return NULL;
+ }
+ *end = '\0';
+ if (endPtr != NULL) {
+ *endPtr = end;
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * CastOutPowersOf2 --
+ *
+ * Adjust the factors 'b2', 'm2', and 's2' to cast out common powers
+ * of 2 from numerator and denominator in preparation for the 'bignum'
+ * method of floating point conversion.
+ *
+ *-----------------------------------------------------------------------------
+ */
- rfac2 = bits + 1;
- sfac2 = 1;
- mplusfac2 = bits;
- mminusfac2 = bits;
+inline static void
+CastOutPowersOf2(int* b2, /* Power of 2 to multiply the significand */
+ int* m2, /* Power of 2 to multiply 1/2 ulp */
+ int* s2) /* Power of 2 to multiply the common
+ * denominator */
+{
+ int i;
+ if (*m2 > 0 && *s2 > 0) { /* Find the smallest power of 2 in the
+ * numerator */
+ if (*m2 < *s2) { /* Find the lowest common denominatorr */
+ i = *m2;
} else {
+ i = *s2;
+ }
+ *b2 -= i; /* Reduce to lowest terms */
+ *m2 -= i;
+ *s2 -= i;
+ }
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningInt64Conversion --
+ *
+ * Converts a double-precision number to the shortest string of
+ * digits that reconverts exactly to the given number, or to
+ * 'ilim' digits if that will yield a shorter result. The numerator and
+ * denominator in David Gay's conversion algorithm are known to fit
+ * in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
+ *
+ * Results:
+ * Returns the string of significant decimal digits, in newly
+ * allocated memory
+ *
+ * Side effects:
+ * Stores the location of the decimal point in '*decpt' and the
+ * location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningInt64Conversion(Double* dPtr,
+ /* Original number to convert */
+ int convType,
+ /* Type of conversion (shortest, Steele,
+ E format, F format) */
+ Tcl_WideUInt bw,
+ /* Integer significand */
+ int b2, int b5,
+ /* Scale factor for the significand
+ * in the numerator */
+ int m2plus, int m2minus, int m5,
+ /* Scale factors for 1/2 ulp in
+ * the numerator (will be different if
+ * bw == 1 */
+ int s2, int s5,
+ /* Scale factors for the denominator */
+ int k,
+ /* Number of output digits before the decimal
+ * point */
+ int len,
+ /* Number of digits to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Position of the terminal '\0'
+ * at the end of the returned string */
+{
+
+ char* retval = ckalloc(len + 1);
+ /* Output buffer */
+ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
+ /* Numerator of the fraction being converted */
+ Tcl_WideUInt S = wuipow5[s5] << s2;
+ /* Denominator of the fraction being
+ * converted */
+ Tcl_WideUInt mplus, mminus; /* Ranges for testing whether the result
+ * is within roundoff of being exact */
+ int digit; /* Current output digit */
+ char* s = retval; /* Cursor in the output buffer */
+ int i; /* Current position in the output buffer */
+
+ /* Adjust if the logarithm was guessed wrong */
+
+ if (b < S) {
+ b = 10 * b;
+ ++m2plus; ++m2minus; ++m5;
+ ilim = ilim1;
+ --k;
+ }
+
+ /* Compute roundoff ranges */
+
+ mplus = wuipow5[m5] << m2plus;
+ mminus = wuipow5[m5] << m2minus;
+
+ /* Loop through the digits */
+
+ i = 1;
+ for (;;) {
+ digit = (int)(b / S);
+ if (digit > 10) {
+ Tcl_Panic("wrong digit!");
+ }
+ b = b % S;
+
+ /*
+ * Does the current digit put us on the low side of the exact value
+ * but within within roundoff of being exact?
+ */
+ if (b < mplus
+ || (b == mplus
+ && convType != TCL_DD_STEELE0
+ && (dPtr->w.word1 & 1) == 0)) {
/*
- * If f is equal to the smallest significand, then we need another
- * factor of FLT_RADIX in s to cope with stepping to the next
- * smaller exponent when going to e's predecessor.
+ * Make sure we shouldn't be rounding *up* instead,
+ * in case the next number above is closer
*/
+ if (2 * b > S
+ || (2 * b == S
+ && (digit & 1) != 0)) {
+ ++digit;
+ if (digit == 10) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
+ break;
+ }
+ }
- rfac2 = bits + log2FLT_RADIX + 1;
- sfac2 = 1 + log2FLT_RADIX;
- mplusfac2 = bits + log2FLT_RADIX;
- mminusfac2 = bits;
+ /* Stash the current digit */
+
+ *s++ = '0' + digit;
+ break;
}
- } else {
+
/*
- * v has digits after the binary point
+ * Does one plus the current digit put us within roundoff of the
+ * number?
*/
+ if (b > S - mminus
+ || (b == S - mminus
+ && convType != TCL_DD_STEELE0
+ && (dPtr->w.word1 & 1) == 0)) {
+ if (digit == 9) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
+ break;
+ }
+ ++digit;
+ *s++ = '0' + digit;
+ break;
+ }
- if (e <= DBL_MIN_EXP-DBL_MANT_DIG || !smallestSig) {
- /*
- * Either f isn't the smallest significand or e is the smallest
- * exponent. mplus and mminus will both be 1.
- */
+ /*
+ * Have we converted all the requested digits?
+ */
+ *s++ = '0' + digit;
+ if (i == ilim) {
+ if (2*b > S
+ || (2*b == S && (digit & 1) != 0)) {
+ s = BumpUp(s, retval, &k);
+ }
+ break;
+ }
+
+ /* Advance to the next digit */
+
+ b = 10 * b;
+ mplus = 10 * mplus;
+ mminus = 10 * mminus;
+ ++i;
+ }
+
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
+ */
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictInt64Conversion --
+ *
+ * Converts a double-precision number to a fixed-length string of
+ * 'ilim' digits that reconverts exactly to the given number.
+ * ('ilim' should be replaced with 'ilim1' in the case where
+ * log10(d) has been overestimated). The numerator and
+ * denominator in David Gay's conversion algorithm are known to fit
+ * in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
+ *
+ * Results:
+ * Returns the string of significant decimal digits, in newly
+ * allocated memory
+ *
+ * Side effects:
+ * Stores the location of the decimal point in '*decpt' and the
+ * location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictInt64Conversion(Double* dPtr,
+ /* Original number to convert */
+ int convType,
+ /* Type of conversion (shortest, Steele,
+ E format, F format) */
+ Tcl_WideUInt bw,
+ /* Integer significand */
+ int b2, int b5,
+ /* Scale factor for the significand
+ * in the numerator */
+ int s2, int s5,
+ /* Scale factors for the denominator */
+ int k,
+ /* Number of output digits before the decimal
+ * point */
+ int len,
+ /* Number of digits to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Position of the terminal '\0'
+ * at the end of the returned string */
+{
+
+ char* retval = ckalloc(len + 1);
+ /* Output buffer */
+ Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
+ /* Numerator of the fraction being converted */
+ Tcl_WideUInt S = wuipow5[s5] << s2;
+ /* Denominator of the fraction being
+ * converted */
+ int digit; /* Current output digit */
+ char* s = retval; /* Cursor in the output buffer */
+ int i; /* Current position in the output buffer */
+
+ /* Adjust if the logarithm was guessed wrong */
+
+ if (b < S) {
+ b = 10 * b;
+ ilim = ilim1;
+ --k;
+ }
+
+ /* Loop through the digits */
+
+ i = 1;
+ for (;;) {
+ digit = (int)(b / S);
+ if (digit > 10) {
+ Tcl_Panic("wrong digit!");
+ }
+ b = b % S;
+
+ /*
+ * Have we converted all the requested digits?
+ */
+ *s++ = '0' + digit;
+ if (i == ilim) {
+ if (2*b > S
+ || (2*b == S && (digit & 1) != 0)) {
+ s = BumpUp(s, retval, &k);
+ }
+ break;
+ }
+
+ /* Advance to the next digit */
+
+ b = 10 * b;
+ ++i;
+ }
+
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
+ */
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpPowD --
+ *
+ * Test whether bankers' rounding should round a digit up. Assumption
+ * is made that the denominator of the fraction being tested is
+ * a power of 2**DIGIT_BIT.
+ *
+ * Results:
+ * Returns 1 iff the fraction is more than 1/2, or if the fraction
+ * is exactly 1/2 and the digit is odd.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpPowD(mp_int* b,
+ /* Numerator of the fraction */
+ int sd, /* Denominator is 2**(sd*DIGIT_BIT) */
+ int isodd)
+ /* 1 if the digit is odd, 0 if even */
+{
+ int i;
+ const static mp_digit topbit = (1<<(DIGIT_BIT-1));
+ if (b->used < sd || (b->dp[sd-1] & topbit) == 0) {
+ return 0;
+ }
+ if (b->dp[sd-1] != topbit) {
+ return 1;
+ }
+ for (i = sd-2; i >= 0; --i) {
+ if (b->dp[i] != 0) {
+ return 1;
+ }
+ }
+ return isodd;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpToNextPowD --
+ *
+ * Tests whether bankers' rounding will round down in the
+ * "denominator is a power of 2**MP_DIGIT" case.
+ *
+ * Results:
+ * Returns 1 if the rounding will be performed - which increases the
+ * digit by one - and 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpToNextPowD(mp_int* b,
+ /* Numerator of the fraction */
+ mp_int* m,
+ /* Numerator of the rounding tolerance */
+ int sd,
+ /* Common denominator is 2**(sd*DIGIT_BIT) */
+ int convType,
+ /* Conversion type: STEELE defeats
+ * round-to-even (Not sure why one wants to
+ * do this; I copied it from Gay) FIXME */
+ int isodd,
+ /* 1 if the integer significand is odd */
+ mp_int* temp)
+ /* Work area for the calculation */
+{
+ int i;
+
+ /*
+ * Compare B and S-m -- which is the same as comparing B+m and S --
+ * which we do by computing b+m and doing a bitwhack compare against
+ * 2**(DIGIT_BIT*sd)
+ */
+ mp_add(b, m, temp);
+ if (temp->used <= sd) { /* too few digits to be > S */
+ return 0;
+ }
+ if (temp->used > sd+1 || temp->dp[sd] > 1) {
+ /* >= 2s */
+ return 1;
+ }
+ for (i = sd-1; i >= 0; --i) {
+ /* check for ==s */
+ if (temp->dp[i] != 0) { /* > s */
+ return 1;
+ }
+ }
+ if (convType == TCL_DD_STEELE0) {
+ /* biased rounding */
+ return 0;
+ }
+ return isodd;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningBignumConversionPowD --
+ *
+ * Converts a double-precision number to the shortest string of
+ * digits that reconverts exactly to the given number, or to
+ * 'ilim' digits if that will yield a shorter result. The denominator
+ * in David Gay's conversion algorithm is known to be a power of
+ * 2**DIGIT_BIT, and hence the division in the main loop may be replaced
+ * by a digit shift and mask.
+ *
+ * Results:
+ * Returns the string of significant decimal digits, in newly
+ * allocated memory
+ *
+ * Side effects:
+ * Stores the location of the decimal point in '*decpt' and the
+ * location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningBignumConversionPowD(Double* dPtr,
+ /* Original number to convert */
+ int convType,
+ /* Type of conversion (shortest, Steele,
+ E format, F format) */
+ Tcl_WideUInt bw,
+ /* Integer significand */
+ int b2, int b5,
+ /* Scale factor for the significand
+ * in the numerator */
+ int m2plus, int m2minus, int m5,
+ /* Scale factors for 1/2 ulp in
+ * the numerator (will be different if
+ * bw == 1 */
+ int sd,
+ /* Scale factor for the denominator */
+ int k,
+ /* Number of output digits before the decimal
+ * point */
+ int len,
+ /* Number of digits to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Position of the terminal '\0'
+ * at the end of the returned string */
+{
+
+ char* retval = ckalloc(len + 1);
+ /* Output buffer */
+ mp_int b; /* Numerator of the fraction being converted */
+ mp_int mplus, mminus; /* Bounds for roundoff */
+ mp_digit digit; /* Current output digit */
+ char* s = retval; /* Cursor in the output buffer */
+ int i; /* Index in the output buffer */
+ mp_int temp;
+ int r1;
+
+ /*
+ * b = bw * 2**b2 * 5**b5
+ * mminus = 5**m5
+ */
+
+ TclBNInitBignumFromWideUInt(&b, bw);
+ mp_init_set_int(&mminus, 1);
+ MulPow5(&b, b5, &b);
+ mp_mul_2d(&b, b2, &b);
+
+ /* Adjust if the logarithm was guessed wrong */
+
+ if (b.used <= sd) {
+ mp_mul_d(&b, 10, &b);
+ ++m2plus; ++m2minus; ++m5;
+ ilim = ilim1;
+ --k;
+ }
+
+ /*
+ * mminus = 5**m5 * 2**m2minus
+ * mplus = 5**m5 * 2**m2plus
+ */
+
+ mp_mul_2d(&mminus, m2minus, &mminus);
+ MulPow5(&mminus, m5, &mminus);
+ if (m2plus > m2minus) {
+ mp_init_copy(&mplus, &mminus);
+ mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
+ }
+ mp_init(&temp);
+
+ /* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+ * by mp_digit extraction */
- rfac2 = 1;
- sfac2 = 1 - e * log2FLT_RADIX;
- mplusfac2 = 0;
- mminusfac2 = 0;
+ i = 0;
+ for (;;) {
+ if (b.used <= sd) {
+ digit = 0;
} else {
+ digit = b.dp[sd];
+ if (b.used > sd+1 || digit >= 10) {
+ Tcl_Panic("wrong digit!");
+ }
+ --b.used; mp_clamp(&b);
+ }
+
+ /*
+ * Does the current digit put us on the low side of the exact value
+ * but within within roundoff of being exact?
+ */
+
+ r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
+ if (r1 == MP_LT
+ || (r1 == MP_EQ
+ && convType != TCL_DD_STEELE0
+ && (dPtr->w.word1 & 1) == 0)) {
/*
- * f is the smallest significand, but e is not the smallest
- * exponent. We need to scale by FLT_RADIX again to cope with the
- * fact that v's predecessor has a smaller exponent.
+ * Make sure we shouldn't be rounding *up* instead,
+ * in case the next number above is closer
*/
+ if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+ ++digit;
+ if (digit == 10) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
+ break;
+ }
+ }
+
+ /* Stash the last digit */
+
+ *s++ = '0' + digit;
+ break;
+ }
+
+ /*
+ * Does one plus the current digit put us within roundoff of the
+ * number?
+ */
+
+ if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd,
+ convType, dPtr->w.word1 & 1,
+ &temp)) {
+ if (digit == 9) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
+ break;
+ }
+ ++digit;
+ *s++ = '0' + digit;
+ break;
+ }
- rfac2 = 1 + log2FLT_RADIX;
- sfac2 = 1 + log2FLT_RADIX * (1 - e);
- mplusfac2 = FLT_RADIX;
- mminusfac2 = 0;
+ /*
+ * Have we converted all the requested digits?
+ */
+ *s++ = '0' + digit;
+ if (i == ilim) {
+ if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+ s = BumpUp(s, retval, &k);
+ }
+ break;
+ }
+
+ /* Advance to the next digit */
+
+ mp_mul_d(&b, 10, &b);
+ mp_mul_d(&mminus, 10, &mminus);
+ if (m2plus > m2minus) {
+ mp_mul_2d(&mminus, m2plus-m2minus, &mplus);
}
+ ++i;
}
- /*
- * Estimate the highest power of ten that will be needed to hold the
- * result.
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
+ */
+ if (m2plus > m2minus) {
+ mp_clear(&mplus);
+ }
+ mp_clear_multi(&b, &mminus, &temp, NULL);
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * StrictBignumConversionPowD --
+ *
+ * Converts a double-precision number to a fixed-lengt string of
+ * 'ilim' digits (or 'ilim1' if log10(d) has been overestimated.)
+ * The denominator in David Gay's conversion algorithm is known to
+ * be a power of 2**DIGIT_BIT, and hence the division in the main
+ * loop may be replaced by a digit shift and mask.
+ *
+ * Results:
+ * Returns the string of significant decimal digits, in newly
+ * allocated memory.
+ *
+ * Side effects:
+ * Stores the location of the decimal point in '*decpt' and the
+ * location of the terminal null byte in '*endPtr'.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+StrictBignumConversionPowD(Double* dPtr,
+ /* Original number to convert */
+ int convType,
+ /* Type of conversion (shortest, Steele,
+ E format, F format) */
+ Tcl_WideUInt bw,
+ /* Integer significand */
+ int b2, int b5,
+ /* Scale factor for the significand
+ * in the numerator */
+ int sd,
+ /* Scale factor for the denominator */
+ int k,
+ /* Number of output digits before the decimal
+ * point */
+ int len,
+ /* Number of digits to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Position of the terminal '\0'
+ * at the end of the returned string */
+{
+
+ char* retval = ckalloc(len + 1);
+ /* Output buffer */
+ mp_int b; /* Numerator of the fraction being converted */
+ mp_digit digit; /* Current output digit */
+ char* s = retval; /* Cursor in the output buffer */
+ int i; /* Index in the output buffer */
+ mp_int temp;
+
+ /*
+ * b = bw * 2**b2 * 5**b5
*/
- k = (int) ceil(log(v) / log(10.));
- if (k >= 0) {
- sfac2 += k;
- sfac5 = k;
- } else {
- rfac2 -= k;
- mplusfac2 -= k;
- mminusfac2 -= k;
- rfac5 = -k;
+ TclBNInitBignumFromWideUInt(&b, bw);
+ MulPow5(&b, b5, &b);
+ mp_mul_2d(&b, b2, &b);
+
+ /* Adjust if the logarithm was guessed wrong */
+
+ if (b.used <= sd) {
+ mp_mul_d(&b, 10, &b);
+ ilim = ilim1;
+ --k;
}
+ mp_init(&temp);
- /*
- * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5.
+ /*
+ * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
+ * by mp_digit extraction
*/
- mp_init_set(&mplus, 1);
- for (i=0 ; i<=8 ; ++i) {
- if (rfac5 & (1 << i)) {
- mp_mul(&mplus, pow5+i, &mplus);
+ i = 1;
+ for (;;) {
+ if (b.used <= sd) {
+ digit = 0;
+ } else {
+ digit = b.dp[sd];
+ if (b.used > sd+1 || digit >= 10) {
+ Tcl_Panic("wrong digit!");
+ }
+ --b.used; mp_clamp(&b);
}
- }
- mp_mul(&r, &mplus, &r);
- mp_mul_2d(&r, rfac2, &r);
- mp_init_copy(&mminus, &mplus);
- mp_mul_2d(&mplus, mplusfac2, &mplus);
- mp_mul_2d(&mminus, mminusfac2, &mminus);
- mp_init_set(&s, 1);
- for (i=0 ; i<=8 ; ++i) {
- if (sfac5 & (1 << i)) {
- mp_mul(&s, pow5+i, &s);
+
+ /*
+ * Have we converted all the requested digits?
+ */
+ *s++ = '0' + digit;
+ if (i == ilim) {
+ if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
+ s = BumpUp(s, retval, &k);
+ }
+ break;
}
+
+ /* Advance to the next digit */
+
+ mp_mul_d(&b, 10, &b);
+ ++i;
}
- mp_mul_2d(&s, sfac2, &s);
- /*
- * It is possible for k to be off by one because we used an inexact
- * logarithm.
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
*/
+ mp_clear_multi(&b, &temp, NULL);
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUp --
+ *
+ * Tests whether a digit should be rounded up or down when finishing
+ * bignum-based floating point conversion.
+ *
+ * Results:
+ * Returns 1 if the number needs to be rounded up, 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
- mp_init(&temp);
- mp_add(&r, &mplus, &temp);
- i = mp_cmp_mag(&temp, &s);
- if (i>0 || (highOK && i==0)) {
- mp_mul_d(&s, 10, &s);
- k++;
- } else {
- mp_mul_d(&temp, 10, &temp);
- i = mp_cmp_mag(&temp, &s);
- if (i<0 || (highOK && i==0)) {
- mp_mul_d(&r, 10, &r);
- mp_mul_d(&mplus, 10, &mplus);
- mp_mul_d(&mminus, 10, &mminus);
- k--;
+inline static int
+ShouldBankerRoundUp(mp_int* twor,
+ /* 2x the remainder from thd division that
+ * produced the last digit */
+ mp_int* S, /* Denominator */
+ int isodd) /* Flag == 1 if the last digit is odd */
+{
+ int r = mp_cmp_mag(twor, S);
+ switch (r) {
+ case MP_LT:
+ return 0;
+ case MP_EQ:
+ return isodd;
+ case MP_GT:
+ return 1;
+ }
+ Tcl_Panic("in ShouldBankerRoundUp, trichotomy fails!");
+ return 0;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShouldBankerRoundUpToNext --
+ *
+ * Tests whether the remainder is great enough to force rounding
+ * to the next higher digit.
+ *
+ * Results:
+ * Returns 1 if the number should be rounded up, 0 otherwise.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static int
+ShouldBankerRoundUpToNext(mp_int* b,
+ /* Remainder from the division that produced
+ * the last digit. */
+ mp_int* m,
+ /* Numerator of the rounding tolerance */
+ mp_int* S,
+ /* Denominator */
+ int convType,
+ /* Conversion type: STEELE0 defeats
+ * round-to-even. (Not sure why one would
+ * want this; I coped it from Gay. FIXME */
+ int isodd,
+ /* 1 if the integer significand is odd */
+ mp_int* temp)
+ /* Work area needed for the calculation */
+{
+ int r;
+ /* Compare b and S-m: this is the same as comparing B+m and S. */
+ mp_add(b, m, temp);
+ r = mp_cmp_mag(temp, S);
+ switch(r) {
+ case MP_LT:
+ return 0;
+ case MP_EQ:
+ if (convType == TCL_DD_STEELE0) {
+ return 0;
+ } else {
+ return isodd;
}
+ case MP_GT:
+ return 1;
}
+ Tcl_Panic("in ShouldBankerRoundUpToNext, trichotomy fails!");
+ return 0;
+}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * ShorteningBignumConversion --
+ *
+ * Convert a floating point number to a variable-length digit string
+ * using the multiprecision method.
+ *
+ * Results:
+ * Returns the string of digits.
+ *
+ * Side effects:
+ * Stores the position of the decimal point in *decpt.
+ * Stores a pointer to the end of the number in *endPtr.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+inline static char*
+ShorteningBignumConversion(Double* dPtr,
+ /* Original number being converted */
+ int convType,
+ /* Conversion type */
+ Tcl_WideUInt bw,
+ /* Integer significand and exponent */
+ int b2,
+ /* Scale factor for the significand */
+ int m2plus, int m2minus,
+ /* Scale factors for 1/2 ulp in numerator */
+ int s2, int s5,
+ /* Scale factors for denominator */
+ int k,
+ /* Guessed position of the decimal point */
+ int len,
+ /* Size of the digit buffer to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Pointer to the end of the number */
+{
+ char* retval = ckalloc(len+1);
+ /* Buffer of digits to return */
+ char* s = retval; /* Cursor in the return value */
+ mp_int b; /* Numerator of the result */
+ mp_int mminus; /* 1/2 ulp below the result */
+ mp_int mplus; /* 1/2 ulp above the result */
+ mp_int S; /* Denominator of the result */
+ mp_int dig; /* Current digit of the result */
+ int digit; /* Current digit of the result */
+ mp_int temp; /* Work area */
+ int minit = 1; /* Fudge factor for when we misguess k */
+ int i;
+ int r1;
/*
- * At this point, k contains the power of ten by which we're scaling the
- * result. r/s is at least 1/10 and strictly less than ten, and v = r/s *
- * 10**k. mplus and mminus give the rounding limits.
+ * b = bw * 2**b2 * 5**b5
+ * S = 2**s2 * 5*s5
*/
- for (;;) {
- int tc1, tc2;
+ 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);
- mp_mul_d(&r, 10, &r);
- mp_div(&r, &s, &temp, &r); /* temp = 10r / s; r = 10r mod s */
- i = temp.dp[0];
- mp_mul_d(&mplus, 10, &mplus);
- mp_mul_d(&mminus, 10, &mminus);
- tc1 = mp_cmp_mag(&r, &mminus);
- if (lowOK) {
- tc1 = (tc1 <= 0);
- } else {
- tc1 = (tc1 < 0);
+ /*
+ * Handle the case where we guess the position of the decimal point
+ * wrong.
+ */
+
+ if (mp_cmp_mag(&b, &S) == MP_LT) {
+ mp_mul_d(&b, 10, &b);
+ minit = 10;
+ ilim =ilim1;
+ --k;
+ }
+
+ /* mminus = 2**m2minus * 5**m5 */
+
+ mp_init_set_int(&mminus, minit);
+ mp_mul_2d(&mminus, m2minus, &mminus);
+ if (m2plus > m2minus) {
+ mp_init_copy(&mplus, &mminus);
+ mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
+ }
+ mp_init(&temp);
+
+ /* Loop through the digits */
+
+ mp_init(&dig);
+ i = 1;
+ for (;;) {
+ mp_div(&b, &S, &dig, &b);
+ if (dig.used > 1 || dig.dp[0] >= 10) {
+ Tcl_Panic("wrong digit!");
}
- mp_add(&r, &mplus, &temp);
- tc2 = mp_cmp_mag(&temp, &s);
- if (highOK) {
- tc2 = (tc2 >= 0);
- } else {
- tc2 = (tc2 > 0);
+ digit = dig.dp[0];
+
+ /*
+ * Does the current digit leave us with a remainder small enough to
+ * round to it?
+ */
+
+ r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
+ if (r1 == MP_LT
+ || (r1 == MP_EQ
+ && convType != TCL_DD_STEELE0
+ && (dPtr->w.word1 & 1) == 0)) {
+ mp_mul_2d(&b, 1, &b);
+ if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+ ++digit;
+ if (digit == 10) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
+ break;
+ }
+ }
+ *s++ = '0' + digit;
+ break;
}
- if (!tc1) {
- if (!tc2) {
- *buffer++ = '0' + i;
- } else {
- c = (char) (i + '1');
+
+ /*
+ * Does the current digit leave us with a remainder large enough
+ * to commit to rounding up to the next higher digit?
+ */
+
+ if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType,
+ dPtr->w.word1 & 1, &temp)) {
+ ++digit;
+ if (digit == 10) {
+ *s++ = '9';
+ s = BumpUp(s, retval, &k);
break;
}
- } else {
- if (!tc2) {
- c = (char) (i + '0');
- } else {
- mp_mul_2d(&r, 1, &r);
- n = mp_cmp_mag(&r, &s);
- if (n < 0) {
- c = (char) (i + '0');
- } else {
- c = (char) (i + '1');
- }
+ *s++ = '0' + digit;
+ break;
+ }
+
+ /* Have we converted all the requested digits? */
+
+ *s++ = '0' + digit;
+ if (i == ilim) {
+ mp_mul_2d(&b, 1, &b);
+ if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+ s = BumpUp(s, retval, &k);
}
break;
}
- };
- *buffer++ = c;
- *buffer++ = '\0';
- /*
- * Free memory, and return.
+ /* Advance to the next digit */
+
+ if (s5 > 0) {
+
+ /* Can possibly shorten the denominator */
+ mp_mul_2d(&b, 1, &b);
+ mp_mul_2d(&mminus, 1, &mminus);
+ if (m2plus > m2minus) {
+ mp_mul_2d(&mplus, 1, &mplus);
+ }
+ mp_div_d(&S, 5, &S, NULL);
+ --s5;
+ /*
+ * TODO: It might possibly be a win to fall back to
+ * int64 arithmetic here if S < 2**64/10. But it's
+ * a win only for a fairly narrow range of magnitudes
+ * so perhaps not worth bothering. We already know that
+ * we shorten the denominator by at least 1 mp_digit, perhaps
+ * 2. as we do the conversion for 17 digits of significance.
+ * Possible savings:
+ * 10**26 1 trip through loop before fallback possible
+ * 10**27 1 trip
+ * 10**28 2 trips
+ * 10**29 3 trips
+ * 10**30 4 trips
+ * 10**31 5 trips
+ * 10**32 6 trips
+ * 10**33 7 trips
+ * 10**34 8 trips
+ * 10**35 9 trips
+ * 10**36 10 trips
+ * 10**37 11 trips
+ * 10**38 12 trips
+ * 10**39 13 trips
+ * 10**40 14 trips
+ * 10**41 15 trips
+ * 10**42 16 trips
+ * thereafter no gain.
+ */
+ } else {
+ mp_mul_d(&b, 10, &b);
+ mp_mul_d(&mminus, 10, &mminus);
+ if (m2plus > m2minus) {
+ mp_mul_2d(&mplus, 10, &mplus);
+ }
+ }
+
+ ++i;
+ }
+
+
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
*/
+ if (m2plus > m2minus) {
+ mp_clear(&mplus);
+ }
+ mp_clear_multi(&b, &mminus, &temp, NULL);
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
- mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL);
- return k;
}
-
+
/*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
*
- * AbsoluteValue --
+ * StrictBignumConversion --
*
- * Splits a 'double' into its absolute value and sign.
+ * Convert a floating point number to a fixed-length digit string
+ * using the multiprecision method.
*
* Results:
- * Returns the absolute value.
+ * Returns the string of digits.
*
* Side effects:
- * Stores the signum in '*signum'.
+ * Stores the position of the decimal point in *decpt.
+ * Stores a pointer to the end of the number in *endPtr.
*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
*/
-static double
-AbsoluteValue(
- double v, /* Number to split */
- int *signum) /* (Output) Sign of the number 1=-, 0=+ */
+inline static char*
+StrictBignumConversion(Double* dPtr,
+ /* Original number being converted */
+ int convType,
+ /* Conversion type */
+ Tcl_WideUInt bw,
+ /* Integer significand and exponent */
+ int b2, /* Scale factor for the significand */
+ int s2, int s5,
+ /* Scale factors for denominator */
+ int k, /* Guessed position of the decimal point */
+ int len, /* Size of the digit buffer to allocate */
+ int ilim,
+ /* Number of digits to convert if b >= s */
+ int ilim1,
+ /* Number of digits to convert if b < s */
+ int* decpt,
+ /* OUTPUT: Position of the decimal point */
+ char** endPtr)
+ /* OUTPUT: Pointer to the end of the number */
{
+ char* retval = ckalloc(len+1);
+ /* Buffer of digits to return */
+ char* s = retval; /* Cursor in the return value */
+ mp_int b; /* Numerator of the result */
+ mp_int S; /* Denominator of the result */
+ mp_int dig; /* Current digit of the result */
+ int digit; /* Current digit of the result */
+ mp_int temp; /* Work area */
+ int g; /* Size of the current digit groun */
+ int i, j;
+
/*
- * Take the absolute value of the number, and report the number's sign.
- * Take special steps to preserve signed zeroes in IEEE floating point.
- * (We can't use fpclassify, because that's a C9x feature and we still
- * have to build on C89 compilers.)
+ * b = bw * 2**b2 * 5**b5
+ * S = 2**s2 * 5*s5
*/
-#ifndef IEEE_FLOATING_POINT
- if (v >= 0.0) {
- *signum = 0;
- } else {
- *signum = 1;
- v = -v;
- }
-#else
- union {
- Tcl_WideUInt iv;
- double dv;
- } bitwhack;
- bitwhack.dv = v;
- if (n770_fp) {
- bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
+ TclBNInitBignumFromWideUInt(&b, bw);
+ mp_mul_2d(&b, b2, &b);
+ mp_init_set_int(&S, 1);
+ MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
+
+ /*
+ * Handle the case where we guess the position of the decimal point
+ * wrong.
+ */
+
+ if (mp_cmp_mag(&b, &S) == MP_LT) {
+ mp_mul_d(&b, 10, &b);
+ ilim =ilim1;
+ --k;
}
- if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
- *signum = 1;
- bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
- if (n770_fp) {
- bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
+ mp_init(&temp);
+
+ /* Convert the leading digit */
+
+ mp_init(&dig);
+ i = 0;
+ mp_div(&b, &S, &dig, &b);
+ if (dig.used > 1 || dig.dp[0] >= 10) {
+ Tcl_Panic("wrong digit!");
+ }
+ digit = dig.dp[0];
+
+ /* Is a single digit all that was requested? */
+
+ *s++ = '0' + digit;
+ if (++i >= ilim) {
+ mp_mul_2d(&b, 1, &b);
+ if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+ s = BumpUp(s, retval, &k);
}
- v = bitwhack.dv;
} else {
- *signum = 0;
+
+ for (;;) {
+
+ /* Shift by a group of digits. */
+
+ g = ilim - i;
+ if (g > DIGIT_GROUP) {
+ g = DIGIT_GROUP;
+ }
+ if (s5 >= g) {
+ mp_div_d(&S, dpow5[g], &S, NULL);
+ s5 -= g;
+ } else if (s5 > 0) {
+ mp_div_d(&S, dpow5[s5], &S, NULL);
+ mp_mul_d(&b, dpow5[g - s5], &b);
+ s5 = 0;
+ } else {
+ mp_mul_d(&b, dpow5[g], &b);
+ }
+ mp_mul_2d(&b, g, &b);
+
+ /*
+ * As with the shortening bignum conversion, it's possible at
+ * this point that we will have reduced the denominator to
+ * less than 2**64/10, at which point it would be possible to
+ * fall back to to int64 arithmetic. But the potential payoff
+ * is tremendously less - unless we're working in F format -
+ * because we know that three groups of digits will always
+ * suffice for %#.17e, the longest format that doesn't introduce
+ * empty precision.
+ */
+
+ /* Extract the next digit */
+
+ mp_div(&b, &S, &dig, &b);
+ if (dig.used > 1) {
+ Tcl_Panic("wrong digit!");
+ }
+ digit = dig.dp[0];
+ for (j = g-1; j >= 0; --j) {
+ int t = itens[j];
+ *s++ = digit / t + '0';
+ digit %= t;
+ }
+ i += g;
+
+ /* Have we converted all the requested digits? */
+
+ if (i == ilim) {
+ mp_mul_2d(&b, 1, &b);
+ if (ShouldBankerRoundUp(&b, &S, digit&1)) {
+ s = BumpUp(s, retval, &k);
+ }
+ break;
+ }
+ }
}
-#endif
- return v;
+ /*
+ * Endgame - store the location of the decimal point and the end of the
+ * string.
+ */
+ mp_clear_multi(&b, &temp, NULL);
+ *s = '\0';
+ *decpt = k;
+ if (endPtr) {
+ *endPtr = s;
+ }
+ return retval;
+
}
/*
- *----------------------------------------------------------------------
+ *-----------------------------------------------------------------------------
*
- * GetIntegerTimesPower --
+ * TclDoubleDigits --
*
- * Converts a floating point number to an exact integer times a power of
- * the floating point radix.
+ * Core of Tcl's conversion of double-precision floating point numbers
+ * to decimal.
*
* Results:
- * Returns 1 if it converted the smallest significand, 0 otherwise.
+ * Returns a newly-allocated string of digits.
*
* Side effects:
- * Initializes the integer value (does not just assign it), and stores
- * the exponent.
+ * Sets *decpt to the index of the character in the string before the
+ * place that the decimal point should go. If 'endPtr' is not NULL,
+ * sets endPtr to point to the terminating '\0' byte of the string.
+ * Sets *sign to 1 if a minus sign should be printed with the number,
+ * or 0 if a plus sign (or no sign) should appear.
*
- *----------------------------------------------------------------------
+ * This function is a service routine that produces the string of digits
+ * for floating-point-to-decimal conversion. It can do a number of things
+ * according to the 'flags' argument. Valid values for 'flags' include:
+ * TCL_DD_SHORTEST - This is the default for floating point conversion
+ * if ::tcl_precision is 0. It constructs the shortest string
+ * of digits that will reconvert to the given number when scanned.
+ * For floating point numbers that are exactly between two
+ * decimal numbers, it resolves using the 'round to even' rule.
+ * With this value, the 'ndigits' parameter is ignored.
+ * TCL_DD_STEELE - This value is not recommended and may be removed
+ * in the future. It follows the conversion algorithm outlined
+ * in "How to Print Floating-Point Numbers Accurately" by
+ * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90,
+ * pp. 112-126]. This rule has the effect of rendering 1e23
+ * as 9.9999999999999999e22 - which is a 'better' approximation
+ * in the sense that it will reconvert correctly even if
+ * a subsequent input conversion is 'round up' or 'round down'
+ * rather than 'round to nearest', but is surprising otherwise.
+ * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e
+ * format conversion (or for default floating->string if
+ * tcl_precision is not 0). It constructs a string of at most
+ * 'ndigits' digits, choosing the one that is closest to the
+ * given number (and resolving ties with 'round to even').
+ * It is allowed to return fewer than 'ndigits' if the number
+ * converts exactly; if the TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG
+ * is supplied instead, it is also allowed to return fewer digits
+ * if the shorter string will still reconvert to the given
+ * input number.
+ * TCL_DD_F_FORMAT - This value is used to prepare numbers for %f
+ * format conversion. It requests that conversion proceed until
+ * 'ndigits' digits after the decimal point have been converted.
+ * It is possible for this format to result in a zero-length
+ * string if the number is sufficiently small. Again, it
+ * is permissible for TCL_DD_F_FORMAT to return fewer digits
+ * for a number that converts exactly, and changing the
+ * argument to TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow
+ * the routine also to return fewer digits if the shorter string
+ * will still reconvert without loss to the given input number.
+ *
+ * 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.
+ *
+ *-----------------------------------------------------------------------------
*/
-
-static int
-GetIntegerTimesPower(
- double v, /* Value to convert */
- mp_int *rPtr, /* (Output) Integer value */
- int *ePtr) /* (Output) Power of FLT_RADIX by which r must
- * be multiplied to yield v*/
+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 a, f;
- int e, i, n;
+ int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK);
+ /* Type of conversion being performed
+ * TCL_DD_SHORTEST0
+ * TCL_DD_STEELE0
+ * TCL_DD_E_FORMAT
+ * TCL_DD_F_FORMAT */
+ Double d; /* Union for deconstructing doubles */
+ Tcl_WideUInt bw; /* Integer significand */
+ int be; /* Power of 2 by which b must be multiplied */
+ int bbits; /* Number of bits needed to represent b */
+ int denorm; /* Flag == 1 iff the input number was
+ * denormalized */
+ int k; /* Estimate of floor(log10(d)) */
+ int k_check; /* Flag == 1 if d is near enough to a
+ * power of ten that k must be checked */
+ int b2, b5, s2, s5; /* Powers of 2 and 5 in the numerator and
+ * denominator of intermediate results */
+ int ilim, ilim1;
+ char* retval; /* Return value from this function */
+ int i;
- /*
- * Develop f and e such that v = f * FLT_RADIX**e, with
- * 1.0/FLT_RADIX <= f < 1.
+ /* Put the input number into a union for bit-whacking */
+
+ d.d = dv;
+
+ /*
+ * Handle the cases of negative numbers (by taking the absolute value:
+ * this includes -Inf and -NaN!), infinity, Not a Number, and zero.
*/
- f = frexp(v, &e);
-#if FLT_RADIX > 2
- n = e % log2FLT_RADIX;
- if (n > 0) {
- n -= log2FLT_RADIX;
- e += 1;
- f *= ldexp(1.0, n);
+ TakeAbsoluteValue(&d, sign);
+ if ((d.w.word0 & EXP_MASK) == EXP_MASK) {
+ return FormatInfAndNaN(&d, decpt, endPtr);
}
- e = (e - n) / log2FLT_RADIX;
-#endif
- if (f == 1.0) {
- f = 1.0 / FLT_RADIX;
- e += 1;
+ if (d.d == 0.0) {
+ return FormatZero(decpt, endPtr);
}
+ /*
+ * Unpack the floating point into a wide integer and an exponent.
+ * Determine the number of bits that the big integer requires, and
+ * compute a quick approximation (which may be one too high) of
+ * ceil(log10(d.d)).
+ */
+ denorm = ((d.w.word0 & EXP_MASK) == 0);
+ DoubleToExpAndSig(d.d, &bw, &be, &bbits);
+ k = ApproximateLog10(bw, be, bbits);
+ k = BetterLog10(d.d, k, &k_check);
+
+ /* At this point, we have:
+ * d is the number to convert.
+ * bw are significand and exponent: d == bw*2**be,
+ * bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits
+ * k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0
+ * if we know that k is exactly ceil(log10(d)) and 1 if we need to
+ * check.
+ * We want a rational number
+ * r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5),
+ * with b2, b5, s2, s5 >= 0. Note that the most significant decimal
+ * digit is floor(r) and that successive digits can be obtained
+ * by setting r <- 10*floor(r) (or b <= 10 * (b % S)).
+ * Find appropriate b2, b5, s2, s5.
+ */
+
+ ComputeScale(be, k, &b2, &b5, &s2, &s5);
+
/*
- * If the original number was denormalized, adjust e and f to be denormal
- * as well.
+ * Correct an incorrect caller-supplied 'ndigits'.
+ * Also determine:
+ * i = The maximum number of decimal digits that will be returned in the
+ * formatted string. This is k + 1 + ndigits for F format, 18 for
+ * shortest and Steele, and ndigits for E format.
+ * ilim = The number of significant digits to convert if
+ * k has been guessed correctly. This is -1 for shortest and Steele
+ * (which stop when all significance has been lost), 'ndigits'
+ * for E format, and 'k + 1 + ndigits' for F format.
+ * ilim1 = The minimum number of significant digits to convert if
+ * k has been guessed 1 too high. This, too, is -1 for shortest
+ * and Steele, and 'ndigits' for E format, but it's 'ndigits-1'
+ * for F format.
*/
- if (e < DBL_MIN_EXP) {
- n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX;
- f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX);
- e = DBL_MIN_EXP;
- n = (n + DIGIT_BIT - 1) / DIGIT_BIT;
- } else {
- n = mantDIGIT;
+ SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1);
+
+ /*
+ * Try to do low-precision conversion in floating point rather
+ * than resorting to expensive multiprecision arithmetic
+ */
+ if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) {
+ if ((retval = QuickConversion(d.d, k, k_check, flags,
+ i, ilim, ilim1,
+ decpt, endPtr)) != NULL) {
+ return retval;
+ }
}
- /*
- * Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
- * integer r. Preserve the invariant v = r * 2**rfac2 * FLT_RADIX**e by
- * adjusting e.
- */
-
- a = f;
- n = mantDIGIT;
- mp_init_size(rPtr, n);
- rPtr->used = n;
- rPtr->sign = MP_ZPOS;
- i = (mantBits % DIGIT_BIT);
- if (i == 0) {
- i = DIGIT_BIT;
- }
- while (n > 0) {
- a *= ldexp(1.0, i);
- i = DIGIT_BIT;
- rPtr->dp[--n] = (mp_digit) a;
- a -= (mp_digit) a;
- }
- *ePtr = e - DBL_MANT_DIG;
- return (f == 1.0 / FLT_RADIX);
+ /*
+ * For shortening conversions, determine the upper and lower bounds
+ * for the remainder at which we can stop.
+ * m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the
+ * high side, and
+ * m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the
+ * low side.
+ * We may need to increase s2 to put m2plus, m2minus, b2 over a
+ * common denominator.
+ */
+
+ if (flags & TCL_DD_SHORTEN_FLAG) {
+ int m2minus = b2;
+ int m2plus;
+ int m5 = b5;
+ int len = i;
+
+ /*
+ * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5)
+ * is 1/2 unit in the least significant place of the floating
+ * point number.
+ */
+ if (denorm) {
+ i = be + EXPONENT_BIAS + (FP_PRECISION-1);
+ } else {
+ i = 1 + FP_PRECISION - bbits;
+ }
+ b2 += i;
+ s2 += i;
+
+ /*
+ * Reduce the fractions to lowest terms, since the above calculation
+ * may have left excess powers of 2 in numerator and denominator
+ */
+ CastOutPowersOf2(&b2, &m2minus, &s2);
+
+ /*
+ * In the special case where bw==1, the nearest floating point number
+ * to it on the low side is 1/4 ulp below it. Adjust accordingly.
+ */
+ m2plus = m2minus;
+ if (!denorm && bw == 1) {
+ ++b2;
+ ++s2;
+ ++m2plus;
+ }
+
+ if (s5+1 < N_LOG2POW5
+ && s2+1 + log2pow5[s5+1] <= 64) {
+ /*
+ * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
+ * word, then all our intermediate calculations can be done
+ * using exact 64-bit arithmetic with no need for expensive
+ * multiprecision operations. (This will be true for all numbers
+ * in the range [1.0e-3 .. 1.0e+24]).
+ */
+
+ return ShorteningInt64Conversion(&d, convType, bw, b2, b5,
+ m2plus, m2minus, m5,
+ s2, s5, k, len, ilim, ilim1,
+ decpt, endPtr);
+ } else if (s5 == 0) {
+ /*
+ * The denominator is a power of 2, so we can replace division
+ * by digit shifts. First we round up s2 to a multiple of
+ * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+ * into a version of the comparison that's specialized for
+ * the 'power of mp_digit in the denominator' case.
+ */
+ if (s2 % DIGIT_BIT != 0) {
+ int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
+ b2 += delta;
+ m2plus += delta;
+ m2minus += delta;
+ s2 += delta;
+ }
+ return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5,
+ m2plus, m2minus, m5,
+ s2/DIGIT_BIT, k, len,
+ ilim, ilim1, decpt, endPtr);
+ } else {
+
+ /*
+ * Alas, there's no helpful special case; use full-up
+ * bignum arithmetic for the conversion
+ */
+
+ return ShorteningBignumConversion(&d, convType, bw,
+ b2, m2plus, m2minus,
+ s2, s5, k, len,
+ ilim, ilim1, decpt, endPtr);
+
+ }
+
+ } else {
+
+ /* Non-shortening conversion */
+
+ int len = i;
+
+ /* Reduce numerator and denominator to lowest terms */
+
+ if (b2 >= s2 && s2 > 0) {
+ b2 -= s2; s2 = 0;
+ } else if (s2 >= b2 && b2 > 0) {
+ s2 -= b2; b2 = 0;
+ }
+
+ if (s5+1 < N_LOG2POW5
+ && s2+1 + log2pow5[s5+1] <= 64) {
+ /*
+ * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
+ * word, then all our intermediate calculations can be done
+ * using exact 64-bit arithmetic with no need for expensive
+ * multiprecision operations.
+ */
+
+ return StrictInt64Conversion(&d, convType, bw, b2, b5,
+ s2, s5, k, len, ilim, ilim1,
+ decpt, endPtr);
+
+ } else if (s5 == 0) {
+ /*
+ * The denominator is a power of 2, so we can replace division
+ * by digit shifts. First we round up s2 to a multiple of
+ * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
+ * into a version of the comparison that's specialized for
+ * the 'power of mp_digit in the denominator' case.
+ */
+ if (s2 % DIGIT_BIT != 0) {
+ int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
+ b2 += delta;
+ s2 += delta;
+ }
+ return StrictBignumConversionPowD(&d, convType, bw, b2, b5,
+ s2/DIGIT_BIT, k, len,
+ ilim, ilim1, decpt, endPtr);
+ } else {
+ /*
+ * There are no helpful special cases, but at least we know
+ * in advance how many digits we will convert. We can run the
+ * conversion in steps of DIGIT_GROUP digits, so as to
+ * have many fewer mp_int divisions.
+ */
+ return StrictBignumConversion(&d, convType, bw, b2, s2, s5,
+ k, len, ilim, ilim1, decpt, endPtr);
+ }
+ }
}
+
/*
*----------------------------------------------------------------------
@@ -2237,6 +4380,11 @@ TclInitDoubleConversion(void)
for (i=0; i<8; ++i) {
mp_sqr(pow5+i, pow5+i+1);
}
+ mp_init_set_int(pow5_13, 1220703125);
+ for (i = 1; i < 5; ++i) {
+ mp_init(pow5_13 + i);
+ mp_sqr(pow5_13 + i - 1, pow5_13 + i);
+ }
/*
* Determine the number of decimal digits to the left and right of the
@@ -2293,7 +4441,7 @@ TclFinalizeDoubleConversion(void)
{
int i;
- Tcl_Free((char *) pow10_wide);
+ ckfree((char *) pow10_wide);
for (i=0; i<9; ++i) {
mp_clear(pow5 + i);
}
@@ -2433,6 +4581,20 @@ TclBignumToDouble(
return -r;
}
}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclCeil --
+ *
+ * Computes the smallest floating point number that is at least the
+ * mp_int argument.
+ *
+ * Results:
+ * Returns the floating point number.
+ *
+ *-----------------------------------------------------------------------------
+ */
double
TclCeil(
@@ -2476,6 +4638,20 @@ TclCeil(
mp_clear(&b);
return r;
}
+
+/*
+ *-----------------------------------------------------------------------------
+ *
+ * TclFloor --
+ *
+ * Computes the largest floating point number less than or equal to
+ * the mp_int argument.
+ *
+ * Results:
+ * Returns the floating point value.
+ *
+ *-----------------------------------------------------------------------------
+ */
double
TclFloor(
generated by cgit v0.12 at 2024-11-16 14:54:02 (GMT)