/* *---------------------------------------------------------------------- * * tclStrToD.c -- * * This file contains a TclStrToD procedure that handles conversion of * string to double, with correct rounding even where extended precision * is needed to achieve that. It also contains a TclDoubleDigits * procedure that handles conversion of double to string (at least the * significand), and several utility functions for interconverting * 'double' and the integer types. * * Copyright (c) 2005 by Kevin B. Kenny. All rights reserved. * * See the file "license.terms" for information on usage and redistribution * of this file, and for a DISCLAIMER OF ALL WARRANTIES. * * RCS: @(#) $Id: tclStrToD.c,v 1.8 2005/08/24 15:15:45 kennykb Exp $ * *---------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include /* * The stuff below is a bit of a hack so that this file can be used in * environments that include no UNIX, i.e. no errno: just arrange to use the * errno from tclExecute.c here. */ #ifdef TCL_GENERIC_ONLY # define NO_ERRNO_H #endif #ifdef NO_ERRNO_H extern int errno; /* Use errno from tclExecute.c. */ # define ERANGE 34 #endif #if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024) # define IEEE_FLOATING_POINT #endif /* * gcc on x86 needs access to rounding controls. It is tempting to include * fpu_control.h, but that file exists only on Linux; it is missing on Cygwin * and MinGW. Most gcc-isms and ix86-isms are factored out here. */ #if defined(__GNUC__) && defined(__i386) typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__))); # define _FPU_GETCW(cw) __asm__ ("fnstcw %0" : "=m" (*&cw)) # define _FPU_SETCW(cw) __asm__ ("fldcw %0" : : "m" (*&cw)) # define FPU_IEEE_ROUNDING 0x027f # define ADJUST_FPU_CONTROL_WORD #endif /* * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN. * Everyone else uses 7ff8000000000000. (Why, HP, why?) */ #ifdef __hppa # define NAN_START 0x7ff4 # define NAN_MASK (((Tcl_WideUInt) 1) << 50) #else # define NAN_START 0x7ff8 # define NAN_MASK (((Tcl_WideUInt) 1) << 51) #endif /* * There now follows a lot of static variables that are shared across all * threads but which are not guarded by mutexes. This is OK, because they are * only ever assigned _ONCE_ during Tcl's library initialization sequence. */ static const double pow_10_2_n[] = { /* Inexact higher powers of ten */ 1.0, 100.0, 10000.0, 1.0e+8, 1.0e+16, 1.0e+32, 1.0e+64, 1.0e+128, 1.0e+256 }; #define MAXPOW 22 /* Num of exactly representable powers of 10 */ static double pow10[MAXPOW+1]; /* The powers of ten that can be represented * exactly as IEEE754 doubles. */ static int mmaxpow; /* Largest power of ten that can be * represented exactly in a 'double'. */ static int log2FLT_RADIX; /* Logarithm of the floating point radix. */ static int mantBits; /* Number of bits in a double's significand. */ static mp_int pow5[9]; /* Table of powers of 5**(2**n), up to * 5**256. */ static double tiny; /* The smallest representable double. */ static int maxDigits; /* The maximum number of digits to the left of * the decimal point of a double. */ static int minDigits; /* The maximum number of digits to the right * of the decimal point in a double. */ static int mantDIGIT; /* Number of mp_digit's needed to hold the * significand of a double. */ /* Static functions defined in this file */ static double RefineResult(double approx, CONST char *start, int nDigits, long exponent); static double ParseNaN(int signum, CONST char **end); static double SafeLdExp(double fraction, int exponent); /* *---------------------------------------------------------------------- * * TclStrToD -- * * Scans a double from a string. * * Results: * Returns the scanned number. In the case of underflow, returns an * appropriately signed zero; in the case of overflow, returns an * appropriately signed HUGE_VAL. * * Side effects: * Stores a pointer to the end of the scanned number in '*endPtr', if * endPtr is not NULL. If '*endPtr' is equal to 's' on return from this * function, it indicates that the input string could not be recognized * as a number. In the case of underflow or overflow, 'errno' is set to * ERANGE. * *------------------------------------------------------------------------ */ double TclStrToD(CONST char *s, /* String to scan. */ CONST char **endPtr) /* Pointer to the end of the scanned number. */ { const char *p = s; const char *startOfSignificand = NULL; /* Start of the significand in the string. */ int signum = 0; /* Sign of the significand. */ double exactSignificand = 0.0; /* Significand, represented exactly as a * floating-point number. */ int seenDigit = 0; /* Flag == 1 if a digit has been seen. */ int nSigDigs = 0; /* Number of significant digits presented. */ int nDigitsAfterDp = 0; /* Number of digits after the decimal point. */ int nTrailZero = 0; /* Number of trailing zeros in the * significand. */ long exponent = 0; /* Exponent. */ int seenDp = 0; /* Flag == 1 if decimal point has been seen. */ char c; /* One character extracted from the input. */ volatile double v; /* Scanned value; must be 'volatile double' on * gc-ix86 to force correct rounding to IEEE * double and not Intel double-extended. */ int machexp; /* Exponent of the machine rep of the scanned * value. */ int expt2; /* Exponent for computing first approximation * to the true value. */ int i, j; /* * With gcc on x86, the floating point rounding mode is double-extended. * This causes the result of double-precision calculations to be rounded * twice: once to the precision of double-extended and then again to the * precision of double. Double-rounding introduces gratuitous errors of * one ulp, so we need to change rounding mode to 53-bits. */ #ifdef ADJUST_FPU_CONTROL_WORD fpu_control_t roundTo53Bits = FPU_IEEE_ROUNDING; fpu_control_t oldRoundingMode; _FPU_GETCW(oldRoundingMode); _FPU_SETCW(roundTo53Bits); # define RestoreRoundingMode() _FPU_SETCW(oldRoundingMode) #else # define RestoreRoundingMode() (void) 0 /* Do nothing */ #endif /* * Discard leading whitespace from input. */ while (isspace(UCHAR(*p))) { ++p; } /* * Determine the sign of the significand. */ switch (*p) { case '-': signum = 1; /* FALLTHROUGH */ case '+': ++p; } /* * Discard leading zeroes from input. */ while (*p == '0') { seenDigit = 1; ++p; } /* * Scan digits from the significand. Simultaneously, keep track of the * number of digits after the decimal point. Maintain a pointer to the * start of the significand. Keep "exactSignificand" equal to the * conversion of the DBL_DIG most significant digits. */ for (;;) { c = *p; if (c == '.' && !seenDp) { seenDp = 1; ++p; } else if (isdigit(UCHAR(c))) { if (c == '0') { if (startOfSignificand != NULL) { ++nTrailZero; } } else { if (startOfSignificand == NULL) { startOfSignificand = p; } else if (nTrailZero) { if (nTrailZero + nSigDigs < DBL_DIG) { exactSignificand *= pow10[nTrailZero]; } else if (nSigDigs < DBL_DIG) { exactSignificand *= pow10[DBL_DIG - nSigDigs]; } nSigDigs += nTrailZero; } if (nSigDigs < DBL_DIG) { exactSignificand = 10. * exactSignificand + (c - '0'); } ++nSigDigs; nTrailZero = 0; } if (seenDp) { ++nDigitsAfterDp; } seenDigit = 1; ++p; } else { break; } } /* * At this point, we've scanned the significand, and p points to the * character beyond it. "startOfSignificand" is the first non-zero * character in the significand. "nSigDigs" is the number of significant * digits of the significand, not including any trailing zeroes. * "exactSignificand" is a floating point number that represents, without * loss of precision, the first min(DBL_DIG,n) digits of the significand. * "nDigitsAfterDp" is the number of digits after the decimal point, again * excluding trailing zeroes. * * Now scan 'E' format */ exponent = 0; if (seenDigit && (*p == 'e' || *p == 'E')) { const char* stringSave = p; ++p; c = *p; if (isdigit(UCHAR(c)) || c == '+' || c == '-') { errno = 0; exponent = strtol(p, (char**)&p, 10); if (errno == ERANGE) { if (exponent > 0) { v = HUGE_VAL; } else { v = 0.0; } *endPtr = p; goto returnValue; } } if (p == stringSave+1) { p = stringSave; exponent = 0; } } exponent += nTrailZero - nDigitsAfterDp; /* * If we come here with no significant digits, we might still be looking * at Inf or NaN. Go parse them. */ if (!seenDigit) { /* * Test for Inf or Infinity (in any case). */ if (c == 'I' || c == 'i') { if ((p[1] == 'N' || p[1] == 'n') && (p[2] == 'F' || p[2] == 'f')) { p += 3; if ((p[0] == 'I' || p[0] == 'i') && (p[1] == 'N' || p[1] == 'n') && (p[2] == 'I' || p[2] == 'i') && (p[3] == 'T' || p[3] == 't') && (p[4] == 'Y' || p[1] == 'y')) { p += 5; } errno = ERANGE; v = HUGE_VAL; if (endPtr != NULL) { *endPtr = p; } goto returnValue; } #ifdef IEEE_FLOATING_POINT /* * Only IEEE floating point supports NaN */ } else if ((c == 'N' || c == 'n') && (sizeof(Tcl_WideUInt) == sizeof(double))) { if ((p[1] == 'A' || p[1] == 'a') && (p[2] == 'N' || p[2] == 'n')) { p += 3; if (endPtr != NULL) { *endPtr = p; } /* * Restore FPU mode word. */ RestoreRoundingMode(); return ParseNaN(signum, endPtr); } #endif } goto error; } /* * We've successfully scanned; update the end-of-element pointer. */ if (endPtr != NULL) { *endPtr = p; } /* * Test for zero. */ if (nSigDigs == 0) { v = 0.0; goto returnValue; } /* * The easy cases are where we have an exact significand and the exponent * is small enough that we can compute the value with only one roundoff. * In addition to the cases where we can multiply or divide an * exact-integer significand by an exact-integer power of 10, there is * also David Gay's case where we can scale the significand by a power of * 10 (still keeping it exact) and then multiply by an exact power of 10. * The last case enables combinations like 83e25 that would otherwise * require high precision arithmetic. */ if (nSigDigs <= DBL_DIG) { if (exponent >= 0) { if (exponent <= mmaxpow) { v = exactSignificand * pow10[exponent]; goto returnValue; } else { int diff = DBL_DIG - nSigDigs; if (exponent - diff <= mmaxpow) { volatile double factor = exactSignificand * pow10[diff]; v = factor * pow10[exponent - diff]; goto returnValue; } } } else if (exponent >= -mmaxpow) { v = exactSignificand / pow10[-exponent]; goto returnValue; } } /* * We don't have one of the easy cases, so we can't compute the scanned * number exactly, and have to do it in multiple precision. Begin by * testing for obvious overflows and underflows. */ if (nSigDigs + exponent - 1 > maxDigits) { v = HUGE_VAL; errno = ERANGE; goto returnValue; } if (nSigDigs + exponent - 1 < minDigits) { errno = ERANGE; v = 0.; goto returnValue; } /* * Nothing exceeds the boundaries of the tables, at least. Compute an * approximate value for the number, with no possibility of overflow * because we manage the exponent separately. */ if (nSigDigs > DBL_DIG) { expt2 = exponent + nSigDigs - DBL_DIG; } else { expt2 = exponent; } v = frexp(exactSignificand, &machexp); if (expt2 > 0) { v = frexp(v * pow10[expt2 & 0xf], &j); machexp += j; for (i=4 ; i<9 ; ++i) { if (expt2 & (1 << i)) { v = frexp(v * pow_10_2_n[i], &j); machexp += j; } } } else { v = frexp(v / pow10[(-expt2) & 0xf], &j); machexp += j; for (i=4 ; i<9 ; ++i) { if ((-expt2) & (1 << i)) { v = frexp(v / pow_10_2_n[i], &j); machexp += j; } } } /* * A first approximation is that the result will be v * 2 ** machexp. v is * greater than or equal to 0.5 and less than 1. If machexp > * DBL_MAX_EXP*log2(FLT_RADIX), there is an overflow. Constrain the result * to the smallest representible number to avoid premature underflow. */ if (machexp > DBL_MAX_EXP * log2FLT_RADIX) { v = HUGE_VAL; errno = ERANGE; goto returnValue; } v = SafeLdExp(v, machexp); if (v < tiny) { v = tiny; } /* * We have a first approximation in v. Now we need to refine it. */ v = RefineResult(v, startOfSignificand, nSigDigs, exponent); /* * In a very few cases, a second iteration is needed. e.g., 457e-102 */ v = RefineResult(v, startOfSignificand, nSigDigs, exponent); /* * Handle underflow. */ returnValue: if (nSigDigs != 0 && v == 0.0) { errno = ERANGE; } /* * Return a number with correct sign. */ if (signum) { v = -v; } /* * Restore FPU mode word and return. */ RestoreRoundingMode(); return v; /* * Come here on an invalid input. */ error: if (endPtr != NULL) { *endPtr = s; } /* * Restore FPU mode word and return. */ RestoreRoundingMode(); return 0.0; } /* *---------------------------------------------------------------------- * * RefineResult -- * * Given a poor approximation to a floating point number, returns a * better one. (The better approximation is correct to within 1 ulp, and * is entirely correct if the poor approximation is correct to 1 ulp.) * * Results: * Returns the improved result. * *---------------------------------------------------------------------- */ static double RefineResult(double approxResult, /* Approximate result of conversion. */ CONST char* sigStart, /* Pointer to start of significand in input * string. */ int nSigDigs, /* Number of significant digits. */ long exponent) /* Power of ten to multiply by significand. */ { int M2, M5; /* Powers of 2 and of 5 needed to put the * decimal and binary numbers over a common * denominator. */ double significand; /* Sigificand of the binary number. */ int binExponent; /* Exponent of the binary number. */ int msb; /* Most significant bit position of an * intermediate result. */ int nDigits; /* Number of mp_digit's in an intermediate * result. */ mp_int twoMv; /* Approx binary value expressed as an exact * integer scaled by the multiplier 2M. */ mp_int twoMd; /* Exact decimal value expressed as an exact * integer scaled by the multiplier 2M. */ int scale; /* Scale factor for M. */ int multiplier; /* Power of two to scale M. */ double num, den; /* Numerator and denominator of the correction * term. */ double quot; /* Correction term. */ double minincr; /* Lower bound on the absolute value of the * correction term. */ int i; const char* p; /* * The first approximation is always low. If we find that it's HUGE_VAL, * we're done. */ if (approxResult == HUGE_VAL) { return approxResult; } /* * Find a common denominator for the decimal and binary fractions. The * common denominator will be 2**M2 + 5**M5. */ significand = frexp(approxResult, &binExponent); i = mantBits - binExponent; if (i < 0) { M2 = 0; } else { M2 = i; } if (exponent > 0) { M5 = 0; } else { M5 = -exponent; if ((M5-1) > M2) { M2 = M5-1; } } /* * The floating point number is significand*2**binExponent. The 2**-1 bit * of the significand (the most significant) corresponds to the * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold * that quantity, then convert the significand to a large integer, scaled * appropriately. Then multiply by the appropriate power of 5. */ msb = binExponent + M2; /* 1008 */ nDigits = msb / DIGIT_BIT + 1; mp_init_size(&twoMv, nDigits); i = (msb % DIGIT_BIT + 1); twoMv.used = nDigits; significand *= SafeLdExp(1.0, i); while (--nDigits >= 0) { twoMv.dp[nDigits] = (mp_digit) significand; significand -= (mp_digit) significand; significand = SafeLdExp(significand, DIGIT_BIT); } for (i=0 ; i<=8 ; ++i) { if (M5 & (1 << i)) { mp_mul(&twoMv, pow5+i, &twoMv); } } /* * Collect the decimal significand as a high precision integer. The least * significant bit corresponds to bit M2+exponent+1 so it will need to be * shifted left by that many bits after being multiplied by * 5**(M5+exponent). */ mp_init(&twoMd); mp_zero(&twoMd); i = nSigDigs; for (p=sigStart ;; ++p) { char c = *p; if (isdigit(UCHAR(c))) { mp_mul_d(&twoMd, (unsigned) 10, &twoMd); mp_add_d(&twoMd, (unsigned) (c - '0'), &twoMd); --i; if (i == 0) { break; } } } for (i=0 ; i<=8 ; ++i) { if ((M5+exponent) & (1 << i)) { mp_mul(&twoMd, pow5+i, &twoMd); } } mp_mul_2d(&twoMd, M2+exponent+1, &twoMd); mp_sub(&twoMd, &twoMv, &twoMd); /* * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction * term. Because 2M may well overflow a double, we need to scale the * denominator by a factor of 2**binExponent-mantBits */ scale = binExponent - mantBits - 1; mp_set(&twoMv, 1); for (i=0 ; i<=8 ; ++i) { if (M5 & (1 << i)) { mp_mul(&twoMv, pow5+i, &twoMv); } } multiplier = M2 + scale + 1; if (multiplier > 0) { mp_mul_2d(&twoMv, multiplier, &twoMv); } else if (multiplier < 0) { mp_div_2d(&twoMv, -multiplier, &twoMv, NULL); } /* * If the result is less than unity, the error is less than 1/2 unit in * the last place, so there's no correction to make. */ if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) { mp_clear(&twoMd); mp_clear(&twoMv); return approxResult; } /* * Convert the numerator and denominator of the corrector term accurately * to floating point numbers. */ num = TclBignumToDouble(&twoMd); den = TclBignumToDouble(&twoMv); quot = SafeLdExp(num/den, scale); minincr = SafeLdExp(1.0, binExponent - mantBits); if (quot<0. && quot>-minincr) { quot = -minincr; } else if (quot>0. && quot= 'A' && c <= 'F') { c -= 'A' + 10; } else if (c >= 'a' && c <= 'f') { c -= 'a' + 10; } else { theNaN.iv = (((Tcl_WideUInt) NAN_START) << 48) | (((Tcl_WideUInt) signum) << 63); return theNaN.dv; } theNaN.iv = (theNaN.iv << 4) | c; } } /* * Mask the hex number down to the least significant 51 bits. */ theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1; if (signum) { theNaN.iv |= ((Tcl_WideUInt) 0xfff8) << 48; } else { theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48; } *endPtr = p; return theNaN.dv; } /* *---------------------------------------------------------------------- * * TclDoubleDigits -- * * Converts a double to a string of digits. * * 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. * * Side effects: * Stores the digits in the given buffer and sets 'signum' according to * the sign of the number. * *---------------------------------------------------------------------- */ int TclDoubleDigits(char * strPtr, /* 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. */ { double f; /* Significand of v. */ 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 */ 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; double a; char c; int i, k, n; /* * 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.) */ #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 (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { *signum = 1; bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63); v = bitwhack.dv; } else { *signum = 0; } #endif /* * Handle zero specially. */ if (v == 0.0) { *strPtr++ = '0'; *strPtr++ = '\0'; return 1; } /* * Develop f and e such that v = f * FLT_RADIX**e, with * 1.0/FLT_RADIX <= f < 1. */ f = frexp(v, &e); n = e % log2FLT_RADIX; if (n > 0) { n -= log2FLT_RADIX; e += 1; } f *= ldexp(1.0, n); e = (e - n) / log2FLT_RADIX; if (f == 1.0) { f = 1.0 / FLT_RADIX; e += 1; } /* * If the original number was denormalized, adjust e and f to be denormal * as well. */ if (e < DBL_MIN_EXP) { n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX; f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX); e = DBL_MIN_EXP; n = (n + DIGIT_BIT - 1) / DIGIT_BIT; } else { n = mantDIGIT; } /* * 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(&r, n); r.used = n; r.sign = MP_ZPOS; i = (mantBits % DIGIT_BIT); if (i == 0) { i = DIGIT_BIT; } while (n > 0) { a *= ldexp(1.0, i); i = DIGIT_BIT; r.dp[--n] = (mp_digit) a; a -= (mp_digit) a; } e -= DBL_MANT_DIG; 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 (f != 1.0/FLT_RADIX) { /* * 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 || f != 1.0/FLT_RADIX) { /* * 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; } } /* * Estimate the highest power of ten that will be needed to hold the * result. */ k = (int) ceil(log(v) / log(10.)); if (k >= 0) { sfac2 += k; sfac5 = k; } else { rfac2 -= k; mplusfac2 -= k; mminusfac2 -= k; rfac5 = -k; } /* * 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); } } 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); } } mp_mul_2d(&s, sfac2, &s); /* * It is possible for k to be off by one because we used an inexact * logarithm. */ mp_init(&temp); mp_add(&r, &mplus, &temp); i = mp_cmp_mag(&temp, &s); if (i>0 || (highOK && i==0)) { mp_mul_d(&s, 10, &s); ++k; } else { mp_mul_d(&temp, 10, &temp); i = mp_cmp_mag(&temp, &s); if (i<0 || (highOK && i==0)) { mp_mul_d(&r, 10, &r); mp_mul_d(&mplus, 10, &mplus); mp_mul_d(&mminus, 10, &mminus); --k; } } /* * 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; 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) { *strPtr++ = '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; } }; *strPtr++ = c; *strPtr++ = '\0'; /* * Free memory, and return. */ mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL); return k; } /* *---------------------------------------------------------------------- * * TclInitDoubleConversion -- * * Initializes constants that are needed for conversions to and from * 'double' * * Results: * None. * * Side effects: * The log base 2 of the floating point radix, the number of bits in a * double mantissa, and a table of the powers of five and ten are * computed and stored. * *---------------------------------------------------------------------- */ void TclInitDoubleConversion(void) { int i; int x; double d; if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) { Tcl_Panic("This code doesn't work on a decimal machine!"); } --log2FLT_RADIX; mantBits = DBL_MANT_DIG * log2FLT_RADIX; d = 1.0; x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0)); if (x < MAXPOW) { mmaxpow = x; } else { mmaxpow = MAXPOW; } for (i=0 ; i<=mmaxpow ; ++i) { pow10[i] = d; d *= 10.0; } for (i=0 ; i<9 ; ++i) { mp_init(pow5 + i); } mp_set(pow5, 5); for (i=0 ; i<8 ; ++i) { mp_sqr(pow5+i, pow5+i+1); } tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits); maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX) + log(10.)/2) / log(10.)); minDigits = (int) floor((DBL_MIN_EXP-DBL_MANT_DIG)*log((double)FLT_RADIX)/log(10.)); mantDIGIT = (mantBits + DIGIT_BIT - 1) / DIGIT_BIT; } /* *---------------------------------------------------------------------- * * TclFinalizeDoubleConversion -- * * Cleans up this file on exit. * * Results: * None * * Side effects: * Memory allocated by TclInitDoubleConversion is freed. * *---------------------------------------------------------------------- */ void TclFinalizeDoubleConversion() { int i; for (i=0 ; i<9 ; ++i) { mp_clear(pow5 + i); } } /* *---------------------------------------------------------------------- * * TclBignumToDouble -- * * Convert an arbitrary-precision integer to a native floating point * number. * * Results: * Returns the converted number. Sets errno to ERANGE if the number is * too large to convert. * *---------------------------------------------------------------------- */ double TclBignumToDouble(mp_int *a) /* Integer to convert. */ { mp_int b; int bits; int shift; int i; double r; /* * Determine how many bits we need, and extract that many from the input. * Round to nearest unit in the last place. */ bits = mp_count_bits(a); if (bits > DBL_MAX_EXP*log2FLT_RADIX) { errno = ERANGE; return HUGE_VAL; } shift = mantBits + 1 - bits; mp_init(&b); if (shift > 0) { mp_mul_2d(a, shift, &b); } else if (shift < 0) { mp_div_2d(a, -shift, &b, NULL); } else { mp_copy(a, &b); } mp_add_d(&b, 1, &b); mp_div_2d(&b, 1, &b, NULL); /* * Accumulate the result, one mp_digit at a time. */ r = 0.0; for (i=b.used-1 ; i>=0 ; --i) { r = ldexp(r, DIGIT_BIT) + b.dp[i]; } mp_clear(&b); /* * Scale the result to the correct number of bits. */ r = ldexp(r, bits - mantBits); /* * Return the result with the appropriate sign. */ if (a->sign == MP_ZPOS) { return r; } else { return -r; } } /* *---------------------------------------------------------------------- * * SafeLdExp -- * * Do an 'ldexp' operation, but handle denormals gracefully. * * Results: * Returns the appropriately scaled value. * * On some platforms, 'ldexp' fails when presented with a number too * small to represent as a normalized double. This routine does 'ldexp' * in two steps for those numbers, to return correctly denormalized * values. * *---------------------------------------------------------------------- */ static double SafeLdExp(double fract, int expt) { int minexpt = DBL_MIN_EXP * log2FLT_RADIX; volatile double a, b, retval; if (expt < minexpt) { a = ldexp(fract, expt - mantBits - minexpt); b = ldexp(1.0, mantBits + minexpt); retval = a * b; } else { retval = ldexp(fract, expt); } return retval; } /* *---------------------------------------------------------------------- * * TclFormatNaN -- * * Makes the string representation of a "Not a Number" * * Results: * None. * * Side effects: * Stores the string representation in the supplied buffer, which must be * at least TCL_DOUBLE_SPACE characters. * *---------------------------------------------------------------------- */ void TclFormatNaN(double value, /* The Not-a-Number to format. */ char *buffer) /* String representation. */ { #ifndef IEEE_FLOATING_POINT strcpy(buffer, "NaN"); return; #else union { double dv; Tcl_WideUInt iv; } bitwhack; bitwhack.dv = value; if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) { bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63); *buffer++ = '-'; } *buffer++ = 'N'; *buffer++ = 'a'; *buffer++ = 'N'; bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1; if (bitwhack.iv != 0) { sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv); } else { *buffer = '\0'; } #endif /* IEEE_FLOATING_POINT */ }