diff options
Diffstat (limited to 'generic/tclStrToD.c')
| -rwxr-xr-x | generic/tclStrToD.c | 1361 |
1 files changed, 1361 insertions, 0 deletions
diff --git a/generic/tclStrToD.c b/generic/tclStrToD.c new file mode 100755 index 0000000..2d622a7 --- /dev/null +++ b/generic/tclStrToD.c @@ -0,0 +1,1361 @@ +/* + *---------------------------------------------------------------------- + * + * 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.2 2005/05/10 18:34:49 kennykb Exp $ + * + *---------------------------------------------------------------------- + */ + +#include <tclInt.h> +#include <stdio.h> +#include <stdlib.h> +#include <float.h> +#include <limits.h> +#include <math.h> +#include <ctype.h> +#include <tommath.h> + +/* + * 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. + */ + +#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)) +#endif + +/* The powers of ten that can be represented exactly as IEEE754 doubles. */ + +#define MAXPOW 22 +static double pow10 [MAXPOW+1]; + +static int mmaxpow; /* Largest power of ten that can be + * represented exactly in a 'double'. */ + +/* Inexact higher powers of ten */ + +static CONST double pow_10_2_n [] = { + 1.0, + 100.0, + 10000.0, + 1.0e+8, + 1.0e+16, + 1.0e+32, + 1.0e+64, + 1.0e+128, + 1.0e+256 +}; + +/* Logarithm of the floating point radix. */ + +static int log2FLT_RADIX; + +/* Number of bits in a double's significand */ + +static int mantBits; + +/* Table of powers of 5**(2**n), up to 5**256 */ + +static mp_int pow5[9]; + +/* The smallest representable double */ + +static double tiny; + +/* The maximum number of digits to the left of the decimal point of a + * double. */ + +static int maxDigits; + +/* The maximum number of digits to the right of the decimal point in a + * double. */ + +static int minDigits; + +/* Number of mp_digit's needed to hold the significand of a double */ + +static int mantDIGIT; + +/* Static functions defined in this file */ + +static double RefineResult _ANSI_ARGS_((double approx, CONST char* start, + int nDigits, long exponent)); +static double ParseNaN _ANSI_ARGS_(( int signum, CONST char** end )); +static double SafeLdExp _ANSI_ARGS_(( 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 */ + + /* + * v must be 'volatile double' on gc-ix86 to force correct rounding + * to IEEE double and not Intel double-extended. + */ + + volatile double v; /* Scanned value */ + 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 + * 1 ulp, so we need to change rounding mode to 53-bits. + */ + +#if defined(__GNUC__) && defined(__i386) + fpu_control_t roundTo53Bits = 0x027f; + fpu_control_t oldRoundingMode; + _FPU_GETCW( oldRoundingMode ); + _FPU_SETCW( roundTo53Bits ); +#endif + + /* Discard leading whitespace */ + + while ( isspace( *p ) ) { + ++p; + } + + /* Determine the sign of the significand */ + + switch( *p ) { + case '-': + signum = 1; + /* FALLTHROUGH */ + case '+': + ++p; + } + + /* Discard leading zeroes */ + + 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(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( 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 = 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 */ + + 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 + + /* 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 */ + +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW( oldRoundingMode ); +#endif + 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 */ + +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW( oldRoundingMode ); +#endif + + return v; + + /* Come here on an invalid input */ + + error: + if ( endPtr != NULL ) { + *endPtr = s; + } + + /* Restore FPU mode word */ + +#if defined(__GNUC__) && defined(__i386) + _FPU_SETCW( oldRoundingMode ); +#endif + 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( 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 ) { + 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 < minincr ) { + quot = minincr; + } + + mp_clear( &twoMd ); + mp_clear( &twoMv ); + + + return approxResult + quot; +} + +/* + *---------------------------------------------------------------------- + * + * ParseNaN -- + * + * Parses a "not a number" from an input string, and returns the + * double precision NaN corresponding to it. + * + * Side effects: + * Advances endPtr to follow any (hex) in the input string. + * + * If the NaN is followed by a left paren, a string of spaes + * and hexadecimal digits, and a right paren, endPtr is advanced + * to follow it. + * + * The string of hexadecimal digits is OR'ed into the resulting + * NaN, and the signum is set as well. Note that a signalling NaN + * is never returned. + * + *---------------------------------------------------------------------- + */ + +double +ParseNaN( int signum, /* Flag == 1 if minus sign has been + * seen in front of NaN */ + CONST char** endPtr ) + /* Pointer-to-pointer to char following "NaN" + * in the input string */ +{ + CONST char* p = *endPtr; + char c; + union { + Tcl_WideUInt iv; + double dv; + } theNaN; + + /* Scan off a hex number in parentheses. Embedded blanks are ok. */ + + theNaN.iv = 0; + if ( *p == '(' ) { + ++p; + for ( ; ; ) { + c = *p++; + if ( isspace(c) ) { + continue; + } else if ( c == ')' ) { + *endPtr = p; + break; + } else if ( isdigit(c) ) { + c -= '0'; + } else if ( c >= 'A' && c <= 'F' ) { + c = c - 'A' + 10; + } else if ( c >= 'a' && c <= 'f' ) { + c = c - 'a' + 10; + } else { + theNaN.iv = ( ((Tcl_WideUInt) 0x7ff8) << 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) 0x7ff8) << 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 low_ok; + int high_ok; + + 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; + + low_ok = high_ok = ( 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 || ( high_ok && 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 || ( high_ok && 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 ( low_ok ) { + tc1 = ( tc1 <= 0 ); + } else { + tc1 = ( tc1 < 0 ); + } + mp_add( &r, &mplus, &temp ); + tc2 = mp_cmp_mag( &temp, &s ); + if ( high_ok ) { + 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 */ + + 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) + + 0.5 * log(10.)) + / 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 +} |