summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--ChangeLog28
-rw-r--r--generic/tclInt.decls7
-rw-r--r--generic/tclInt.h29
-rw-r--r--generic/tclIntDecls.h8
-rwxr-xr-xgeneric/tclStrToD.c2772
-rw-r--r--generic/tclStubInit.c5
-rw-r--r--generic/tclTest.c105
-rw-r--r--generic/tclTomMath.decls8
-rw-r--r--generic/tclTomMathDecls.h14
-rw-r--r--generic/tclUtil.c167
-rw-r--r--tests/util.test803
-rw-r--r--unix/Makefile.in15
-rw-r--r--win/Makefile.in4
-rw-r--r--win/makefile.vc4
14 files changed, 3559 insertions, 410 deletions
diff --git a/ChangeLog b/ChangeLog
index 2af363c..a0acecc 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,31 @@
+2010-11-29 Kevin B. Kenny <kennykb@acm.org>
+
+ * generic/tclInt.decls:
+ * generic/tclInt.h:
+ * generic/tclStrToD.c:
+ * generic/tclTest.c:
+ * generic/tclTomMath.decls:
+ * generic/tclUtil.c:
+ * tests/util.test:
+ * unix/Makefile.in:
+ * win/Makefile.in:
+ * win/makefile.vc: Rewrite of Tcl_PrintDouble and TclDoubleDigits
+ that (a) fixes a severe performance problem with floating point
+ shimmering reported by Karl Lehenbauer, (b) allows TclDoubleDigits
+ to generate the digit strings for 'e' and 'f' format, so that it
+ can be used for tcl_precision != 0 (and possibly later for [format]),
+ (c) fixes [Bug 3120139] by making TclPrintDouble inherently
+ locale-independent, (d) adds test cases to util.test for
+ correct rounding in difficult cases of TclDoubleDigits where fixed-
+ precision results are requested. (e) adds test cases to util.test for
+ the controversial aspects of [Bug 3105247]. As a side effect, two
+ more modules from libtommath (bn_mp_set_int.c and bn_mp_init_set_int.c)
+ are brought into the build, since the new code uses them.
+
+ * generic/tclIntDecls.h:
+ * generic/tclStubInit.c:
+ * generic/tclTomMathDecls.h: Regenerated.
+
2010-11-24 Donal K. Fellows <dkf@users.sf.net>
* tests/chanio.test, tests/iogt.test, tests/ioTrans.test: Convert more
diff --git a/generic/tclInt.decls b/generic/tclInt.decls
index 1a114f1..75780d5 100644
--- a/generic/tclInt.decls
+++ b/generic/tclInt.decls
@@ -13,7 +13,7 @@
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
-# RCS: @(#) $Id: tclInt.decls,v 1.150 2010/10/02 00:23:44 hobbs Exp $
+# RCS: @(#) $Id: tclInt.decls,v 1.151 2010/11/28 23:20:10 kennykb Exp $
library tcl
@@ -996,6 +996,11 @@ declare 248 {
int TclCopyChannel(Tcl_Interp *interp, Tcl_Channel inChan,
Tcl_Channel outChan, Tcl_WideInt toRead, Tcl_Obj *cmdPtr)
}
+
+declare 249 {
+ char* TclDoubleDigits(double dv, int ndigits, int flags,
+ int* decpt, int* signum, char** endPtr)
+}
##############################################################################
diff --git a/generic/tclInt.h b/generic/tclInt.h
index 218c40b..8c5d863 100644
--- a/generic/tclInt.h
+++ b/generic/tclInt.h
@@ -15,7 +15,7 @@
* See the file "license.terms" for information on usage and redistribution of
* this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclInt.h,v 1.486 2010/11/15 21:34:54 andreas_kupries Exp $
+ * RCS: @(#) $Id: tclInt.h,v 1.487 2010/11/28 23:20:11 kennykb Exp $
*/
#ifndef _TCLINT
@@ -2795,6 +2795,32 @@ struct Tcl_LoadHandle_ {
/* Procedure that unloads a loaded module */
};
+/* Flags for conversion of doubles to digit strings */
+
+#define TCL_DD_SHORTEST 0x4
+ /* Use the shortest possible string */
+#define TCL_DD_STEELE 0x5
+ /* Use the original Steele&White algorithm */
+#define TCL_DD_E_FORMAT 0x2
+ /* Use a fixed-length string of digits,
+ * suitable for E format*/
+#define TCL_DD_F_FORMAT 0x3
+ /* Use a fixed number of digits after the
+ * decimal point, suitable for F format */
+
+#define TCL_DD_SHORTEN_FLAG 0x4
+ /* Allow return of a shorter digit string
+ * if it converts losslessly */
+#define TCL_DD_NO_QUICK 0x8
+ /* Debug flag: forbid quick FP conversion */
+
+#define TCL_DD_CONVERSION_TYPE_MASK 0x3
+ /* Mask to isolate the conversion type */
+#define TCL_DD_STEELE0 0x1
+ /* 'Steele&White' after masking */
+#define TCL_DD_SHORTEST0 0x0
+ /* 'Shortest possible' after masking */
+
/*
*----------------------------------------------------------------
* Procedures shared among Tcl modules but not used by the outside world:
@@ -2843,7 +2869,6 @@ MODULE_SCOPE void TclContinuationsEnterDerived(Tcl_Obj *objPtr,
MODULE_SCOPE ContLineLoc *TclContinuationsGet(Tcl_Obj *objPtr);
MODULE_SCOPE void TclContinuationsCopy(Tcl_Obj *objPtr,
Tcl_Obj *originObjPtr);
-MODULE_SCOPE int TclDoubleDigits(char *buf, double value, int *signum);
MODULE_SCOPE void TclDeleteNamespaceVars(Namespace *nsPtr);
/* TIP #280 - Modified token based evulation, with line information. */
MODULE_SCOPE int TclEvalEx(Tcl_Interp *interp, const char *script,
diff --git a/generic/tclIntDecls.h b/generic/tclIntDecls.h
index 259229b..ae4046d 100644
--- a/generic/tclIntDecls.h
+++ b/generic/tclIntDecls.h
@@ -11,7 +11,7 @@
* See the file "license.terms" for information on usage and redistribution
* of this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclIntDecls.h,v 1.144 2010/10/02 00:23:45 hobbs Exp $
+ * RCS: @(#) $Id: tclIntDecls.h,v 1.145 2010/11/28 23:20:11 kennykb Exp $
*/
#ifndef _TCLINTDECLS
@@ -596,6 +596,9 @@ EXTERN void TclResetRewriteEnsemble(Tcl_Interp *interp,
EXTERN int TclCopyChannel(Tcl_Interp *interp,
Tcl_Channel inChan, Tcl_Channel outChan,
Tcl_WideInt toRead, Tcl_Obj *cmdPtr);
+/* 249 */
+EXTERN char* TclDoubleDigits(double dv, int ndigits, int flags,
+ int*decpt, int*signum, char**endPtr);
typedef struct TclIntStubs {
int magic;
@@ -850,6 +853,7 @@ typedef struct TclIntStubs {
int (*tclInitRewriteEnsemble) (Tcl_Interp *interp, int numRemoved, int numInserted, Tcl_Obj *const *objv); /* 246 */
void (*tclResetRewriteEnsemble) (Tcl_Interp *interp, int isRootEnsemble); /* 247 */
int (*tclCopyChannel) (Tcl_Interp *interp, Tcl_Channel inChan, Tcl_Channel outChan, Tcl_WideInt toRead, Tcl_Obj *cmdPtr); /* 248 */
+ char* (*tclDoubleDigits) (double dv, int ndigits, int flags, int*decpt, int*signum, char**endPtr); /* 249 */
} TclIntStubs;
#ifdef __cplusplus
@@ -1269,6 +1273,8 @@ extern const TclIntStubs *tclIntStubsPtr;
(tclIntStubsPtr->tclResetRewriteEnsemble) /* 247 */
#define TclCopyChannel \
(tclIntStubsPtr->tclCopyChannel) /* 248 */
+#define TclDoubleDigits \
+ (tclIntStubsPtr->tclDoubleDigits) /* 249 */
#endif /* defined(USE_TCL_STUBS) */
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(
diff --git a/generic/tclStubInit.c b/generic/tclStubInit.c
index b804b9a..d308f34 100644
--- a/generic/tclStubInit.c
+++ b/generic/tclStubInit.c
@@ -8,7 +8,7 @@
* See the file "license.terms" for information on usage and redistribution
* of this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclStubInit.c,v 1.197 2010/09/16 14:49:37 nijtmans Exp $
+ * RCS: @(#) $Id: tclStubInit.c,v 1.198 2010/11/28 23:20:11 kennykb Exp $
*/
#include "tclInt.h"
@@ -305,6 +305,7 @@ static const TclIntStubs tclIntStubs = {
TclInitRewriteEnsemble, /* 246 */
TclResetRewriteEnsemble, /* 247 */
TclCopyChannel, /* 248 */
+ TclDoubleDigits, /* 249 */
};
static const TclIntPlatStubs tclIntPlatStubs = {
@@ -460,6 +461,8 @@ const TclTomMathStubs tclTomMathStubs = {
TclBN_s_mp_mul_digs, /* 58 */
TclBN_s_mp_sqr, /* 59 */
TclBN_s_mp_sub, /* 60 */
+ TclBN_mp_init_set_int, /* 61 */
+ TclBN_mp_set_int, /* 62 */
};
static const TclStubHooks tclStubHooks = {
diff --git a/generic/tclTest.c b/generic/tclTest.c
index b845d72..80dbbae 100644
--- a/generic/tclTest.c
+++ b/generic/tclTest.c
@@ -14,9 +14,11 @@
* See the file "license.terms" for information on usage and redistribution of
* this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclTest.c,v 1.154 2010/09/27 19:42:38 msofer Exp $
+ * RCS: @(#) $Id: tclTest.c,v 1.155 2010/11/28 23:20:11 kennykb Exp $
*/
+#include <math.h>
+
#undef STATIC_BUILD
#ifndef USE_TCL_STUBS
# define USE_TCL_STUBS
@@ -233,6 +235,9 @@ static int TestdelCmd(ClientData dummy,
Tcl_Interp *interp, int argc, const char **argv);
static int TestdelassocdataCmd(ClientData dummy,
Tcl_Interp *interp, int argc, const char **argv);
+static int TestdoubledigitsObjCmd(ClientData dummy,
+ Tcl_Interp* interp,
+ int objc, Tcl_Obj* const objv[]);
static int TestdstringCmd(ClientData dummy,
Tcl_Interp *interp, int argc, const char **argv);
static int TestencodingObjCmd(ClientData dummy,
@@ -569,6 +574,8 @@ Tcltest_Init(
Tcl_CreateCommand(interp, "testdel", TestdelCmd, NULL, NULL);
Tcl_CreateCommand(interp, "testdelassocdata", TestdelassocdataCmd,
NULL, NULL);
+ Tcl_CreateObjCommand(interp, "testdoubledigits", TestdoubledigitsObjCmd,
+ NULL, NULL);
Tcl_DStringInit(&dstring);
Tcl_CreateCommand(interp, "testdstring", TestdstringCmd, NULL,
NULL);
@@ -1594,6 +1601,102 @@ TestdelassocdataCmd(
}
/*
+ *-----------------------------------------------------------------------------
+ *
+ * TestdoubledigitsCmd --
+ *
+ * This procedure implements the 'testdoubledigits' command. It is
+ * used to test the low-level floating-point formatting primitives
+ * in Tcl.
+ *
+ * Usage:
+ * testdoubledigits fpval ndigits type ?shorten"
+ *
+ * Parameters:
+ * fpval - Floating-point value to format.
+ * ndigits - Digit count to request from Tcl_DoubleDigits
+ * type - One of 'shortest', 'Steele', 'e', 'f'
+ * shorten - Indicates that the 'shorten' flag should be passed in.
+ *
+ *-----------------------------------------------------------------------------
+ */
+
+static int
+TestdoubledigitsObjCmd(ClientData unused,
+ /* NULL */
+ Tcl_Interp* interp,
+ /* Tcl interpreter */
+ int objc,
+ /* Parameter count */
+ Tcl_Obj* const objv[])
+ /* Parameter vector */
+{
+ static const char* options[] = {
+ "shortest",
+ "Steele",
+ "e",
+ "f",
+ NULL
+ };
+ static const int types[] = {
+ TCL_DD_SHORTEST,
+ TCL_DD_STEELE,
+ TCL_DD_E_FORMAT,
+ TCL_DD_F_FORMAT
+ };
+
+ const Tcl_ObjType* doubleType;
+ double d;
+ int status;
+ int ndigits;
+ int type;
+ int decpt;
+ int signum;
+ char* str;
+ char* endPtr;
+ Tcl_Obj* strObj;
+ Tcl_Obj* retval;
+
+ if (objc < 4 || objc > 5) {
+ Tcl_WrongNumArgs(interp, 1, objv, "fpval ndigits type ?shorten?");
+ return TCL_ERROR;
+ }
+ status = Tcl_GetDoubleFromObj(interp, objv[1], &d);
+ if (status != TCL_OK) {
+ doubleType = Tcl_GetObjType("double");
+ if (objv[1]->typePtr == doubleType
+ || TclIsNaN(objv[1]->internalRep.doubleValue)) {
+ status = TCL_OK;
+ memcpy(&d, &(objv[1]->internalRep.doubleValue), sizeof(double));
+ }
+ }
+ if (status != TCL_OK
+ || Tcl_GetIntFromObj(interp, objv[2], &ndigits) != TCL_OK
+ || Tcl_GetIndexFromObj(interp, objv[3], options, "conversion type",
+ TCL_EXACT, &type) != TCL_OK) {
+ fprintf(stderr, "bad value? %g\n", d);
+ return TCL_ERROR;
+ }
+ type = types[type];
+ if (objc > 4) {
+ if (strcmp(Tcl_GetString(objv[4]), "shorten")) {
+ Tcl_SetObjResult(interp, Tcl_NewStringObj("bad flag", -1));
+ return TCL_ERROR;
+ }
+ type |= TCL_DD_SHORTEN_FLAG;
+ }
+ str = TclDoubleDigits(d, ndigits, type, &decpt, &signum, &endPtr);
+ strObj = Tcl_NewStringObj(str, endPtr-str);
+ ckfree(str);
+ retval = Tcl_NewListObj(1, &strObj);
+ Tcl_ListObjAppendElement(NULL, retval, Tcl_NewIntObj(decpt));
+ strObj = Tcl_NewStringObj(signum ? "-" : "+", 1);
+ Tcl_ListObjAppendElement(NULL, retval, strObj);
+ Tcl_SetObjResult(interp, retval);
+ return TCL_OK;
+}
+
+/*
*----------------------------------------------------------------------
*
* TestdstringCmd --
diff --git a/generic/tclTomMath.decls b/generic/tclTomMath.decls
index 12fa7cd..0bd5ca5 100644
--- a/generic/tclTomMath.decls
+++ b/generic/tclTomMath.decls
@@ -13,7 +13,7 @@
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
-# RCS: @(#) $Id: tclTomMath.decls,v 1.8 2010/09/15 07:33:54 nijtmans Exp $
+# RCS: @(#) $Id: tclTomMath.decls,v 1.9 2010/11/28 23:20:11 kennykb Exp $
library tcl
@@ -214,3 +214,9 @@ declare 59 {
declare 60 {
int TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c)
}
+declare 61 {
+ int TclBN_mp_init_set_int(mp_int* a, unsigned long i)
+}
+declare 62 {
+ int TclBN_mp_set_int(mp_int* a, unsigned long i)
+} \ No newline at end of file
diff --git a/generic/tclTomMathDecls.h b/generic/tclTomMathDecls.h
index 61e9c07..4e55ce1 100644
--- a/generic/tclTomMathDecls.h
+++ b/generic/tclTomMathDecls.h
@@ -11,7 +11,7 @@
* See the file "license.terms" for information on usage and redistribution
* of this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclTomMathDecls.h,v 1.13 2010/08/19 04:26:04 nijtmans Exp $
+ * RCS: @(#) $Id: tclTomMathDecls.h,v 1.14 2010/11/28 23:20:11 kennykb Exp $
*/
#ifndef _TCLTOMMATHDECLS
@@ -79,6 +79,7 @@
#define mp_init_copy TclBN_mp_init_copy
#define mp_init_multi TclBN_mp_init_multi
#define mp_init_set TclBN_mp_init_set
+#define mp_init_set_int TclBN_mp_init_set_int
#define mp_init_size TclBN_mp_init_size
#define mp_karatsuba_mul TclBN_mp_karatsuba_mul
#define mp_karatsuba_sqr TclBN_mp_karatsuba_sqr
@@ -96,6 +97,7 @@
#define mp_rshd TclBN_mp_rshd
#define mp_s_rmap TclBNMpSRmap
#define mp_set TclBN_mp_set
+#define mp_set_int TclBN_mp_set_int
#define mp_shrink TclBN_mp_shrink
#define mp_sqr TclBN_mp_sqr
#define mp_sqrt TclBN_mp_sqrt
@@ -268,6 +270,10 @@ EXTERN int TclBN_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c,
EXTERN int TclBN_s_mp_sqr(mp_int *a, mp_int *b);
/* 60 */
EXTERN int TclBN_s_mp_sub(mp_int *a, mp_int *b, mp_int *c);
+/* 61 */
+EXTERN int TclBN_mp_init_set_int(mp_int*a, unsigned long i);
+/* 62 */
+EXTERN int TclBN_mp_set_int(mp_int*a, unsigned long i);
typedef struct TclTomMathStubs {
int magic;
@@ -334,6 +340,8 @@ typedef struct TclTomMathStubs {
int (*tclBN_s_mp_mul_digs) (mp_int *a, mp_int *b, mp_int *c, int digs); /* 58 */
int (*tclBN_s_mp_sqr) (mp_int *a, mp_int *b); /* 59 */
int (*tclBN_s_mp_sub) (mp_int *a, mp_int *b, mp_int *c); /* 60 */
+ int (*tclBN_mp_init_set_int) (mp_int*a, unsigned long i); /* 61 */
+ int (*tclBN_mp_set_int) (mp_int*a, unsigned long i); /* 62 */
} TclTomMathStubs;
#ifdef __cplusplus
@@ -472,6 +480,10 @@ extern const TclTomMathStubs *tclTomMathStubsPtr;
(tclTomMathStubsPtr->tclBN_s_mp_sqr) /* 59 */
#define TclBN_s_mp_sub \
(tclTomMathStubsPtr->tclBN_s_mp_sub) /* 60 */
+#define TclBN_mp_init_set_int \
+ (tclTomMathStubsPtr->tclBN_mp_init_set_int) /* 61 */
+#define TclBN_mp_set_int \
+ (tclTomMathStubsPtr->tclBN_mp_set_int) /* 62 */
#endif /* defined(USE_TCL_STUBS) */
diff --git a/generic/tclUtil.c b/generic/tclUtil.c
index 45cc8a1..a04c29c 100644
--- a/generic/tclUtil.c
+++ b/generic/tclUtil.c
@@ -11,7 +11,7 @@
* See the file "license.terms" for information on usage and redistribution of
* this file, and for a DISCLAIMER OF ALL WARRANTIES.
*
- * RCS: @(#) $Id: tclUtil.c,v 1.118 2010/10/01 12:52:50 dkf Exp $
+ * RCS: @(#) $Id: tclUtil.c,v 1.119 2010/11/28 23:20:11 kennykb Exp $
*/
#include "tclInt.h"
@@ -2234,129 +2234,100 @@ Tcl_PrintDouble(
char *p, c;
int exponent;
int signum;
- char buffer[TCL_DOUBLE_SPACE];
+ char* digits;
+ char* end;
Tcl_UniChar ch;
int *precisionPtr = Tcl_GetThreadData(&precisionKey, (int)sizeof(int));
/*
- * If *precisionPtr == 0, then use TclDoubleDigits to develop a decimal
- * significand and exponent, then format it in E or F format as
- * appropriate. If *precisionPtr != 0, use the native sprintf and then add
- * a trailing ".0" if there is no decimal point in the rep.
+ * Handle NaN.
*/
+
+ if (TclIsNaN(value)) {
+ TclFormatNaN(value, dst);
+ return;
+ }
- if (*precisionPtr == 0) {
+ /*
+ * Handle infinities.
+ */
+
+ if (TclIsInfinite(value)) {
/*
- * Handle NaN.
+ * Remember to copy the terminating NUL too.
*/
-
- if (TclIsNaN(value)) {
- TclFormatNaN(value, dst);
- return;
+
+ if (value < 0) {
+ memcpy(dst, "-Inf", 5);
+ } else {
+ memcpy(dst, "Inf", 4);
}
+ return;
+ }
+ /*
+ * Ordinary (normal and denormal) values.
+ */
+
+ if (*precisionPtr == 0) {
+ digits = TclDoubleDigits(value, -1, TCL_DD_SHORTEST,
+ &exponent, &signum, &end);
+ } else {
+ digits = TclDoubleDigits(value, *precisionPtr, TCL_DD_E_FORMAT,
+ &exponent, &signum, &end);
+ }
+ if (signum) {
+ *dst++ = '-';
+ }
+ p = digits;
+ if (exponent < -4 || exponent > 16) {
/*
- * Handle infinities.
+ * E format for numbers < 1e-3 or >= 1e17.
*/
-
- if (TclIsInfinite(value)) {
- /*
- * Remember to copy the terminating NUL too.
- */
-
- if (value < 0) {
- memcpy(dst, "-Inf", 5);
- } else {
- memcpy(dst, "Inf", 4);
+
+ *dst++ = *p++;
+ c = *p;
+ if (c != '\0') {
+ *dst++ = '.';
+ while (c != '\0') {
+ *dst++ = c;
+ c = *++p;
}
- return;
}
-
+ sprintf(dst, "e%+d", exponent);
+ } else {
/*
- * Ordinary (normal and denormal) values.
+ * F format for others.
*/
-
- exponent = TclDoubleDigits(buffer, value, &signum);
- if (signum) {
- *dst++ = '-';
+
+ if (exponent < 0) {
+ *dst++ = '0';
}
- p = buffer;
- if (exponent < -3 || exponent > 17) {
- /*
- * E format for numbers < 1e-3 or >= 1e17.
- */
-
- *dst++ = *p++;
- c = *p;
+ c = *p;
+ while (exponent-- >= 0) {
if (c != '\0') {
- *dst++ = '.';
- while (c != '\0') {
- *dst++ = c;
- c = *++p;
- }
- }
- sprintf(dst, "e%+d", exponent-1);
- } else {
- /*
- * F format for others.
- */
-
- if (exponent <= 0) {
- *dst++ = '0';
- }
- c = *p;
- while (exponent-- > 0) {
- if (c != '\0') {
- *dst++ = c;
- c = *++p;
- } else {
- *dst++ = '0';
- }
- }
- *dst++ = '.';
- if (c == '\0') {
- *dst++ = '0';
+ *dst++ = c;
+ c = *++p;
} else {
- while (++exponent < 0) {
- *dst++ = '0';
- }
- while (c != '\0') {
- *dst++ = c;
- c = *++p;
- }
+ *dst++ = '0';
}
- *dst++ = '\0';
}
- } else {
- /*
- * tcl_precision is supplied, pass it to the native sprintf.
- */
-
- sprintf(dst, "%.*g", *precisionPtr, value);
-
- /*
- * If the ASCII result looks like an integer, add ".0" so that it
- * doesn't look like an integer anymore. This prevents floating-point
- * values from being converted to integers unintentionally. Check for
- * ASCII specifically to speed up the function.
- */
-
- for (p = dst; *p != 0;) {
- if (UCHAR(*p) < 0x80) {
- c = *p++;
- } else {
- p += Tcl_UtfToUniChar(p, &ch);
- c = UCHAR(ch);
+ *dst++ = '.';
+ if (c == '\0') {
+ *dst++ = '0';
+ } else {
+ while (++exponent < -1) {
+ *dst++ = '0';
}
- if ((c == '.') || isalpha(UCHAR(c))) { /* INTL: ISO only. */
- return;
+ while (c != '\0') {
+ *dst++ = c;
+ c = *++p;
}
}
- p[0] = '.';
- p[1] = '0';
- p[2] = 0;
+ *dst++ = '\0';
}
+ ckfree(digits);
}
/*
diff --git a/tests/util.test b/tests/util.test
index 994fc0f..1d31609 100644
--- a/tests/util.test
+++ b/tests/util.test
@@ -7,13 +7,14 @@
# See the file "license.terms" for information on usage and redistribution
# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
-# RCS: @(#) $Id: util.test,v 1.20 2008/10/14 16:35:44 dgp Exp $
+# RCS: @(#) $Id: util.test,v 1.21 2010/11/28 23:20:11 kennykb Exp $
if {[lsearch [namespace children] ::tcltest] == -1} {
package require tcltest
namespace import -force ::tcltest::*
}
+testConstraint controversialNaN 1
testConstraint testdstring [llength [info commands testdstring]]
testConstraint testconcatobj [llength [info commands testconcatobj]]
@@ -43,6 +44,10 @@ proc testIEEE {} {
ieeeValues(+Infinity)
binary scan \x00\x00\x00\x00\x00\x00\xf8\x7f d \
ieeeValues(NaN)
+ binary scan \x00\x00\x00\x00\x00\x00\xf8\xff d \
+ ieeeValues(-NaN)
+ binary scan \xef\xcd\xab\x89\x67\x45\xfb\xff d \
+ ieeeValues(-NaN(3456789abcdef))
set ieeeValues(littleEndian) 1
return 1
}
@@ -65,6 +70,10 @@ proc testIEEE {} {
ieeeValues(+Infinity)
binary scan \x7f\xf8\x00\x00\x00\x00\x00\x00 d \
ieeeValues(NaN)
+ binary scan \xff\xf8\x00\x00\x00\x00\x00\x00 d \
+ ieeeValues(-NaN)
+ binary scan \xff\xfb\x45\x67\x89\xab\xcd\xef d \
+ ieeeValues(-NaN(3456789abcdef))
set ieeeValues(littleEndian) 0
return 1
}
@@ -85,6 +94,30 @@ proc convertDouble { x } {
return $result
}
+proc verdonk_test {sig binexp shouldbe exp} {
+ regexp {([-+]?)([0-9a-f]+)} $sig -> signum sig
+ scan $sig %llx sig
+ if {$signum eq {-}} {
+ set signum [expr 1<<63]
+ } else {
+ set signum 0
+ }
+ regexp {E([-+]?[0-9]+)} $binexp -> binexp
+ set word [expr {$signum | (($binexp + 0x3ff)<<52)|($sig & ~(1<<52))}]
+ binary scan [binary format w $word] q double
+ regexp {([-+])(\d+)_(\d+)\&} $shouldbe -> signum digits1 digits2
+ regexp {E([-+]\d+)} $exp -> decexp
+ incr decexp [expr {[string length $digits1] - 1}]
+ lassign [testdoubledigits $double [string length $digits1] e] \
+ outdigits decpt outsign
+ if {[string index $digits2 0] >= 5} {
+ incr digits1
+ }
+ if {$outsign != $signum || $outdigits != $digits1 || $decpt != $decexp} {
+ return -code error "result is ${outsign}0.${outdigits}E$decpt\
+ should be ${signum}0.${digits1}E$decexp"
+ }
+}
test util-1.1 {TclFindElement procedure - binary element in middle of list} {
lindex {0 foo\x00help 1} 1
@@ -1106,6 +1139,774 @@ test util-11.23 {Tcl_PrintDouble - scaling} {
expr 1.1e17
} {1.1e+17}
+test util-12.1 {TclDoubleDigits - Inf} ieeeFloatingPoint {
+ testdoubledigits Inf -1 shortest
+} {Infinity 9999 +}
+test util-12.2 {TclDoubleDigits - -Inf} ieeeFloatingPoint {
+ testdoubledigits -Inf -1 shortest
+} {Infinity 9999 -}
+test util-12.3 {TclDoubleDigits - NaN} ieeeFloatingPoint {
+ testdoubledigits $ieeeValues(NaN) -1 shortest
+} {NaN 9999 +}
+test util-12.4 {TclDoubleDigits - NaN} {*}{
+ -constraints {ieeeFloatingPoint && controversialNaN}
+ -body {
+ testdoubledigits -NaN -1 shortest
+ }
+ -result {NaN 9999 -}
+}
+test util-12.5 {TclDoubleDigits - 0} {
+ testdoubledigits 0.0 -1 shortest
+} {0 0 +}
+test util-12.6 {TclDoubleDigits - -0} {
+ testdoubledigits -0.0 -1 shortest
+} {0 0 -}
+
+# Verdonk test vectors
+
+test util-13.1 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1754e31cd072da E+1008 +4_000000000000000000& E+303
+ }
+ -result {}
+}
+test util-13.2 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1afcef51f0fb5f E+265 -1_000000000000000000& E+80
+ }
+ -result {}
+}
+test util-13.3 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1754e31cd072da E+1006 +1_000000000000000000& E+303
+ }
+ -result {}
+}
+test util-13.4 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1754e31cd072da E+1007 -2_000000000000000000& E+303
+ }
+ -result {}
+}
+test util-13.5 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1e07b27dd78b14 E-848 +1_00000000000000000& E-255
+ }
+ -result {}
+}
+test util-13.6 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1e29e9c56687fe E-709 -7_00000000000000000& E-214
+ }
+ -result {}
+}
+test util-13.7 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1be03d0bf225c7 E-137 +1_00000000000000000& E-41
+ }
+ -result {}
+}
+test util-13.8 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1a2fe76a3f9475 E-499 -1_00000000000000000& E-150
+ }
+ -result {}
+}
+test util-13.9 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 19a2028368022e E+1019 +8_999999999999999999& E+306
+ }
+ -result {}
+}
+test util-13.10 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1317e5ef3ab327 E+509 -1_999999999999999999& E+153
+ }
+ -result {}
+}
+test util-13.11 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1317e5ef3ab327 E+510 +3_99999999999999999& E+153
+ }
+ -result {}
+}
+test util-13.12 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1317e5ef3ab327 E+511 -7_99999999999999999& E+153
+ }
+ -result {}
+}
+test util-13.13 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1eb8e84fa0b278 E-1008 +6_999999999999999999& E-304
+ }
+ -result {}
+}
+test util-13.14 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -13339131c46f8b E-1004 -6_999999999999999999& E-303
+ }
+ -result {}
+}
+test util-13.15 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1c0f92a6276c9d E-162 +2_999999999999999999& E-49
+ }
+ -result {}
+}
+test util-13.16 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -15ce1f143d7ad2 E-443 -5_99999999999999999& E-134
+ }
+ -result {}
+}
+test util-13.17 {just over exact - 2 digits} {*}{
+ -body {
+ verdonk_test 1c0794d9d40e96 E-301 +43_000000000000000000& E-92
+ }
+ -result {}
+}
+test util-13.18 {just over exact - 2 digits} {*}{
+ -body {
+ verdonk_test -1c0794d9d40e96 E-300 -86_000000000000000000& E-92
+ }
+ -result {}
+}
+test util-13.19 {just over exact - 2 digits} {*}{
+ -body {
+ verdonk_test 1cd5bee57763e6 E-241 +51_000000000000000000& E-74
+ }
+ -result {}
+}
+test util-13.20 {just under exact - 2 digits} {*}{
+ -body {
+ verdonk_test 1d1c26db7d0dae E+651 +16_999999999999999999& E+195
+ }
+ -result {}
+}
+test util-13.21 {just under exact - 2 digits} {*}{
+ -body {
+ verdonk_test -13f7ced916872b E-5 -38_999999999999999999& E-3
+ }
+ -result {}
+}
+test util-13.22 {just over exact - 3 digits} {*}{
+ -body {
+ verdonk_test 17d93193f78fc6 E+588 +151_0000000000000000000& E+175
+ }
+ -result {}
+}
+test util-13.23 {just over exact - 3 digits} {*}{
+ -body {
+ verdonk_test -1a82a1631eeb30 E-625 -119_000000000000000000& E-190
+ }
+ -result {}
+}
+test util-13.24 {just under exact - 3 digits} {*}{
+ -body {
+ verdonk_test -16c309024bab4b E+290 -282_999999999999999999& E+85
+ }
+ -result {}
+}
+test util-13.25 {just over exact - 8 digits} {*}{
+ -body {
+ verdonk_test 1dbbac6f83a821 E-800 +27869147_0000000000000000000& E-248
+ }
+ -result {}
+}
+test util-13.26 {just under exact - 9 digits} {*}{
+ -body {
+ verdonk_test -1c569e968e0944 E+430 -491080653_9999999999999999999& E+121
+ }
+ -result {}
+}
+test util-13.27 {just under exact - 9 digits} {*}{
+ -body {
+ verdonk_test 1c569e968e0944 E+429 +245540326_9999999999999999999& E+121
+ }
+ -result {}
+}
+test util-13.28 {just over exact - 10 digits} {*}{
+ -body {
+ verdonk_test -1fc575867314ee E-330 -9078555839_0000000000000000000& E-109
+ }
+ -result {}
+}
+test util-13.29 {just under exact - 10 digits} {*}{
+ -body {
+ verdonk_test -1c569e968e0944 E+428 -1227701634_9999999999999999999& E+120
+ }
+ -result {}
+}
+test util-13.30 {just over exact - 11 digits} {*}{
+ -body {
+ verdonk_test 1fc575867314ee E-329 +18157111678_0000000000000000000& E-109
+ }
+ -result {}
+}
+test util-13.31 {just over exact - 14 digits} {*}{
+ -body {
+ verdonk_test -18bf7e7fa6f02a E-196 -15400733123779_0000000000000000000& E-72
+ }
+ -result {}
+}
+test util-13.32 {just over exact - 17 digits} {*}{
+ -body {
+ verdonk_test -13de005bd620df E+217 -26153245263757307_0000000000000000000& E+49
+ }
+ -result {}
+}
+test util-13.33 {just over exact - 18 digits} {*}{
+ -body {
+ verdonk_test 1f92bacb3cb40c E+718 +272104041512242479_0000000000000000000& E+199
+ }
+ -result {}
+}
+test util-13.34 {just over exact - 18 digits} {*}{
+ -body {
+ verdonk_test -1f92bacb3cb40c E+719 -544208083024484958_0000000000000000000& E+199
+ }
+ -result {}
+}
+test util-13.35 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 142dbf25096cf5 E+148 +4_500000000000000000& E+44
+ }
+ -result {}
+}
+test util-13.36 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1afcef51f0fb5f E+263 -2_500000000000000000& E+79
+ }
+ -result {}
+}
+test util-13.37 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 102498ea6df0c4 E+145 +4_500000000000000000& E+43
+ }
+ -result {}
+}
+test util-13.38 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1754e31cd072da E+1004 -2_500000000000000000& E+302
+ }
+ -result {}
+}
+test util-13.39 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 12deac01e2b4f7 E-557 +2_50000000000000000& E-168
+ }
+ -result {}
+}
+test util-13.40 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1b1df536c13eee E-307 -6_50000000000000000& E-93
+ }
+ -result {}
+}
+test util-13.41 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 10711fed5b19a4 E-154 +4_50000000000000000& E-47
+ }
+ -result {}
+}
+test util-13.42 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -148d67e8b1e00d E-151 -4_50000000000000000& E-46
+ }
+ -result {}
+}
+test util-13.43 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 1c8c574c0c6be7 E+187 +3_49999999999999999& E+56
+ }
+ -result {}
+}
+test util-13.44 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1756183c147514 E+206 -1_49999999999999999& E+62
+ }
+ -result {}
+}
+test util-13.45 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 12ab469676c410 E+203 +1_49999999999999999& E+61
+ }
+ -result {}
+}
+test util-13.46 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1539684e774b48 E+246 -1_49999999999999999& E+74
+ }
+ -result {}
+}
+test util-13.47 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 12e5f5dfa4fe9d E-286 +9_499999999999999999& E-87
+ }
+ -result {}
+}
+test util-13.48 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1bdc2417bf7787 E-838 -9_499999999999999999& E-253
+ }
+ -result {}
+}
+test util-13.49 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 1eb8e84fa0b278 E-1009 +3_499999999999999999& E-304
+ }
+ -result {}
+}
+test util-13.50 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1e3cbc9907fdc8 E-290 -9_499999999999999999& E-88
+ }
+ -result {}
+}
+test util-13.51 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test 10ad836f269a17 E-324 +30_500000000000000000& E-99
+ }
+ -result {}
+}
+test util-13.52 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test -1b39ae1909c31b E-687 -26_500000000000000000& E-208
+ }
+ -result {}
+}
+test util-13.53 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test 1b2ab18615fcc6 E-576 +686_500000000000000000& E-176
+ }
+ -result {}
+}
+test util-13.54 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test -13e1f90a573064 E-624 -178_500000000000000000& E-190
+ }
+ -result {}
+}
+test util-13.55 {just under half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test 16c309024bab4b E+289 +141_499999999999999999& E+85
+ }
+ -result {}
+}
+test util-13.56 {just under half ulp - 4 digits} {*}{
+ -body {
+ verdonk_test -159bd3ad46e346 E+193 -1695_499999999999999999& E+55
+ }
+ -result {}
+}
+test util-13.57 {just under half ulp - 4 digits} {*}{
+ -body {
+ verdonk_test 1df4170f0fdecc E+124 +3981_499999999999999999& E+34
+ }
+ -result {}
+}
+test util-13.58 {just over half ulp - 6 digits} {*}{
+ -body {
+ verdonk_test 17e1e0f1c7a4ac E+415 +126300_5000000000000000000& E+120
+ }
+ -result {}
+}
+test util-13.59 {just over half ulp - 6 digits} {*}{
+ -body {
+ verdonk_test -1dda592e398dd7 E+418 -126300_5000000000000000000& E+121
+ }
+ -result {}
+}
+test util-13.60 {just under half ulp - 7 digits} {*}{
+ -body {
+ verdonk_test -1e597c0b94b7ae E+453 -4411845_499999999999999999& E+130
+ }
+ -result {}
+}
+test util-13.61 {just under half ulp - 9 digits} {*}{
+ -body {
+ verdonk_test 1c569e968e0944 E+427 +613850817_4999999999999999999& E+120
+ }
+ -result {}
+}
+test util-13.62 {just under half ulp - 9 digits} {*}{
+ -body {
+ verdonk_test -1c569e968e0944 E+428 -122770163_49999999999999999999& E+121
+ }
+ -result {}
+}
+test util-13.63 {just over half ulp - 18 digits} {*}{
+ -body {
+ verdonk_test 17ae0c186d8709 E+719 +408156062268363718_5000000000000000000& E+199
+ }
+ -result {}
+}
+test util-13.64 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 152d02c7e14af7 E+76 +1_0000000000000000& E+23
+ }
+ -result {}
+}
+test util-13.65 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -19d971e4fe8402 E+89 -1_0000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.66 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 19d971e4fe8402 E+90 +2_0000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.67 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -19d971e4fe8402 E+91 -4_0000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.68 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 15798ee2308c3a E-27 +1_0000000000000000& E-8
+ }
+ -result {}
+}
+test util-13.69 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -15798ee2308c3a E-26 -2_0000000000000000& E-8
+ }
+ -result {}
+}
+test util-13.70 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test 15798ee2308c3a E-25 +4_0000000000000000& E-8
+ }
+ -result {}
+}
+test util-13.71 {just over exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1ef2d0f5da7dd9 E-84 -1_0000000000000000& E-25
+ }
+ -result {}
+}
+test util-13.72 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 1a784379d99db4 E+78 +4_9999999999999999& E+23
+ }
+ -result {}
+}
+test util-13.73 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1a784379d99db4 E+80 -1_9999999999999999& E+24
+ }
+ -result {}
+}
+test util-13.74 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 13da329b633647 E+81 +2_9999999999999999& E+24
+ }
+ -result {}
+}
+test util-13.75 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -1cf389cd46047d E+85 -6_9999999999999999& E+25
+ }
+ -result {}
+}
+test util-13.76 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 19999999999999 E-3 +1_99999999999999999& E-1
+ }
+ -result {}
+}
+test util-13.77 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -13333333333333 E-2 -2_99999999999999999& E-1
+ }
+ -result {}
+}
+test util-13.78 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test 16849b86a12b9b E-48 +4_99999999999999999& E-15
+ }
+ -result {}
+}
+test util-13.79 {just under exact - 1 digits} {*}{
+ -body {
+ verdonk_test -16849b86a12b9b E-46 -1_99999999999999999& E-14
+ }
+ -result {}
+}
+test util-13.80 {just over exact - 2 digits} {*}{
+ -body {
+ verdonk_test 17ccfc73126788 E-71 +63_00000000000000000& E-23
+ }
+ -result {}
+}
+test util-13.81 {just over exact - 2 digits} {*}{
+ -body {
+ verdonk_test -1dc03b8fd7016a E-68 -63_00000000000000000& E-22
+ }
+ -result {}
+}
+test util-13.82 {just under exact - 2 digits} {*}{
+ -body {
+ verdonk_test 13f7ced916872b E-5 +38_999999999999999999& E-3
+ }
+ -result {}
+}
+test util-13.83 {just over exact - 3 digits} {*}{
+ -body {
+ verdonk_test 1b297cad9f70b6 E+97 +269_000000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.84 {just over exact - 3 digits} {*}{
+ -body {
+ verdonk_test -1b297cad9f70b6 E+98 -538_00000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.85 {just over exact - 3 digits} {*}{
+ -body {
+ verdonk_test 1cdc06b20ef183 E-82 +373_00000000000000000& E-27
+ }
+ -result {}
+}
+test util-13.86 {just over exact - 4 digits} {*}{
+ -body {
+ verdonk_test 1b297cad9f70b6 E+96 +1345_00000000000000000& E+26
+ }
+ -result {}
+}
+# this one is not 4 digits, it is 3, and it is covered above.
+test util-13.87 {just over exact - 4 digits} {*}{
+ -constraints knownBadTest
+ -body {
+ verdonk_test -1b297cad9f70b6 E+97 -2690_00000000000000000& E+26
+ }
+ -result {}
+}
+test util-13.88 {just over exact - 5 digits} {*}{
+ -body {
+ verdonk_test -150a246ecd44f3 E-63 -14257_00000000000000000& E-23
+ }
+ -result {}
+}
+test util-13.89 {just under exact - 6 digits} {*}{
+ -body {
+ verdonk_test -119b96f36ec68b E-19 -209900_999999999999999999& E-11
+ }
+ -result {}
+}
+test util-13.90 {just over exact - 11 digits} {*}{
+ -body {
+ verdonk_test 1c06d366394441 E-35 +50980203373_000000000000000000& E-21
+ }
+ -result {}
+}
+test util-13.91 {just under exact - 12 digits} {*}{
+ -body {
+ verdonk_test -1f58ac4db68c90 E+122 -104166211810_99999999999999999& E+26
+ }
+ -result {}
+}
+test util-13.92 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 19d971e4fe8402 E+87 +2_5000000000000000& E+26
+ }
+ -result {}
+}
+test util-13.93 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1dc74be914d16b E+81 -4_500000000000000& E+24
+ }
+ -result {}
+}
+test util-13.94 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 14adf4b7320335 E+84 +2_500000000000000& E+25
+ }
+ -result {}
+}
+test util-13.95 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1ae22487c1042b E+85 -6_5000000000000000& E+25
+ }
+ -result {}
+}
+test util-13.96 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 187fe49aab41e0 E-54 +8_5000000000000000& E-17
+ }
+ -result {}
+}
+test util-13.97 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1f5c05e4b23fd7 E-61 -8_5000000000000000& E-19
+ }
+ -result {}
+}
+test util-13.98 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 1faa7ab552a552 E-42 +4_5000000000000000& E-13
+ }
+ -result {}
+}
+test util-13.99 {just over half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1b7cdfd9d7bdbb E-36 -2_5000000000000000& E-11
+ }
+ -result {}
+}
+test util-13.100 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 13da329b633647 E+80 +1_4999999999999999& E+24
+ }
+ -result {}
+}
+test util-13.101 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1cf389cd46047d E+84 -3_49999999999999999& E+25
+ }
+ -result {}
+}
+test util-13.102 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 1f04ef12cb04cf E+85 +7_4999999999999999& E+25
+ }
+ -result {}
+}
+test util-13.103 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -1f04ef12cb04cf E+86 -1_4999999999999999& E+26
+ }
+ -result {}
+}
+test util-13.104 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 13333333333333 E-3 +1_49999999999999999& E-1
+ }
+ -result {}
+}
+test util-13.105 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -107e1fe91b0b70 E-36 -1_49999999999999999& E-11
+ }
+ -result {}
+}
+test util-13.106 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test 149da7e361ce4c E-33 +1_49999999999999999& E-10
+ }
+ -result {}
+}
+test util-13.107 {just under half ulp - 1 digits} {*}{
+ -body {
+ verdonk_test -19c511dc3a41df E-30 -1_49999999999999999& E-9
+ }
+ -result {}
+}
+test util-13.108 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test -1aa83d74267822 E+93 -16_5000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.109 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test 18f1d5969453de E+89 +96_5000000000000000& E+25
+ }
+ -result {}
+}
+test util-13.110 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test 11d9bd564dcda6 E-70 +94_50000000000000000& E-23
+ }
+ -result {}
+}
+test util-13.111 {just over half ulp - 2 digits} {*}{
+ -body {
+ verdonk_test -1a58973ecbede6 E-48 -58_50000000000000000& E-16
+ }
+ -result {}
+}
+test util-13.112 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test 1b297cad9f70b6 E+95 +672_50000000000000000& E+26
+ }
+ -result {}
+}
+test util-13.113 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test -1b297cad9f70b6 E+96 -134_500000000000000000& E+27
+ }
+ -result {}
+}
+test util-13.114 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test 1cdc06b20ef183 E-83 +186_50000000000000000& E-27
+ }
+ -result {}
+}
+test util-13.115 {just over half ulp - 3 digits} {*}{
+ -body {
+ verdonk_test -136071dcae4565 E-47 -860_50000000000000000& E-17
+ }
+ -result {}
+}
+test util-13.116 {just over half ulp - 6 digits} {*}{
+ -body {
+ verdonk_test 1cb968d297dde8 E+99 +113788_50000000000000000& E+25
+ }
+ -result {}
+}
+test util-13.117 {just over half ulp - 6 digits} {*}{
+ -body {
+ verdonk_test -11f3e1839eeab1 E+103 -113788_50000000000000000& E+26
+ }
+ -result {}
+}
+test util-13.118 {just under half ulp - 9 digits} {*}{
+ -body {
+ verdonk_test 1e9cec176c96f8 E+117 +317903333_49999999999999999& E+27
+ }
+ -result {}
+}
+test util-13.119 {just over half ulp - 11 digits} {*}{
+ -body {
+ verdonk_test 1c06d366394441 E-36 +25490101686_500000000000000000& E-21
+ }
+ -result {}
+}
+test util-13.120 {just under half ulp - 11 digits} {*}{
+ -body {
+ verdonk_test 1f58ac4db68c90 E+121 +52083105905_49999999999999999& E+26
+ }
+ -result {}
+}
+
+test util-14.1 {funky NaN} {*}{
+ -constraints {ieeeFloatingPoint && controversialNaN}
+ -body {
+ set ieeeValues(-NaN)
+ }
+ -result -NaN
+}
+
+test util-14.2 {funky NaN} {*}{
+ -constraints {ieeeFloatingPoint && controversialNaN}
+ -body {
+ set ieeeValues(-NaN(3456789abcdef))
+ }
+ -result -NaN(3456789abcdef)
+}
+
# cleanup
::tcltest::cleanupTests
return
+
+# Local Variables:
+# mode: tcl
+# End: \ No newline at end of file
diff --git a/unix/Makefile.in b/unix/Makefile.in
index fad07e6..8eebff2 100644
--- a/unix/Makefile.in
+++ b/unix/Makefile.in
@@ -4,7 +4,7 @@
# "./configure", which is a configuration script generated by the "autoconf"
# program (constructs like "@foo@" will get replaced in the actual Makefile.
#
-# RCS: @(#) $Id: Makefile.in,v 1.310 2010/11/04 15:55:45 dkf Exp $
+# RCS: @(#) $Id: Makefile.in,v 1.311 2010/11/28 23:20:11 kennykb Exp $
VERSION = @TCL_VERSION@
MAJOR_VERSION = @TCL_MAJOR_VERSION@
@@ -322,12 +322,13 @@ TOMMATH_OBJS = bncore.o bn_reverse.o bn_fast_s_mp_mul_digs.o \
bn_mp_div_2d.o bn_mp_div_3.o \
bn_mp_exch.o bn_mp_expt_d.o bn_mp_grow.o bn_mp_init.o \
bn_mp_init_copy.o bn_mp_init_multi.o bn_mp_init_set.o \
- bn_mp_init_size.o bn_mp_karatsuba_mul.o \
+ bn_mp_init_set_int.o bn_mp_init_size.o bn_mp_karatsuba_mul.o \
bn_mp_karatsuba_sqr.o \
bn_mp_lshd.o bn_mp_mod.o bn_mp_mod_2d.o bn_mp_mul.o bn_mp_mul_2.o \
bn_mp_mul_2d.o bn_mp_mul_d.o bn_mp_neg.o bn_mp_or.o \
bn_mp_radix_size.o bn_mp_radix_smap.o \
- bn_mp_read_radix.o bn_mp_rshd.o bn_mp_set.o bn_mp_shrink.o \
+ bn_mp_read_radix.o bn_mp_rshd.o bn_mp_set.o bn_mp_set_int.o \
+ bn_mp_shrink.o \
bn_mp_sqr.o bn_mp_sqrt.o bn_mp_sub.o bn_mp_sub_d.o \
bn_mp_to_unsigned_bin.o bn_mp_to_unsigned_bin_n.o \
bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_toradix_n.o \
@@ -496,6 +497,7 @@ TOMMATH_SRCS = \
$(TOMMATH_DIR)/bn_mp_init_copy.c \
$(TOMMATH_DIR)/bn_mp_init_multi.c \
$(TOMMATH_DIR)/bn_mp_init_set.c \
+ $(TOMMATH_DIR)/bn_mp_init_set_int.c \
$(TOMMATH_DIR)/bn_mp_init_size.c \
$(TOMMATH_DIR)/bn_mp_karatsuba_mul.c \
$(TOMMATH_DIR)/bn_mp_karatsuba_sqr.c \
@@ -513,6 +515,7 @@ TOMMATH_SRCS = \
$(TOMMATH_DIR)/bn_mp_read_radix.c \
$(TOMMATH_DIR)/bn_mp_rshd.c \
$(TOMMATH_DIR)/bn_mp_set.c \
+ $(TOMMATH_DIR)/bn_mp_set_int.c \
$(TOMMATH_DIR)/bn_mp_shrink.c \
$(TOMMATH_DIR)/bn_mp_sqr.c \
$(TOMMATH_DIR)/bn_mp_sqrt.c \
@@ -1380,6 +1383,9 @@ bn_mp_init_multi.o: $(TOMMATH_DIR)/bn_mp_init_multi.c $(MATHHDRS)
bn_mp_init_set.o: $(TOMMATH_DIR)/bn_mp_init_set.c $(MATHHDRS)
$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_set.c
+bn_mp_init_set_int.o: $(TOMMATH_DIR)/bn_mp_init_set_int.c $(MATHHDRS)
+ $(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_set_int.c
+
bn_mp_init_size.o:$(TOMMATH_DIR)/bn_mp_init_size.c $(MATHHDRS)
$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_init_size.c
@@ -1431,6 +1437,9 @@ bn_mp_rshd.o: $(TOMMATH_DIR)/bn_mp_rshd.c $(MATHHDRS)
bn_mp_set.o: $(TOMMATH_DIR)/bn_mp_set.c $(MATHHDRS)
$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_set.c
+bn_mp_set_int.o: $(TOMMATH_DIR)/bn_mp_set_int.c $(MATHHDRS)
+ $(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_set_int.c
+
bn_mp_shrink.o: $(TOMMATH_DIR)/bn_mp_shrink.c $(MATHHDRS)
$(CC) -c $(CC_SWITCHES) $(TOMMATH_DIR)/bn_mp_shrink.c
diff --git a/win/Makefile.in b/win/Makefile.in
index 3097849..74d006c 100644
--- a/win/Makefile.in
+++ b/win/Makefile.in
@@ -4,7 +4,7 @@
# "./configure", which is a configuration script generated by the "autoconf"
# program (constructs like "@foo@" will get replaced in the actual Makefile.
#
-# RCS: @(#) $Id: Makefile.in,v 1.185 2010/11/04 21:48:23 nijtmans Exp $
+# RCS: @(#) $Id: Makefile.in,v 1.186 2010/11/28 23:20:11 kennykb Exp $
VERSION = @TCL_VERSION@
@@ -318,6 +318,7 @@ TOMMATH_OBJS = \
bn_mp_init_copy.${OBJEXT} \
bn_mp_init_multi.${OBJEXT} \
bn_mp_init_set.${OBJEXT} \
+ bn_mp_init_set_int.${OBJEXT} \
bn_mp_init_size.${OBJEXT} \
bn_mp_karatsuba_mul.${OBJEXT} \
bn_mp_karatsuba_sqr.$(OBJEXT) \
@@ -335,6 +336,7 @@ TOMMATH_OBJS = \
bn_mp_read_radix.${OBJEXT} \
bn_mp_rshd.${OBJEXT} \
bn_mp_set.${OBJEXT} \
+ bn_mp_set_int.${OBJEXT} \
bn_mp_shrink.${OBJEXT} \
bn_mp_sqr.${OBJEXT} \
bn_mp_sqrt.${OBJEXT} \
diff --git a/win/makefile.vc b/win/makefile.vc
index 9cc1a71..988d823 100644
--- a/win/makefile.vc
+++ b/win/makefile.vc
@@ -13,7 +13,7 @@
# Copyright (c) 2003-2008 Pat Thoyts.
#
#------------------------------------------------------------------------------
-# RCS: @(#) $Id: makefile.vc,v 1.216 2010/11/04 21:48:23 nijtmans Exp $
+# RCS: @(#) $Id: makefile.vc,v 1.217 2010/11/28 23:20:11 kennykb Exp $
#------------------------------------------------------------------------------
# Check to see we are configured to build with MSVC (MSDEVDIR or MSVCDIR)
@@ -369,6 +369,7 @@ TOMMATHOBJS = \
$(TMP_DIR)\bn_mp_init_copy.obj \
$(TMP_DIR)\bn_mp_init_multi.obj \
$(TMP_DIR)\bn_mp_init_set.obj \
+ $(TMP_DIR)\bn_mp_init_set_int.obj \
$(TMP_DIR)\bn_mp_init_size.obj \
$(TMP_DIR)\bn_mp_karatsuba_mul.obj \
$(TMP_DIR)\bn_mp_karatsuba_sqr.obj \
@@ -386,6 +387,7 @@ TOMMATHOBJS = \
$(TMP_DIR)\bn_mp_read_radix.obj \
$(TMP_DIR)\bn_mp_rshd.obj \
$(TMP_DIR)\bn_mp_set.obj \
+ $(TMP_DIR)\bn_mp_set_int.obj \
$(TMP_DIR)\bn_mp_shrink.obj \
$(TMP_DIR)\bn_mp_sqr.obj \
$(TMP_DIR)\bn_mp_sqrt.obj \
generated by cgit v0.12 at 2024-06-28 15:20:33 (GMT)