From a70fba24c70c293375aabb150e77a53043eb67b7 Mon Sep 17 00:00:00 2001 From: Kevin B Kenny Date: Sun, 28 Nov 2010 23:20:10 +0000 Subject: 2010-11-29 Kevin B. Kenny * 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. --- ChangeLog | 28 + generic/tclInt.decls | 7 +- generic/tclInt.h | 29 +- generic/tclIntDecls.h | 8 +- generic/tclStrToD.c | 2906 +++++++++++++++++++++++++++++++++++++++------ generic/tclStubInit.c | 5 +- generic/tclTest.c | 105 +- generic/tclTomMath.decls | 8 +- generic/tclTomMathDecls.h | 14 +- generic/tclUtil.c | 167 ++- tests/util.test | 803 ++++++++++++- unix/Makefile.in | 15 +- win/Makefile.in | 4 +- win/makefile.vc | 4 +- 14 files changed, 3626 insertions(+), 477 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 + + * 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 * 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,435 +1927,2383 @@ RefineApproximation( } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * - * TclDoubleDigits -- - * - * Converts a double to a string of digits. + * MultPow5 -- * - * 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. + * Multiply a bignum by a power of 5. * * Side effects: - * Stores the digits in the given buffer and sets 'signum' according to - * the sign of the number. + * Stores base*5**n in result * - *---------------------------------------------------------------------- - + *----------------------------------------------------------------------------- */ -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. */ +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 */ { - 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; - - /* - * Split the number into absolute value and signum. - */ - - v = AbsoluteValue(v, signum); - - /* - * Handle zero specially. - */ - - if (v == 0.0) { - *buffer++ = '0'; - *buffer++ = '\0'; - return 1; + mp_int* p = base; + int n13 = n / 13; + int r = n % 13; + if (r != 0) { + mp_mul_d(p, dpow5[r], result); + p = result; } - - /* - * 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. - */ - - smallestSig = GetIntegerTimesPower(v, &r, &e); - - lowOK = highOK = (mp_iseven(&r)); - - /* - * 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. - */ - - if (e >= 0) { - int bits = e * log2FLT_RADIX; - - if (!smallestSig) { - /* - * Normal case, m+ and m- are both FLT_RADIX**e - */ - - rfac2 = bits + 1; - sfac2 = 1; - mplusfac2 = bits; - mminusfac2 = bits; - } else { - /* - * 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. - */ - - rfac2 = bits + log2FLT_RADIX + 1; - sfac2 = 1 + log2FLT_RADIX; - mplusfac2 = bits + log2FLT_RADIX; - mminusfac2 = bits; - } - } else { - /* - * v has digits after the binary point - */ - - 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. - */ - - rfac2 = 1; - sfac2 = 1 - e * log2FLT_RADIX; - mplusfac2 = 0; - mminusfac2 = 0; - } else { - /* - * 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. - */ - - rfac2 = 1 + log2FLT_RADIX; - sfac2 = 1 + log2FLT_RADIX * (1 - e); - mplusfac2 = FLT_RADIX; - mminusfac2 = 0; + 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 -- + * + * Shifts a number rightward until it is odd (that is, until the + * least significant bit is nonzero. + * + * Results: + * Returns the number of bit positions by which the number was shifted. + * + * Side effects: + * Shifts the number in place; *wPtr is replaced by the shifted number. + * + *----------------------------------------------------------------------------- + */ - /* - * Estimate the highest power of ten that will be needed to hold the - * result. - */ +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. + * + *---------------------------------------------------------------------------- + */ - k = (int) ceil(log(v) / log(10.)); - if (k >= 0) { - sfac2 += k; - sfac5 = k; +static int +RequiredPrecision(Tcl_WideUInt w) + /* Number to interrogate */ +{ + int rv; + unsigned long wi; + if (w & ((Tcl_WideUInt) 0xffffffff << 32)) { + wi = w >> 32; rv = 32; } else { - rfac2 -= k; - mplusfac2 -= k; - mminusfac2 -= k; - rfac5 = -k; + wi = w; rv = 0; } - - /* - * Scale r, s, mplus, mminus by the appropriate powers of 2 and 5. - */ - - mp_init_set(&mplus, 1); - for (i=0 ; i<=8 ; ++i) { - if (rfac5 & (1 << i)) { - mp_mul(&mplus, pow5+i, &mplus); - } + if (wi & 0xffff0000) { + wi >>= 16; rv += 16; } - 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); - } + if (wi & 0xff00) { + wi >>= 8; rv += 8; + } + if (wi & 0xf0) { + wi >>= 4; rv += 4; } - mp_mul_2d(&s, sfac2, &s); + 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'. + * + *----------------------------------------------------------------------------- + */ - /* - * It is possible for k to be off by one because we used an inexact - * logarithm. - */ +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'. + * + *----------------------------------------------------------------------------- + */ - 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++; +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 { - 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--; + 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. + * + *----------------------------------------------------------------------------- + */ - /* - * 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. - */ - - for (;;) { - int tc1, tc2; +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. + * + *----------------------------------------------------------------------------- + */ - 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); - } - mp_add(&r, &mplus, &temp); - tc2 = mp_cmp_mag(&temp, &s); - if (highOK) { - tc2 = (tc2 >= 0); - } else { - tc2 = (tc2 > 0); - } - if (!tc1) { - if (!tc2) { - *buffer++ = '0' + i; - } else { - c = (char) (i + '1'); - 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'); - } - } - break; - } - }; - *buffer++ = c; - *buffer++ = '\0'; +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; /* - * Free memory, and return. + * 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)) */ - mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL); + 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; } /* - *---------------------------------------------------------------------- - * - * AbsoluteValue -- + *----------------------------------------------------------------------------- * - * Splits a 'double' into its absolute value and sign. + * BetterLog10 -- * - * Results: - * Returns the absolute value. + * Improves the result of ApproximateLog10 for numbers in the range + * 1 .. 10**(TEN_PMAX)-1 * * Side effects: - * Stores the signum in '*signum'. + * Sets k_check to 0 if the new result is known to be exact, and to + * 1 if it may still be one too high. * - *---------------------------------------------------------------------- + * Results: + * Returns the improved approximation to log10(d) + * + *----------------------------------------------------------------------------- */ -static double -AbsoluteValue( - double v, /* Number to split */ - int *signum) /* (Output) Sign of the number 1=-, 0=+ */ +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 */ { - /* - * 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.) + /* + * Performance hack. If k is in the range 0..TEN_PMAX, then we can + * use a powers-of-ten table to check it. */ - -#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); - } - if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { - *signum = 1; - bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); - if (n770_fp) { - bitwhack.iv = Nokia770Twiddle(bitwhack.iv); + if (k >= 0 && k <= TEN_PMAX) { + if (d < tens[k]) { + k--; } - v = bitwhack.dv; + *k_check = 0; } else { - *signum = 0; + *k_check = 1; } -#endif - return v; + return k; } -/* - *---------------------------------------------------------------------- +/* + *----------------------------------------------------------------------------- * - * GetIntegerTimesPower -- + * ComputeScale -- * - * Converts a floating point number to an exact integer times a power of - * the floating point radix. + * Prepares to format a floating-point number as decimal. * - * Results: - * Returns 1 if it converted the smallest significand, 0 otherwise. + * Parameters: + * floor(log10*x) is k (or possibly k-1). floor(log2(x) is i. + * The significand of x requires bbits bits to represent. * - * Side effects: - * Initializes the integer value (does not just assign it), and stores - * the exponent. + * Results: + * Determines integers b2, b5, s2, s5 so that sig*2**b2*5**b5/2**s2*2**s5 + * exactly represents the value of the x/10**k. This value will lie + * in the range [1 .. 10), and allows for computing successive digits + * by multiplying sig%10 by 10. * - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- */ -static 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*/ +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 */ { - double a, f; - int e, i, n; - /* - * Develop f and e such that v = f * FLT_RADIX**e, with - * 1.0/FLT_RADIX <= f < 1. + /* + * Scale numerator and denominator powers of 2 so that the + * input binary number is the ratio of integers */ - - 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); - } - e = (e - n) / log2FLT_RADIX; -#endif - if (f == 1.0) { - f = 1.0 / FLT_RADIX; - e += 1; + if (be <= 0) { + *b2 = 0; + *s2 = -be; + } else { + *b2 = be; + *s2 = 0; } - /* - * If the original number was denormalized, adjust e and f to be denormal - * as well. + /* + * Scale numerator and denominator so that the output decimal number + * is the ratio of integers */ - - 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; + if (k >= 0) { + *b5 = 0; + *s5 = k; + *s2 += k; } else { - n = mantDIGIT; + *b2 -= k; + *b5 = -k; + *s5 = 0; } - - /* - * 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); } /* - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- * - * TclInitDoubleConversion -- + * SetPrecisionLimits -- * - * Initializes constants that are needed for conversions to and from - * 'double' + * Determines how many digits of significance should be computed + * (and, hence, how much memory need be allocated) for formatting a + * floating point number. * - * Results: - * None. + * 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: - * The log base 2 of the floating point radix, the number of bits in a - * double mantissa, and a table of the powers of five and ten are - * computed and stored. + * 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. * - *---------------------------------------------------------------------- + *----------------------------------------------------------------------------- */ -void +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 */ + + /* + * Bring d into the range [1 .. 10) + */ + ieps = AdjustRange(&d, k); + + /* + * 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; + } + + /* + * Compute estimated roundoff error + */ + eps.d = ieps * d + 7.; + eps.w.word0 -= (FP_PRECISION-1) << EXP_SHIFT; + + /* + * Handle the peculiar case where the result has no significant + * digits. + */ + retval = ckalloc(len + 1); + if (ilim == 0) { + d -= 5.; + if (d > eps.d) { + *retval = '1'; + *decpt = k; + return retval; + } else if (d < -eps.d) { + *decpt = k; + return retval; + } else { + ckfree(retval); + return NULL; + } + } + + /* Format the digit string */ + + if (flags & TCL_DD_SHORTEN_FLAG) { + end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt); + } else { + end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt); + } + if (end == NULL) { + ckfree(retval); + return NULL; + } + *end = '\0'; + if (endPtr != NULL) { + *endPtr = end; + } + return retval; +} + +/* + *----------------------------------------------------------------------------- + * + * 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. + * + *----------------------------------------------------------------------------- + */ + +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)) { + /* + * Make sure we shouldn't be rounding *up* instead, + * in case the next number above is closer + */ + if (2 * b > S + || (2 * b == S + && (digit & 1) != 0)) { + ++digit; + if (digit == 10) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + } + + /* Stash the current digit */ + + *s++ = '0' + digit; + break; + } + + /* + * Does one plus the current digit put us within roundoff of the + * number? + */ + if (b > S - mminus + || (b == S - mminus + && convType != TCL_DD_STEELE0 + && (dPtr->w.word1 & 1) == 0)) { + if (digit == 9) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + ++digit; + *s++ = '0' + digit; + break; + } + + /* + * Have we converted all the requested digits? + */ + *s++ = '0' + digit; + if (i == ilim) { + if (2*b > S + || (2*b == S && (digit & 1) != 0)) { + s = BumpUp(s, retval, &k); + } + break; + } + + /* Advance to the next digit */ + + b = 10 * b; + mplus = 10 * mplus; + mminus = 10 * mminus; + ++i; + } + + /* + * Endgame - store the location of the decimal point and the end of the + * string. + */ + *s = '\0'; + *decpt = k; + if (endPtr) { + *endPtr = s; + } + return retval; +} + +/* + *----------------------------------------------------------------------------- + * + * 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 */ + + i = 0; + for (;;) { + if (b.used <= sd) { + digit = 0; + } else { + digit = b.dp[sd]; + if (b.used > sd+1 || digit >= 10) { + Tcl_Panic("wrong digit!"); + } + --b.used; mp_clamp(&b); + } + + /* + * Does the current digit put us on the low side of the exact value + * but within within roundoff of being exact? + */ + + r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus); + if (r1 == MP_LT + || (r1 == MP_EQ + && convType != TCL_DD_STEELE0 + && (dPtr->w.word1 & 1) == 0)) { + /* + * Make sure we shouldn't be rounding *up* instead, + * in case the next number above is closer + */ + if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { + ++digit; + if (digit == 10) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + } + + /* Stash the last digit */ + + *s++ = '0' + digit; + break; + } + + /* + * Does one plus the current digit put us within roundoff of the + * number? + */ + + if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd, + convType, dPtr->w.word1 & 1, + &temp)) { + if (digit == 9) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + ++digit; + *s++ = '0' + digit; + break; + } + + /* + * Have we converted all the requested digits? + */ + *s++ = '0' + digit; + if (i == ilim) { + if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { + s = BumpUp(s, retval, &k); + } + break; + } + + /* Advance to the next digit */ + + mp_mul_d(&b, 10, &b); + mp_mul_d(&mminus, 10, &mminus); + if (m2plus > m2minus) { + mp_mul_2d(&mminus, m2plus-m2minus, &mplus); + } + ++i; + } + + /* + * Endgame - store the location of the decimal point and the end of the + * string. + */ + if (m2plus > m2minus) { + mp_clear(&mplus); + } + mp_clear_multi(&b, &mminus, &temp, NULL); + *s = '\0'; + *decpt = k; + if (endPtr) { + *endPtr = s; + } + return retval; +} + +/* + *----------------------------------------------------------------------------- + * + * 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 + */ + + TclBNInitBignumFromWideUInt(&b, bw); + MulPow5(&b, b5, &b); + mp_mul_2d(&b, b2, &b); + + /* Adjust if the logarithm was guessed wrong */ + + if (b.used <= sd) { + mp_mul_d(&b, 10, &b); + ilim = ilim1; + --k; + } + mp_init(&temp); + + /* + * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT) + * by mp_digit extraction + */ + + i = 1; + for (;;) { + if (b.used <= sd) { + digit = 0; + } else { + digit = b.dp[sd]; + if (b.used > sd+1 || digit >= 10) { + Tcl_Panic("wrong digit!"); + } + --b.used; mp_clamp(&b); + } + + /* + * Have we converted all the requested digits? + */ + *s++ = '0' + digit; + if (i == ilim) { + if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) { + s = BumpUp(s, retval, &k); + } + break; + } + + /* Advance to the next digit */ + + mp_mul_d(&b, 10, &b); + ++i; + } + + /* + * Endgame - store the location of the decimal point and the end of the + * string. + */ + mp_clear_multi(&b, &temp, NULL); + *s = '\0'; + *decpt = k; + if (endPtr) { + *endPtr = s; + } + return retval; +} + +/* + *----------------------------------------------------------------------------- + * + * 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. + * + *----------------------------------------------------------------------------- + */ + +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; + + /* + * b = bw * 2**b2 * 5**b5 + * S = 2**s2 * 5*s5 + */ + + 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); + 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!"); + } + digit = dig.dp[0]; + + /* + * Does the current digit leave us with a remainder small enough to + * round to it? + */ + + r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus); + if (r1 == MP_LT + || (r1 == MP_EQ + && convType != TCL_DD_STEELE0 + && (dPtr->w.word1 & 1) == 0)) { + mp_mul_2d(&b, 1, &b); + if (ShouldBankerRoundUp(&b, &S, digit&1)) { + ++digit; + if (digit == 10) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + } + *s++ = '0' + digit; + break; + } + + /* + * Does the current digit leave us with a remainder large enough + * to commit to rounding up to the next higher digit? + */ + + if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType, + dPtr->w.word1 & 1, &temp)) { + ++digit; + if (digit == 10) { + *s++ = '9'; + s = BumpUp(s, retval, &k); + break; + } + *s++ = '0' + digit; + break; + } + + /* Have we converted all the requested digits? */ + + *s++ = '0' + digit; + if (i == ilim) { + mp_mul_2d(&b, 1, &b); + if (ShouldBankerRoundUp(&b, &S, digit&1)) { + s = BumpUp(s, retval, &k); + } + break; + } + + /* Advance to the next digit */ + + if (s5 > 0) { + + /* Can possibly shorten the denominator */ + mp_mul_2d(&b, 1, &b); + mp_mul_2d(&mminus, 1, &mminus); + if (m2plus > m2minus) { + mp_mul_2d(&mplus, 1, &mplus); + } + mp_div_d(&S, 5, &S, NULL); + --s5; + /* + * 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; + +} + +/* + *----------------------------------------------------------------------------- + * + * StrictBignumConversion -- + * + * Convert a floating point number to a fixed-length digit string + * using the multiprecision method. + * + * Results: + * Returns the string of digits. + * + * Side effects: + * Stores the position of the decimal point in *decpt. + * Stores a pointer to the end of the number in *endPtr. + * + *----------------------------------------------------------------------------- + */ + +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; + + /* + * b = bw * 2**b2 * 5**b5 + * S = 2**s2 * 5*s5 + */ + + 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; + } + 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); + } + } else { + + for (;;) { + + /* Shift by a group of digits. */ + + g = ilim - i; + if (g > DIGIT_GROUP) { + g = DIGIT_GROUP; + } + if (s5 >= g) { + mp_div_d(&S, dpow5[g], &S, NULL); + s5 -= g; + } else if (s5 > 0) { + mp_div_d(&S, dpow5[s5], &S, NULL); + mp_mul_d(&b, dpow5[g - s5], &b); + s5 = 0; + } else { + mp_mul_d(&b, dpow5[g], &b); + } + mp_mul_2d(&b, g, &b); + + /* + * As with the shortening bignum conversion, it's possible at + * this point that we will have reduced the denominator to + * less than 2**64/10, at which point it would be possible to + * fall back to to int64 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; + } + } + } + /* + * 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; + +} + +/* + *----------------------------------------------------------------------------- + * + * TclDoubleDigits -- + * + * Core of Tcl's conversion of double-precision floating point numbers + * to decimal. + * + * Results: + * Returns a newly-allocated string of digits. + * + * Side effects: + * Sets *decpt to the index of the character in the string before the + * place that the decimal point should go. If 'endPtr' is not NULL, + * sets endPtr to point to the terminating '\0' byte of the string. + * Sets *sign to 1 if a minus sign should be printed with the number, + * or 0 if a plus sign (or no sign) should appear. + * + * This function is a service routine that produces the string of digits + * for floating-point-to-decimal conversion. It can do a number of things + * according to the 'flags' argument. Valid values for 'flags' include: + * TCL_DD_SHORTEST - This is the default for floating point conversion + * if ::tcl_precision is 0. It constructs the shortest string + * of digits that will reconvert to the given number when scanned. + * For floating point numbers that are exactly between two + * decimal numbers, it resolves using the 'round to even' rule. + * With this value, the 'ndigits' parameter is ignored. + * TCL_DD_STEELE - This value is not recommended and may be removed + * in the future. It follows the conversion algorithm outlined + * in "How to Print Floating-Point Numbers Accurately" by + * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, + * pp. 112-126]. This rule has the effect of rendering 1e23 + * as 9.9999999999999999e22 - which is a 'better' approximation + * in the sense that it will reconvert correctly even if + * a subsequent input conversion is 'round up' or 'round down' + * rather than 'round to nearest', but is surprising otherwise. + * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e + * format conversion (or for default floating->string if + * tcl_precision is not 0). It constructs a string of at most + * 'ndigits' digits, choosing the one that is closest to the + * given number (and resolving ties with 'round to even'). + * It is allowed to return fewer than 'ndigits' if the number + * converts exactly; if the TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG + * is supplied instead, it 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. + * + *----------------------------------------------------------------------------- + */ +char* +TclDoubleDigits(double dv, /* Number to convert */ + int ndigits, /* Number of digits requested */ + int flags, /* Conversion flags */ + int* decpt, /* OUTPUT: Position of the decimal point */ + int* sign, /* OUTPUT: 1 if the result is negative */ + char** endPtr) /* OUTPUT: If not NULL, receives a pointer + * to one character beyond the end + * of the returned string */ +{ + int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK); + /* Type of conversion being performed + * TCL_DD_SHORTEST0 + * TCL_DD_STEELE0 + * TCL_DD_E_FORMAT + * 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; + + /* Put the input number into a union for bit-whacking */ + + d.d = dv; + + /* + * Handle the cases of negative numbers (by taking the absolute value: + * this includes -Inf and -NaN!), infinity, Not a Number, and zero. + */ + + TakeAbsoluteValue(&d, sign); + if ((d.w.word0 & EXP_MASK) == EXP_MASK) { + return FormatInfAndNaN(&d, decpt, endPtr); + } + if (d.d == 0.0) { + return FormatZero(decpt, endPtr); + } + + /* + * Unpack the floating point into a wide integer and an exponent. + * Determine the number of bits that the big integer requires, and + * compute a quick approximation (which may be one too high) of + * ceil(log10(d.d)). + */ + denorm = ((d.w.word0 & EXP_MASK) == 0); + DoubleToExpAndSig(d.d, &bw, &be, &bbits); + k = ApproximateLog10(bw, be, bbits); + k = BetterLog10(d.d, k, &k_check); + + /* At this point, we have: + * d is the number to convert. + * bw are significand and exponent: d == bw*2**be, + * bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits + * k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0 + * if we know that k is exactly ceil(log10(d)) and 1 if we need to + * check. + * We want a rational number + * r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5), + * with b2, b5, s2, s5 >= 0. Note that the most significant decimal + * digit is floor(r) and that successive digits can be obtained + * by setting r <- 10*floor(r) (or b <= 10 * (b % S)). + * Find appropriate b2, b5, s2, s5. + */ + + ComputeScale(be, k, &b2, &b5, &s2, &s5); + + /* + * Correct an incorrect caller-supplied 'ndigits'. + * Also determine: + * i = The maximum number of decimal digits that will be returned in the + * formatted string. This is k + 1 + ndigits for F format, 18 for + * shortest and Steele, and ndigits for E format. + * ilim = The number of significant digits to convert if + * k has been guessed correctly. This is -1 for shortest and Steele + * (which stop when all significance has been lost), 'ndigits' + * for E format, and 'k + 1 + ndigits' for F format. + * ilim1 = The minimum number of significant digits to convert if + * k has been guessed 1 too high. This, too, is -1 for shortest + * and Steele, and 'ndigits' for E format, but it's 'ndigits-1' + * for F format. + */ + + SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1); + + /* + * Try to do low-precision conversion in floating point rather + * than resorting to expensive multiprecision arithmetic + */ + if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) { + if ((retval = QuickConversion(d.d, k, k_check, flags, + i, ilim, ilim1, + decpt, endPtr)) != NULL) { + return retval; + } + } + + /* + * For shortening conversions, determine the upper and lower bounds + * for the remainder at which we can stop. + * m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the + * high side, and + * m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the + * low side. + * We may need to increase s2 to put m2plus, m2minus, b2 over a + * common denominator. + */ + + 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); + } + } +} + + +/* + *---------------------------------------------------------------------- + * + * TclInitDoubleConversion -- + * + * Initializes constants that are needed for conversions to and from + * 'double' + * + * Results: + * None. + * + * Side effects: + * The log base 2 of the floating point radix, the number of bits in a + * double mantissa, and a table of the powers of five and ten are + * computed and stored. + * + *---------------------------------------------------------------------- + */ + +void TclInitDoubleConversion(void) { int i; @@ -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 + #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 \ -- cgit v0.12