Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

source: downloads/tcl8.5.2/generic/tclStrToD.c @ 25

Last change on this file since 25 was 25, checked in by landauf, 18 years ago
added tcl to libs
File size: 68.8 KB

Rev	Line
[25]	1	/*
	2	*----------------------------------------------------------------------
	3	*
	4	* tclStrToD.c --
	5	*
	6	* This file contains a collection of procedures for managing conversions
	7	* to/from floating-point in Tcl. They include TclParseNumber, which
	8	* parses numbers from strings; TclDoubleDigits, which formats numbers
	9	* into strings of digits, and procedures for interconversion among
	10	* 'double' and 'mp_int' types.
	11	*
	12	* Copyright (c) 2005 by Kevin B. Kenny. All rights reserved.
	13	*
	14	* See the file "license.terms" for information on usage and redistribution of
	15	* this file, and for a DISCLAIMER OF ALL WARRANTIES.
	16	*
	17	* RCS: @(#) $Id: tclStrToD.c,v 1.33 2008/03/13 17:14:19 dgp Exp $
	18	*
	19	*----------------------------------------------------------------------
	20	*/
	21
	22	#include <tclInt.h>
	23	#include <stdio.h>
	24	#include <stdlib.h>
	25	#include <float.h>
	26	#include <limits.h>
	27	#include <math.h>
	28	#include <ctype.h>
	29	#include <tommath.h>
	30
	31	/*
	32	* Define KILL_OCTAL to suppress interpretation of numbers with leading zero
	33	* as octal. (Ceterum censeo: numeros octonarios delendos esse.)
	34	*/
	35
	36	#undef KILL_OCTAL
	37
	38	/*
	39	* This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754
	40	* floating point; of these, only IEEE-754 can represent NaN. IEEE-754 can be
	41	* uniquely determined by radix and by the widths of significand and exponent.
	42	*/
	43
	44	#if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024)
	45	# define IEEE_FLOATING_POINT
	46	#endif
	47
	48	/*
	49	* gcc on x86 needs access to rounding controls, because of a questionable
	50	* feature where it retains intermediate results as IEEE 'long double' values
	51	* somewhat unpredictably. It is tempting to include fpu_control.h, but that
	52	* file exists only on Linux; it is missing on Cygwin and MinGW. Most gcc-isms
	53	* and ix86-isms are factored out here.
	54	*/
	55
	56	#if defined(__GNUC__) && defined(__i386)
	57	typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
	58	#define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
	59	#define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
	60	# define FPU_IEEE_ROUNDING 0x027f
	61	# define ADJUST_FPU_CONTROL_WORD
	62	#endif
	63
	64	/*
	65	* HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN.
	66	* Everyone else uses 7ff8000000000000. (Why, HP, why?)
	67	*/
	68
	69	#ifdef __hppa
	70	# define NAN_START 0x7ff4
	71	# define NAN_MASK (((Tcl_WideUInt) 1) << 50)
	72	#else
	73	# define NAN_START 0x7ff8
	74	# define NAN_MASK (((Tcl_WideUInt) 1) << 51)
	75	#endif
	76
	77	/*
	78	* Constants used by this file (most of which are only ever calculated at
	79	* runtime).
	80	*/
	81
	82	static int maxpow10_wide; /* The powers of ten that can be represented
	83	* exactly as wide integers. */
	84	static Tcl_WideUInt *pow10_wide;
	85	#define MAXPOW 22
	86	static double pow10vals[MAXPOW+1]; /* The powers of ten that can be represented
	87	* exactly as IEEE754 doubles. */
	88	static int mmaxpow; /* Largest power of ten that can be
	89	* represented exactly in a 'double'. */
	90	static int log10_DIGIT_MAX; /* The number of decimal digits that fit in an
	91	* mp_digit. */
	92	static int log2FLT_RADIX; /* Logarithm of the floating point radix. */
	93	static int mantBits; /* Number of bits in a double's significand */
	94	static mp_int pow5[9]; /* Table of powers of 5(2n), up to
	95	* 5*256 /
	96	static double tiny; /* The smallest representable double */
	97	static int maxDigits; /* The maximum number of digits to the left of
	98	* the decimal point of a double. */
	99	static int minDigits; /* The maximum number of digits to the right
	100	* of the decimal point in a double. */
	101	static int mantDIGIT; /* Number of mp_digit's needed to hold the
	102	* significand of a double. */
	103	static const double pow_10_2_n[] = { /* Inexact higher powers of ten. */
	104	1.0,
	105	100.0,
	106	10000.0,
	107	1.0e+8,
	108	1.0e+16,
	109	1.0e+32,
	110	1.0e+64,
	111	1.0e+128,
	112	1.0e+256
	113	};
	114	static int n770_fp; /* Flag is 1 on Nokia N770 floating point.
	115	* Nokia's floating point has the words
	116	* reversed: if big-endian is 7654 3210,
	117	* and little-endian is 0123 4567,
	118	* then Nokia's FP is 4567 0123;
	119	* little-endian within the 32-bit words
	120	* but big-endian between them. */
	121
	122	/*
	123	* Static functions defined in this file.
	124	*/
	125
	126	static double AbsoluteValue(double v, int *signum);
	127	static int AccumulateDecimalDigit(unsigned, int,
	128	Tcl_WideUInt , mp_int , int);
	129	static double BignumToBiasedFrExp(mp_int big, int machexp);
	130	static int GetIntegerTimesPower(double v, mp_int r, int e);
	131	static double MakeHighPrecisionDouble(int signum,
	132	mp_int *significand, int nSigDigs, int exponent);
	133	static double MakeLowPrecisionDouble(int signum,
	134	Tcl_WideUInt significand, int nSigDigs,
	135	int exponent);
	136	static double MakeNaN(int signum, Tcl_WideUInt tag);
	137	static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
	138	static double Pow10TimesFrExp(int exponent, double fraction,
	139	int *machexp);
	140	static double RefineApproximation(double approx,
	141	mp_int *exactSignificand, int exponent);
	142	static double SafeLdExp(double fraction, int exponent);
	143
	144	/*
	145	*----------------------------------------------------------------------
	146	*
	147	* TclParseNumber --
	148	*
	149	* Scans bytes, interpreted as characters in Tcl's internal encoding, and
	150	* parses the longest prefix that is the string representation of a
	151	* number in a format recognized by Tcl.
	152	*
	153	* The arguments bytes, numBytes, and objPtr are the inputs which
	154	* determine the string to be parsed. If bytes is non-NULL, it points to
	155	* the first byte to be scanned. If bytes is NULL, then objPtr must be
	156	* non-NULL, and the string representation of objPtr will be scanned
	157	* (generated first, if necessary). The numBytes argument determines the
	158	* number of bytes to be scanned. If numBytes is negative, the first NUL
	159	* byte encountered will terminate the scan. If numBytes is non-negative,
	160	* then no more than numBytes bytes will be scanned.
	161	*
	162	* The argument flags is an input that controls the numeric formats
	163	* recognized by the parser. The flag bits are:
	164	*
	165	* - TCL_PARSE_INTEGER_ONLY: accept only integer values; reject
	166	* strings that denote floating point values (or accept only the
	167	* leading portion of them that are integer values).
	168	* - TCL_PARSE_SCAN_PREFIXES: ignore the prefixes 0b and 0o that are
	169	* not part of the [scan] command's vocabulary. Use only in
	170	* combination with TCL_PARSE_INTEGER_ONLY.
	171	* - TCL_PARSE_OCTAL_ONLY: parse only in the octal format, whether
	172	* or not a prefix is present that would lead to octal parsing.
	173	* Use only in combination with TCL_PARSE_INTEGER_ONLY.
	174	* - TCL_PARSE_HEXADECIMAL_ONLY: parse only in the hexadecimal format,
	175	* whether or not a prefix is present that would lead to
	176	* hexadecimal parsing. Use only in combination with
	177	* TCL_PARSE_INTEGER_ONLY.
	178	* - TCL_PARSE_DECIMAL_ONLY: parse only in the decimal format, no
	179	* matter whether a 0 prefix would normally force a different
	180	* base.
	181	* - TCL_PARSE_NO_WHITESPACE: reject any leading/trailing whitespace
	182	*
	183	* The arguments interp and expected are inputs that control error
	184	* message generation. If interp is NULL, no error message will be
	185	* generated. If interp is non-NULL, then expected must also be non-NULL.
	186	* When TCL_ERROR is returned, an error message will be left in the
	187	* result of interp, and the expected argument will appear in the error
	188	* message as the thing TclParseNumber expected, but failed to find in
	189	* the string.
	190	*
	191	* The arguments objPtr and endPtrPtr as well as the return code are the
	192	* outputs.
	193	*
	194	* When the parser cannot find any prefix of the string that matches a
	195	* format it is looking for, TCL_ERROR is returned and an error message
	196	* may be generated and returned as described above. The contents of
	197	* objPtr will not be changed. If endPtrPtr is non-NULL, a pointer to the
	198	* character in the string that terminated the scan will be written to
	199	* *endPtrPtr.
	200	*
	201	* When the parser determines that the entire string matches a format it
	202	* is looking for, TCL_OK is returned, and if objPtr is non-NULL, then
	203	* the internal rep and Tcl_ObjType of objPtr are set to the "canonical"
	204	* numeric value that matches the scanned string. If endPtrPtr is not
	205	* NULL, a pointer to the end of the string will be written to *endPtrPtr
	206	* (that is, either bytes+numBytes or a pointer to a terminating NUL
	207	* byte).
	208	*
	209	* When the parser determines that a partial string matches a format it
	210	* is looking for, the value of endPtrPtr determines what happens:
	211	*
	212	* - If endPtrPtr is NULL, then TCL_ERROR is returned, with error message
	213	* generation as above.
	214	*
	215	* - If endPtrPtr is non-NULL, then TCL_OK is returned and objPtr
	216	* internals are set as above. Also, a pointer to the first
	217	* character following the parsed numeric string is written to
	218	* *endPtrPtr.
	219	*
	220	* In some cases where the string being scanned is the string rep of
	221	* objPtr, this routine can leave objPtr in an inconsistent state where
	222	* its string rep and its internal rep do not agree. In these cases the
	223	* internal rep will be in agreement with only some substring of the
	224	* string rep. This might happen if the caller passes in a non-NULL bytes
	225	* value that points somewhere into the string rep. It might happen if
	226	* the caller passes in a numBytes value that limits the scan to only a
	227	* prefix of the string rep. Or it might happen if a non-NULL value of
	228	* endPtrPtr permits a TCL_OK return from only a partial string match. It
	229	* is the responsibility of the caller to detect and correct such
	230	* inconsistencies when they can and do arise.
	231	*
	232	* Results:
	233	* Returns a standard Tcl result.
	234	*
	235	* Side effects:
	236	* The string representaton of objPtr may be generated.
	237	*
	238	* The internal representation and Tcl_ObjType of objPtr may be changed.
	239	* This may involve allocation and/or freeing of memory.
	240	*
	241	*----------------------------------------------------------------------
	242	*/
	243
	244	int
	245	TclParseNumber(
	246	Tcl_Interp interp, / Used for error reporting. May be NULL. */
	247	Tcl_Obj objPtr, / Object to receive the internal rep. */
	248	const char expected, / Description of the type of number the
	249	* caller expects to be able to parse
	250	* ("integer", "boolean value", etc.). */
	251	const char bytes, / Pointer to the start of the string to
	252	* scan. */
	253	int numBytes, /* Maximum number of bytes to scan, see
	254	* above. */
	255	const char *endPtrPtr, / Place to store pointer to the character
	256	* that terminated the scan. */
	257	int flags) /* Flags governing the parse. */
	258	{
	259	enum State {
	260	INITIAL, SIGNUM, ZERO, ZERO_X,
	261	ZERO_O, ZERO_B, BINARY,
	262	HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL,
	263	LEADING_RADIX_POINT, FRACTION,
	264	EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
	265	sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
	266	#ifdef IEEE_FLOATING_POINT
	267	, sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
	268	#endif
	269	} state = INITIAL;
	270	enum State acceptState = INITIAL;
	271
	272	int signum = 0; /* Sign of the number being parsed */
	273	Tcl_WideUInt significandWide = 0;
	274	/* Significand of the number being parsed (if
	275	* no overflow) */
	276	mp_int significandBig; /* Significand of the number being parsed (if
	277	* it overflows significandWide) */
	278	int significandOverflow = 0;/* Flag==1 iff significandBig is used */
	279	Tcl_WideUInt octalSignificandWide = 0;
	280	/* Significand of an octal number; needed
	281	* because we don't know whether a number with
	282	* a leading zero is octal or decimal until
	283	* we've scanned forward to a '.' or 'e' */
	284	mp_int octalSignificandBig; /* Significand of octal number once
	285	* octalSignificandWide overflows */
	286	int octalSignificandOverflow = 0;
	287	/* Flag==1 if octalSignificandBig is used */
	288	int numSigDigs = 0; /* Number of significant digits in the decimal
	289	* significand */
	290	int numTrailZeros = 0; /* Number of trailing zeroes at the current
	291	* point in the parse. */
	292	int numDigitsAfterDp = 0; /* Number of digits scanned after the decimal
	293	* point */
	294	int exponentSignum = 0; /* Signum of the exponent of a floating point
	295	* number */
	296	long exponent = 0; /* Exponent of a floating point number */
	297	const char p; / Pointer to next character to scan */
	298	size_t len; /* Number of characters remaining after p */
	299	const char acceptPoint; / Pointer to position after last character in
	300	* an acceptable number */
	301	size_t acceptLen; /* Number of characters following that
	302	* point. */
	303	int status = TCL_OK; /* Status to return to caller */
	304	char d = 0; /* Last hexadecimal digit scanned; initialized
	305	* to avoid a compiler warning. */
	306	int shift = 0; /* Amount to shift when accumulating binary */
	307	int explicitOctal = 0;
	308
	309	#define ALL_BITS (~(Tcl_WideUInt)0)
	310	#define MOST_BITS (ALL_BITS >> 1)
	311
	312	/*
	313	* Initialize bytes to start of the object's string rep if the caller
	314	* didn't pass anything else.
	315	*/
	316
	317	if (bytes == NULL) {
	318	bytes = TclGetString(objPtr);
	319	}
	320
	321	p = bytes;
	322	len = numBytes;
	323	acceptPoint = p;
	324	acceptLen = len;
	325	while (1) {
	326	char c = len ? *p : '\0';
	327	switch (state) {
	328
	329	case INITIAL:
	330	/*
	331	* Initial state. Acceptable characters are +, -, digits, period,
	332	* I, N, and whitespace.
	333	*/
	334
	335	if (isspace(UCHAR(c))) {
	336	if (flags & TCL_PARSE_NO_WHITESPACE) {
	337	goto endgame;
	338	}
	339	break;
	340	} else if (c == '+') {
	341	state = SIGNUM;
	342	break;
	343	} else if (c == '-') {
	344	signum = 1;
	345	state = SIGNUM;
	346	break;
	347	}
	348	/* FALLTHROUGH */
	349
	350	case SIGNUM:
	351	/*
	352	* Scanned a leading + or -. Acceptable characters are digits,
	353	* period, I, and N.
	354	*/
	355
	356	if (c == '0') {
	357	if (flags & TCL_PARSE_DECIMAL_ONLY) {
	358	state = DECIMAL;
	359	} else {
	360	state = ZERO;
	361	}
	362	break;
	363	} else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
	364	goto zerox;
	365	} else if (flags & TCL_PARSE_OCTAL_ONLY) {
	366	goto zeroo;
	367	} else if (isdigit(UCHAR(c))) {
	368	significandWide = c - '0';
	369	numSigDigs = 1;
	370	state = DECIMAL;
	371	break;
	372	} else if (flags & TCL_PARSE_INTEGER_ONLY) {
	373	goto endgame;
	374	} else if (c == '.') {
	375	state = LEADING_RADIX_POINT;
	376	break;
	377	} else if (c == 'I' \|\| c == 'i') {
	378	state = sI;
	379	break;
	380	#ifdef IEEE_FLOATING_POINT
	381	} else if (c == 'N' \|\| c == 'n') {
	382	state = sN;
	383	break;
	384	#endif
	385	}
	386	goto endgame;
	387
	388	case ZERO:
	389	/*
	390	* Scanned a leading zero (perhaps with a + or -). Acceptable
	391	* inputs are digits, period, X, and E. If 8 or 9 is encountered,
	392	* the number can't be octal. This state and the OCTAL state
	393	* differ only in whether they recognize 'X'.
	394	*/
	395
	396	acceptState = state;
	397	acceptPoint = p;
	398	acceptLen = len;
	399	if (c == 'x' \|\| c == 'X') {
	400	state = ZERO_X;
	401	break;
	402	}
	403	if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
	404	goto zerox;
	405	}
	406	if (flags & TCL_PARSE_SCAN_PREFIXES) {
	407	goto zeroo;
	408	}
	409	if (c == 'b' \|\| c == 'B') {
	410	state = ZERO_B;
	411	break;
	412	}
	413	if (c == 'o' \|\| c == 'O') {
	414	explicitOctal = 1;
	415	state = ZERO_O;
	416	break;
	417	}
	418	#ifdef KILL_OCTAL
	419	goto decimal;
	420	#endif
	421	/* FALLTHROUGH */
	422
	423	case OCTAL:
	424	/*
	425	* Scanned an optional + or -, followed by a string of octal
	426	* digits. Acceptable inputs are more digits, period, or E. If 8
	427	* or 9 is encountered, commit to floating point.
	428	*/
	429
	430	acceptState = state;
	431	acceptPoint = p;
	432	acceptLen = len;
	433	/* FALLTHROUGH */
	434	case ZERO_O:
	435	zeroo:
	436	if (c == '0') {
	437	++numTrailZeros;
	438	state = OCTAL;
	439	break;
	440	} else if (c >= '1' && c <= '7') {
	441	if (objPtr != NULL) {
	442	shift = 3 * (numTrailZeros + 1);
	443	significandOverflow = AccumulateDecimalDigit(
	444	(unsigned)(c-'0'), numTrailZeros,
	445	&significandWide, &significandBig,
	446	significandOverflow);
	447
	448	if (!octalSignificandOverflow) {
	449	/*
	450	* Shifting by more bits than are in the value being
	451	* shifted is at least de facto nonportable. Check for
	452	* too large shifts first.
	453	*/
	454
	455	if ((octalSignificandWide != 0)
	456	&& (((size_t)shift >=
	457	CHAR_BIT*sizeof(Tcl_WideUInt))
	458	\|\| (octalSignificandWide >
	459	(~(Tcl_WideUInt)0 >> shift)))) {
	460	octalSignificandOverflow = 1;
	461	TclBNInitBignumFromWideUInt(&octalSignificandBig,
	462	octalSignificandWide);
	463	}
	464	}
	465	if (!octalSignificandOverflow) {
	466	octalSignificandWide =
	467	(octalSignificandWide << shift) + (c - '0');
	468	} else {
	469	mp_mul_2d(&octalSignificandBig, shift,
	470	&octalSignificandBig);
	471	mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
	472	&octalSignificandBig);
	473	}
	474	}
	475	if (numSigDigs != 0) {
	476	numSigDigs += numTrailZeros+1;
	477	} else {
	478	numSigDigs = 1;
	479	}
	480	numTrailZeros = 0;
	481	state = OCTAL;
	482	break;
	483	}
	484	/* FALLTHROUGH */
	485
	486	case BAD_OCTAL:
	487	if (explicitOctal) {
	488	/*
	489	* No forgiveness for bad digits in explicitly octal numbers.
	490	*/
	491
	492	goto endgame;
	493	}
	494	if (flags & TCL_PARSE_INTEGER_ONLY) {
	495	/*
	496	* No seeking floating point when parsing only integer.
	497	*/
	498
	499	goto endgame;
	500	}
	501	#ifndef KILL_OCTAL
	502
	503	/*
	504	* Scanned a number with a leading zero that contains an 8, 9,
	505	* radix point or E. This is an invalid octal number, but might
	506	* still be floating point.
	507	*/
	508
	509	if (c == '0') {
	510	++numTrailZeros;
	511	state = BAD_OCTAL;
	512	break;
	513	} else if (isdigit(UCHAR(c))) {
	514	if (objPtr != NULL) {
	515	significandOverflow = AccumulateDecimalDigit(
	516	(unsigned)(c-'0'), numTrailZeros,
	517	&significandWide, &significandBig,
	518	significandOverflow);
	519	}
	520	if (numSigDigs != 0) {
	521	numSigDigs += (numTrailZeros + 1);
	522	} else {
	523	numSigDigs = 1;
	524	}
	525	numTrailZeros = 0;
	526	state = BAD_OCTAL;
	527	break;
	528	} else if (c == '.') {
	529	state = FRACTION;
	530	break;
	531	} else if (c == 'E' \|\| c == 'e') {
	532	state = EXPONENT_START;
	533	break;
	534	}
	535	#endif
	536	goto endgame;
	537
	538	/*
	539	* Scanned 0x. If state is HEXADECIMAL, scanned at least one
	540	* character following the 0x. The only acceptable inputs are
	541	* hexadecimal digits.
	542	*/
	543
	544	case HEXADECIMAL:
	545	acceptState = state;
	546	acceptPoint = p;
	547	acceptLen = len;
	548	/* FALLTHROUGH */
	549
	550	case ZERO_X:
	551	zerox:
	552	if (c == '0') {
	553	++numTrailZeros;
	554	state = HEXADECIMAL;
	555	break;
	556	} else if (isdigit(UCHAR(c))) {
	557	d = (c-'0');
	558	} else if (c >= 'A' && c <= 'F') {
	559	d = (c-'A'+10);
	560	} else if (c >= 'a' && c <= 'f') {
	561	d = (c-'a'+10);
	562	} else {
	563	goto endgame;
	564	}
	565	if (objPtr != NULL) {
	566	shift = 4 * (numTrailZeros + 1);
	567	if (!significandOverflow) {
	568	/*
	569	* Shifting by more bits than are in the value being
	570	* shifted is at least de facto nonportable. Check for too
	571	* large shifts first.
	572	*/
	573
	574	if (significandWide != 0 &&
	575	((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) \|\|
	576	significandWide > (~(Tcl_WideUInt)0 >> shift))) {
	577	significandOverflow = 1;
	578	TclBNInitBignumFromWideUInt(&significandBig,
	579	significandWide);
	580	}
	581	}
	582	if (!significandOverflow) {
	583	significandWide = (significandWide << shift) + d;
	584	} else {
	585	mp_mul_2d(&significandBig, shift, &significandBig);
	586	mp_add_d(&significandBig, (mp_digit) d, &significandBig);
	587	}
	588	}
	589	numTrailZeros = 0;
	590	state = HEXADECIMAL;
	591	break;
	592
	593	case BINARY:
	594	acceptState = state;
	595	acceptPoint = p;
	596	acceptLen = len;
	597	case ZERO_B:
	598	if (c == '0') {
	599	++numTrailZeros;
	600	state = BINARY;
	601	break;
	602	} else if (c != '1') {
	603	goto endgame;
	604	}
	605	if (objPtr != NULL) {
	606	shift = numTrailZeros + 1;
	607	if (!significandOverflow) {
	608	/*
	609	* Shifting by more bits than are in the value being
	610	* shifted is at least de facto nonportable. Check for too
	611	* large shifts first.
	612	*/
	613
	614	if (significandWide != 0 &&
	615	((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) \|\|
	616	significandWide > (~(Tcl_WideUInt)0 >> shift))) {
	617	significandOverflow = 1;
	618	TclBNInitBignumFromWideUInt(&significandBig,
	619	significandWide);
	620	}
	621	}
	622	if (!significandOverflow) {
	623	significandWide = (significandWide << shift) + 1;
	624	} else {
	625	mp_mul_2d(&significandBig, shift, &significandBig);
	626	mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
	627	}
	628	}
	629	numTrailZeros = 0;
	630	state = BINARY;
	631	break;
	632
	633	case DECIMAL:
	634	/*
	635	* Scanned an optional + or - followed by a string of decimal
	636	* digits.
	637	*/
	638
	639	#ifdef KILL_OCTAL
	640	decimal:
	641	#endif
	642	acceptState = state;
	643	acceptPoint = p;
	644	acceptLen = len;
	645	if (c == '0') {
	646	++numTrailZeros;
	647	state = DECIMAL;
	648	break;
	649	} else if (isdigit(UCHAR(c))) {
	650	if (objPtr != NULL) {
	651	significandOverflow = AccumulateDecimalDigit(
	652	(unsigned)(c - '0'), numTrailZeros,
	653	&significandWide, &significandBig,
	654	significandOverflow);
	655	}
	656	numSigDigs += numTrailZeros+1;
	657	numTrailZeros = 0;
	658	state = DECIMAL;
	659	break;
	660	} else if (flags & TCL_PARSE_INTEGER_ONLY) {
	661	goto endgame;
	662	} else if (c == '.') {
	663	state = FRACTION;
	664	break;
	665	} else if (c == 'E' \|\| c == 'e') {
	666	state = EXPONENT_START;
	667	break;
	668	}
	669	goto endgame;
	670
	671	/*
	672	* Found a decimal point. If no digits have yet been scanned, E is
	673	* not allowed; otherwise, it introduces the exponent. If at least
	674	* one digit has been found, we have a possible complete number.
	675	*/
	676
	677	case FRACTION:
	678	acceptState = state;
	679	acceptPoint = p;
	680	acceptLen = len;
	681	if (c == 'E' \|\| c=='e') {
	682	state = EXPONENT_START;
	683	break;
	684	}
	685	/* FALLTHROUGH */
	686
	687	case LEADING_RADIX_POINT:
	688	if (c == '0') {
	689	++numDigitsAfterDp;
	690	++numTrailZeros;
	691	state = FRACTION;
	692	break;
	693	} else if (isdigit(UCHAR(c))) {
	694	++numDigitsAfterDp;
	695	if (objPtr != NULL) {
	696	significandOverflow = AccumulateDecimalDigit(
	697	(unsigned)(c-'0'), numTrailZeros,
	698	&significandWide, &significandBig,
	699	significandOverflow);
	700	}
	701	if (numSigDigs != 0) {
	702	numSigDigs += numTrailZeros+1;
	703	} else {
	704	numSigDigs = 1;
	705	}
	706	numTrailZeros = 0;
	707	state = FRACTION;
	708	break;
	709	}
	710	goto endgame;
	711
	712	case EXPONENT_START:
	713	/*
	714	* Scanned the E at the start of an exponent. Make sure a legal
	715	* character follows before using the C library strtol routine,
	716	* which allows whitespace.
	717	*/
	718
	719	if (c == '+') {
	720	state = EXPONENT_SIGNUM;
	721	break;
	722	} else if (c == '-') {
	723	exponentSignum = 1;
	724	state = EXPONENT_SIGNUM;
	725	break;
	726	}
	727	/* FALLTHROUGH */
	728
	729	case EXPONENT_SIGNUM:
	730	/*
	731	* Found the E at the start of the exponent, followed by a sign
	732	* character.
	733	*/
	734
	735	if (isdigit(UCHAR(c))) {
	736	exponent = c - '0';
	737	state = EXPONENT;
	738	break;
	739	}
	740	goto endgame;
	741
	742	case EXPONENT:
	743	/*
	744	* Found an exponent with at least one digit. Accumulate it,
	745	* making sure to hard-pin it to LONG_MAX on overflow.
	746	*/
	747
	748	acceptState = state;
	749	acceptPoint = p;
	750	acceptLen = len;
	751	if (isdigit(UCHAR(c))) {
	752	if (exponent < (LONG_MAX - 9) / 10) {
	753	exponent = 10 * exponent + (c - '0');
	754	} else {
	755	exponent = LONG_MAX;
	756	}
	757	state = EXPONENT;
	758	break;
	759	}
	760	goto endgame;
	761
	762	/*
	763	* Parse out INFINITY by simply spelling it out. INF is accepted
	764	* as an abbreviation; other prefices are not.
	765	*/
	766
	767	case sI:
	768	if (c == 'n' \|\| c == 'N') {
	769	state = sIN;
	770	break;
	771	}
	772	goto endgame;
	773	case sIN:
	774	if (c == 'f' \|\| c == 'F') {
	775	state = sINF;
	776	break;
	777	}
	778	goto endgame;
	779	case sINF:
	780	acceptState = state;
	781	acceptPoint = p;
	782	acceptLen = len;
	783	if (c == 'i' \|\| c == 'I') {
	784	state = sINFI;
	785	break;
	786	}
	787	goto endgame;
	788	case sINFI:
	789	if (c == 'n' \|\| c == 'N') {
	790	state = sINFIN;
	791	break;
	792	}
	793	goto endgame;
	794	case sINFIN:
	795	if (c == 'i' \|\| c == 'I') {
	796	state = sINFINI;
	797	break;
	798	}
	799	goto endgame;
	800	case sINFINI:
	801	if (c == 't' \|\| c == 'T') {
	802	state = sINFINIT;
	803	break;
	804	}
	805	goto endgame;
	806	case sINFINIT:
	807	if (c == 'y' \|\| c == 'Y') {
	808	state = sINFINITY;
	809	break;
	810	}
	811	goto endgame;
	812
	813	/*
	814	* Parse NaN's.
	815	*/
	816	#ifdef IEEE_FLOATING_POINT
	817	case sN:
	818	if (c == 'a' \|\| c == 'A') {
	819	state = sNA;
	820	break;
	821	}
	822	goto endgame;
	823	case sNA:
	824	if (c == 'n' \|\| c == 'N') {
	825	state = sNAN;
	826	break;
	827	}
	828	goto endgame;
	829	case sNAN:
	830	acceptState = state;
	831	acceptPoint = p;
	832	acceptLen = len;
	833	if (c == '(') {
	834	state = sNANPAREN;
	835	break;
	836	}
	837	goto endgame;
	838
	839	/*
	840	* Parse NaN(hexdigits)
	841	*/
	842	case sNANHEX:
	843	if (c == ')') {
	844	state = sNANFINISH;
	845	break;
	846	}
	847	/* FALLTHROUGH */
	848	case sNANPAREN:
	849	if (isspace(UCHAR(c))) {
	850	break;
	851	}
	852	if (numSigDigs < 13) {
	853	if (c >= '0' && c <= '9') {
	854	d = c - '0';
	855	} else if (c >= 'a' && c <= 'f') {
	856	d = 10 + c - 'a';
	857	} else if (c >= 'A' && c <= 'F') {
	858	d = 10 + c - 'A';
	859	}
	860	significandWide = (significandWide << 4) + d;
	861	state = sNANHEX;
	862	break;
	863	}
	864	goto endgame;
	865	case sNANFINISH:
	866	#endif
	867
	868	case sINFINITY:
	869	acceptState = state;
	870	acceptPoint = p;
	871	acceptLen = len;
	872	goto endgame;
	873	}
	874	++p;
	875	--len;
	876	}
	877
	878	endgame:
	879	if (acceptState == INITIAL) {
	880	/*
	881	* No numeric string at all found.
	882	*/
	883
	884	status = TCL_ERROR;
	885	if (endPtrPtr != NULL) {
	886	*endPtrPtr = p;
	887	}
	888	} else {
	889	/*
	890	* Back up to the last accepting state in the lexer.
	891	*/
	892
	893	p = acceptPoint;
	894	len = acceptLen;
	895	if (!(flags & TCL_PARSE_NO_WHITESPACE)) {
	896	/*
	897	* Accept trailing whitespace.
	898	*/
	899
	900	while (len != 0 && isspace(UCHAR(*p))) {
	901	++p;
	902	--len;
	903	}
	904	}
	905	if (endPtrPtr == NULL) {
	906	if ((len != 0) && ((numBytes > 0) \|\| (*p != '\0'))) {
	907	status = TCL_ERROR;
	908	}
	909	} else {
	910	*endPtrPtr = p;
	911	}
	912	}
	913
	914	/*
	915	* Generate and store the appropriate internal rep.
	916	*/
	917
	918	if (status == TCL_OK && objPtr != NULL) {
	919	TclFreeIntRep(objPtr);
	920	switch (acceptState) {
	921	case SIGNUM:
	922	case BAD_OCTAL:
	923	case ZERO_X:
	924	case ZERO_O:
	925	case ZERO_B:
	926	case LEADING_RADIX_POINT:
	927	case EXPONENT_START:
	928	case EXPONENT_SIGNUM:
	929	case sI:
	930	case sIN:
	931	case sINFI:
	932	case sINFIN:
	933	case sINFINI:
	934	case sINFINIT:
	935	case sN:
	936	case sNA:
	937	case sNANPAREN:
	938	case sNANHEX:
	939	Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'",
	940	acceptState, bytes);
	941
	942	case BINARY:
	943	shift = numTrailZeros;
	944	if (!significandOverflow && significandWide != 0 &&
	945	((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) \|\|
	946	significandWide > (MOST_BITS + signum) >> shift)) {
	947	significandOverflow = 1;
	948	TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	949	}
	950	if (shift) {
	951	if (!significandOverflow) {
	952	significandWide <<= shift;
	953	} else {
	954	mp_mul_2d(&significandBig, shift, &significandBig);
	955	}
	956	}
	957	goto returnInteger;
	958
	959	case HEXADECIMAL:
	960	/*
	961	* Returning a hex integer. Final scaling step.
	962	*/
	963
	964	shift = 4 * numTrailZeros;
	965	if (!significandOverflow && significandWide !=0 &&
	966	((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) \|\|
	967	significandWide > (MOST_BITS + signum) >> shift)) {
	968	significandOverflow = 1;
	969	TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	970	}
	971	if (shift) {
	972	if (!significandOverflow) {
	973	significandWide <<= shift;
	974	} else {
	975	mp_mul_2d(&significandBig, shift, &significandBig);
	976	}
	977	}
	978	goto returnInteger;
	979
	980	case OCTAL:
	981	/*
	982	* Returning an octal integer. Final scaling step
	983	*/
	984
	985	shift = 3 * numTrailZeros;
	986	if (!octalSignificandOverflow && octalSignificandWide != 0 &&
	987	((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) \|\|
	988	octalSignificandWide > (MOST_BITS + signum) >> shift)) {
	989	octalSignificandOverflow = 1;
	990	TclBNInitBignumFromWideUInt(&octalSignificandBig,
	991	octalSignificandWide);
	992	}
	993	if (shift) {
	994	if (!octalSignificandOverflow) {
	995	octalSignificandWide <<= shift;
	996	} else {
	997	mp_mul_2d(&octalSignificandBig, shift,
	998	&octalSignificandBig);
	999	}
	1000	}
	1001	if (!octalSignificandOverflow) {
	1002	if (octalSignificandWide >
	1003	(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
	1004	#ifndef NO_WIDE_TYPE
	1005	if (octalSignificandWide <= (MOST_BITS + signum)) {
	1006	objPtr->typePtr = &tclWideIntType;
	1007	if (signum) {
	1008	objPtr->internalRep.wideValue =
	1009	- (Tcl_WideInt) octalSignificandWide;
	1010	} else {
	1011	objPtr->internalRep.wideValue =
	1012	(Tcl_WideInt) octalSignificandWide;
	1013	}
	1014	break;
	1015	}
	1016	#endif
	1017	TclBNInitBignumFromWideUInt(&octalSignificandBig,
	1018	octalSignificandWide);
	1019	octalSignificandOverflow = 1;
	1020	} else {
	1021	objPtr->typePtr = &tclIntType;
	1022	if (signum) {
	1023	objPtr->internalRep.longValue =
	1024	- (long) octalSignificandWide;
	1025	} else {
	1026	objPtr->internalRep.longValue =
	1027	(long) octalSignificandWide;
	1028	}
	1029	}
	1030	}
	1031	if (octalSignificandOverflow) {
	1032	if (signum) {
	1033	mp_neg(&octalSignificandBig, &octalSignificandBig);
	1034	}
	1035	TclSetBignumIntRep(objPtr, &octalSignificandBig);
	1036	}
	1037	break;
	1038
	1039	case ZERO:
	1040	case DECIMAL:
	1041	significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
	1042	&significandWide, &significandBig, significandOverflow);
	1043	if (!significandOverflow && (significandWide > MOST_BITS+signum)) {
	1044	significandOverflow = 1;
	1045	TclBNInitBignumFromWideUInt(&significandBig, significandWide);
	1046	}
	1047	returnInteger:
	1048	if (!significandOverflow) {
	1049	if (significandWide >
	1050	(Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
	1051	#ifndef NO_WIDE_TYPE
	1052	if (significandWide <= MOST_BITS+signum) {
	1053	objPtr->typePtr = &tclWideIntType;
	1054	if (signum) {
	1055	objPtr->internalRep.wideValue =
	1056	- (Tcl_WideInt) significandWide;
	1057	} else {
	1058	objPtr->internalRep.wideValue =
	1059	(Tcl_WideInt) significandWide;
	1060	}
	1061	break;
	1062	}
	1063	#endif
	1064	TclBNInitBignumFromWideUInt(&significandBig,
	1065	significandWide);
	1066	significandOverflow = 1;
	1067	} else {
	1068	objPtr->typePtr = &tclIntType;
	1069	if (signum) {
	1070	objPtr->internalRep.longValue =
	1071	- (long) significandWide;
	1072	} else {
	1073	objPtr->internalRep.longValue =
	1074	(long) significandWide;
	1075	}
	1076	}
	1077	}
	1078	if (significandOverflow) {
	1079	if (signum) {
	1080	mp_neg(&significandBig, &significandBig);
	1081	}
	1082	TclSetBignumIntRep(objPtr, &significandBig);
	1083	}
	1084	break;
	1085
	1086	case FRACTION:
	1087	case EXPONENT:
	1088
	1089	/*
	1090	* Here, we're parsing a floating-point number. 'significandWide'
	1091	* or 'significandBig' contains the exact significand, according
	1092	* to whether 'significandOverflow' is set. The desired floating
	1093	* point value is significand * 10**k, where
	1094	* k = numTrailZeros+exponent-numDigitsAfterDp.
	1095	*/
	1096
	1097	objPtr->typePtr = &tclDoubleType;
	1098	if (exponentSignum) {
	1099	exponent = - exponent;
	1100	}
	1101	if (!significandOverflow) {
	1102	objPtr->internalRep.doubleValue = MakeLowPrecisionDouble(
	1103	signum, significandWide, numSigDigs,
	1104	(numTrailZeros + exponent - numDigitsAfterDp));
	1105	} else {
	1106	objPtr->internalRep.doubleValue = MakeHighPrecisionDouble(
	1107	signum, &significandBig, numSigDigs,
	1108	(numTrailZeros + exponent - numDigitsAfterDp));
	1109	}
	1110	break;
	1111
	1112	case sINF:
	1113	case sINFINITY:
	1114	if (signum) {
	1115	objPtr->internalRep.doubleValue = -HUGE_VAL;
	1116	} else {
	1117	objPtr->internalRep.doubleValue = HUGE_VAL;
	1118	}
	1119	objPtr->typePtr = &tclDoubleType;
	1120	break;
	1121
	1122	case sNAN:
	1123	case sNANFINISH:
	1124	objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide);
	1125	objPtr->typePtr = &tclDoubleType;
	1126	break;
	1127
	1128	case INITIAL:
	1129	/* This case only to silence compiler warning */
	1130	Tcl_Panic("TclParseNumber: state INITIAL can't happen here");
	1131	}
	1132	}
	1133
	1134	/*
	1135	* Format an error message when an invalid number is encountered.
	1136	*/
	1137
	1138	if (status != TCL_OK) {
	1139	if (interp != NULL) {
	1140	Tcl_Obj *msg;
	1141
	1142	TclNewLiteralStringObj(msg, "expected ");
	1143	Tcl_AppendToObj(msg, expected, -1);
	1144	Tcl_AppendToObj(msg, " but got \"", -1);
	1145	Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, "");
	1146	Tcl_AppendToObj(msg, "\"", -1);
	1147	if (state == BAD_OCTAL) {
	1148	Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1);
	1149	}
	1150	Tcl_SetObjResult(interp, msg);
	1151	}
	1152	}
	1153
	1154	/*
	1155	* Free memory.
	1156	*/
	1157
	1158	if (octalSignificandOverflow) {
	1159	mp_clear(&octalSignificandBig);
	1160	}
	1161	if (significandOverflow) {
	1162	mp_clear(&significandBig);
	1163	}
	1164	return status;
	1165	}
	1166
	1167	/*
	1168	*----------------------------------------------------------------------
	1169	*
	1170	* AccumulateDecimalDigit --
	1171	*
	1172	* Consume a decimal digit in a number being scanned.
	1173	*
	1174	* Results:
	1175	* Returns 1 if the number has overflowed to a bignum, 0 if it still fits
	1176	* in a wide integer.
	1177	*
	1178	* Side effects:
	1179	* Updates either the wide or bignum representation.
	1180	*
	1181	*----------------------------------------------------------------------
	1182	*/
	1183
	1184	static int
	1185	AccumulateDecimalDigit(
	1186	unsigned digit, /* Digit being scanned. */
	1187	int numZeros, /* Count of zero digits preceding the digit
	1188	* being scanned. */
	1189	Tcl_WideUInt wideRepPtr, / Representation of the partial number as a
	1190	* wide integer. */
	1191	mp_int bignumRepPtr, / Representation of the partial number as a
	1192	* bignum. */
	1193	int bignumFlag) /* Flag == 1 if the number overflowed previous
	1194	* to this digit. */
	1195	{
	1196	int i, n;
	1197	Tcl_WideUInt w;
	1198
	1199	/*
	1200	* Try wide multiplication first
	1201	*/
	1202
	1203	if (!bignumFlag) {
	1204	w = *wideRepPtr;
	1205	if (w == 0) {
	1206	/*
	1207	* There's no need to multiply if the multiplicand is zero.
	1208	*/
	1209
	1210	*wideRepPtr = digit;
	1211	return 0;
	1212	} else if (numZeros >= maxpow10_wide
	1213	\|\| w > ((~(Tcl_WideUInt)0)-digit)/pow10_wide[numZeros+1]) {
	1214	/*
	1215	* Wide multiplication will overflow. Expand the
	1216	* number to a bignum and fall through into the bignum case.
	1217	*/
	1218
	1219	TclBNInitBignumFromWideUInt (bignumRepPtr, w);
	1220	} else {
	1221	/*
	1222	* Wide multiplication.
	1223	*/
	1224	wideRepPtr = w pow10_wide[numZeros+1] + digit;
	1225	return 0;
	1226	}
	1227	}
	1228
	1229	/*
	1230	* Bignum multiplication.
	1231	*/
	1232
	1233	if (numZeros < log10_DIGIT_MAX) {
	1234	/*
	1235	* Up to about 8 zeros - single digit multiplication.
	1236	*/
	1237
	1238	mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
	1239	bignumRepPtr);
	1240	mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
	1241	} else {
	1242	/*
	1243	* More than single digit multiplication. Multiply by the appropriate
	1244	* small powers of 5, and then shift. Large strings of zeroes are
	1245	* eaten 256 at a time; this is less efficient than it could be, but
	1246	* seems implausible. We presume that DIGIT_BIT is at least 27. The
	1247	* first multiplication, by up to 10**7, is done with a one-DIGIT
	1248	* multiply (this presumes that DIGIT_BIT >= 24).
	1249	*/
	1250
	1251	n = numZeros + 1;
	1252	mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
	1253	for (i=3; i<=7; ++i) {
	1254	if (n & (1 << i)) {
	1255	mp_mul(bignumRepPtr, pow5+i, bignumRepPtr);
	1256	}
	1257	}
	1258	while (n >= 256) {
	1259	mp_mul(bignumRepPtr, pow5+8, bignumRepPtr);
	1260	n -= 256;
	1261	}
	1262	mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr);
	1263	mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
	1264	}
	1265
	1266	return 1;
	1267	}
	1268
	1269	/*
	1270	*----------------------------------------------------------------------
	1271	*
	1272	* MakeLowPrecisionDouble --
	1273	*
	1274	* Makes the double precision number, signumsignificand10**exponent.
	1275	*
	1276	* Results:
	1277	* Returns the constructed number.
	1278	*
	1279	* Common cases, where there are few enough digits that the number can be
	1280	* represented with at most roundoff, are handled specially here. If the
	1281	* number requires more than one rounded operation to compute, the code
	1282	* promotes the significand to a bignum and calls MakeHighPrecisionDouble
	1283	* to do it instead.
	1284	*
	1285	*----------------------------------------------------------------------
	1286	*/
	1287
	1288	static double
	1289	MakeLowPrecisionDouble(
	1290	int signum, /* 1 if the number is negative, 0 otherwise */
	1291	Tcl_WideUInt significand, /* Significand of the number */
	1292	int numSigDigs, /* Number of digits in the significand */
	1293	int exponent) /* Power of ten */
	1294	{
	1295	double retval; /* Value of the number */
	1296	mp_int significandBig; /* Significand expressed as a bignum */
	1297
	1298	/*
	1299	* With gcc on x86, the floating point rounding mode is double-extended.
	1300	* This causes the result of double-precision calculations to be rounded
	1301	* twice: once to the precision of double-extended and then again to the
	1302	* precision of double. Double-rounding introduces gratuitous errors of 1
	1303	* ulp, so we need to change rounding mode to 53-bits.
	1304	*/
	1305
	1306	#if defined(__GNUC__) && defined(__i386)
	1307	fpu_control_t roundTo53Bits = 0x027f;
	1308	fpu_control_t oldRoundingMode;
	1309	_FPU_GETCW(oldRoundingMode);
	1310	_FPU_SETCW(roundTo53Bits);
	1311	#endif
	1312
	1313	/*
	1314	* Test for the easy cases.
	1315	*/
	1316
	1317	if (numSigDigs <= DBL_DIG) {
	1318	if (exponent >= 0) {
	1319	if (exponent <= mmaxpow) {
	1320	/*
	1321	* The significand is an exact integer, and so is
	1322	* 10**exponent. The product will be correct to within 1/2 ulp
	1323	* without special handling.
	1324	*/
	1325
	1326	retval = (double)(Tcl_WideInt)significand * pow10vals[ exponent ];
	1327	goto returnValue;
	1328	} else {
	1329	int diff = DBL_DIG - numSigDigs;
	1330	if (exponent-diff <= mmaxpow) {
	1331	/*
	1332	* 10**exponent is not an exact integer, but
	1333	* 10**(exponent-diff) is exact, and so is
	1334	* significand10*diff, so we can still compute the value
	1335	* with only one roundoff.
	1336	*/
	1337
	1338	volatile double factor =
	1339	(double)(Tcl_WideInt)significand * pow10vals[diff];
	1340	retval = factor * pow10vals[exponent-diff];
	1341	goto returnValue;
	1342	}
	1343	}
	1344	} else {
	1345	if (exponent >= -mmaxpow) {
	1346	/*
	1347	* 10**-exponent is an exact integer, and so is the
	1348	* significand. Compute the result by one division, again with
	1349	* only one rounding.
	1350	*/
	1351
	1352	retval = (double)(Tcl_WideInt)significand / pow10vals[-exponent];
	1353	goto returnValue;
	1354	}
	1355	}
	1356	}
	1357
	1358	/*
	1359	* All the easy cases have failed. Promote ths significand to bignum and
	1360	* call MakeHighPrecisionDouble to do it the hard way.
	1361	*/
	1362
	1363	TclBNInitBignumFromWideUInt(&significandBig, significand);
	1364	retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
	1365	exponent);
	1366	mp_clear(&significandBig);
	1367
	1368	/*
	1369	* Come here to return the computed value.
	1370	*/
	1371
	1372	returnValue:
	1373	if (signum) {
	1374	retval = -retval;
	1375	}
	1376
	1377	/*
	1378	* On gcc on x86, restore the floating point mode word.
	1379	*/
	1380
	1381	#if defined(__GNUC__) && defined(__i386)
	1382	_FPU_SETCW(oldRoundingMode);
	1383	#endif
	1384
	1385	return retval;
	1386	}
	1387
	1388	/*
	1389	*----------------------------------------------------------------------
	1390	*
	1391	* MakeHighPrecisionDouble --
	1392	*
	1393	* Makes the double precision number, signumsignificand10**exponent.
	1394	*
	1395	* Results:
	1396	* Returns the constructed number.
	1397	*
	1398	* MakeHighPrecisionDouble is used when arbitrary-precision arithmetic is
	1399	* needed to ensure correct rounding. It begins by calculating a
	1400	* low-precision approximation to the desired number, and then refines
	1401	* the answer in high precision.
	1402	*
	1403	*----------------------------------------------------------------------
	1404	*/
	1405
	1406	static double
	1407	MakeHighPrecisionDouble(
	1408	int signum, /* 1=negative, 0=nonnegative */
	1409	mp_int significand, / Exact significand of the number */
	1410	int numSigDigs, /* Number of significant digits */
	1411	int exponent) /* Power of 10 by which to multiply */
	1412	{
	1413	double retval;
	1414	int machexp; /* Machine exponent of a power of 10 */
	1415
	1416	/*
	1417	* With gcc on x86, the floating point rounding mode is double-extended.
	1418	* This causes the result of double-precision calculations to be rounded
	1419	* twice: once to the precision of double-extended and then again to the
	1420	* precision of double. Double-rounding introduces gratuitous errors of 1
	1421	* ulp, so we need to change rounding mode to 53-bits.
	1422	*/
	1423
	1424	#if defined(__GNUC__) && defined(__i386)
	1425	fpu_control_t roundTo53Bits = 0x027f;
	1426	fpu_control_t oldRoundingMode;
	1427	_FPU_GETCW(oldRoundingMode);
	1428	_FPU_SETCW(roundTo53Bits);
	1429	#endif
	1430
	1431	/*
	1432	* Quick checks for over/underflow.
	1433	*/
	1434
	1435	if (numSigDigs+exponent-1 > maxDigits) {
	1436	retval = HUGE_VAL;
	1437	goto returnValue;
	1438	}
	1439	if (numSigDigs+exponent-1 < minDigits) {
	1440	retval = 0;
	1441	goto returnValue;
	1442	}
	1443
	1444	/*
	1445	* Develop a first approximation to the significand. It is tempting simply
	1446	* to force bignum to double, but that will overflow on input numbers like
	1447	* 1.[string repeat 0 1000]1; while this is a not terribly likely
	1448	* scenario, we still have to deal with it. Use fraction and exponent
	1449	* instead. Once we have the significand, multiply by 10**exponent. Test
	1450	* for overflow. Convert back to a double, and test for underflow.
	1451	*/
	1452
	1453	retval = BignumToBiasedFrExp(significand, &machexp);
	1454	retval = Pow10TimesFrExp(exponent, retval, &machexp);
	1455	if (machexp > DBL_MAX_EXP*log2FLT_RADIX) {
	1456	retval = HUGE_VAL;
	1457	goto returnValue;
	1458	}
	1459	retval = SafeLdExp(retval, machexp);
	1460	if (retval < tiny) {
	1461	retval = tiny;
	1462	}
	1463
	1464	/*
	1465	* Refine the result twice. (The second refinement should be necessary
	1466	* only if the best approximation is a power of 2 minus 1/2 ulp).
	1467	*/
	1468
	1469	retval = RefineApproximation(retval, significand, exponent);
	1470	retval = RefineApproximation(retval, significand, exponent);
	1471
	1472	/*
	1473	* Come here to return the computed value.
	1474	*/
	1475
	1476	returnValue:
	1477	if (signum) {
	1478	retval = -retval;
	1479	}
	1480
	1481	/*
	1482	* On gcc on x86, restore the floating point mode word.
	1483	*/
	1484
	1485	#if defined(__GNUC__) && defined(__i386)
	1486	_FPU_SETCW(oldRoundingMode);
	1487	#endif
	1488	return retval;
	1489	}
	1490
	1491	/*
	1492	*----------------------------------------------------------------------
	1493	*
	1494	* MakeNaN --
	1495	*
	1496	* Makes a "Not a Number" given a set of bits to put in the tag bits
	1497	*
	1498	* Note that a signalling NaN is never returned.
	1499	*
	1500	*----------------------------------------------------------------------
	1501	*/
	1502
	1503	#ifdef IEEE_FLOATING_POINT
	1504	static double
	1505	MakeNaN(
	1506	int signum, /* Sign bit (1=negative, 0=nonnegative */
	1507	Tcl_WideUInt tags) /* Tag bits to put in the NaN */
	1508	{
	1509	union {
	1510	Tcl_WideUInt iv;
	1511	double dv;
	1512	} theNaN;
	1513
	1514	theNaN.iv = tags;
	1515	theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
	1516	if (signum) {
	1517	theNaN.iv \|= ((Tcl_WideUInt) (0x8000 \| NAN_START)) << 48;
	1518	} else {
	1519	theNaN.iv \|= ((Tcl_WideUInt) NAN_START) << 48;
	1520	}
	1521	if (n770_fp) {
	1522	theNaN.iv = Nokia770Twiddle(theNaN.iv);
	1523	}
	1524	return theNaN.dv;
	1525	}
	1526	#endif
	1527
	1528	/*
	1529	*----------------------------------------------------------------------
	1530	*
	1531	* RefineApproximation --
	1532	*
	1533	* Given a poor approximation to a floating point number, returns a
	1534	* better one. (The better approximation is correct to within 1 ulp, and
	1535	* is entirely correct if the poor approximation is correct to 1 ulp.)
	1536	*
	1537	* Results:
	1538	* Returns the improved result.
	1539	*
	1540	*----------------------------------------------------------------------
	1541	*/
	1542
	1543	static double
	1544	RefineApproximation(
	1545	double approxResult, /* Approximate result of conversion */
	1546	mp_int exactSignificand, / Integer significand */
	1547	int exponent) /* Power of 10 to multiply by significand */
	1548	{
	1549	int M2, M5; /* Powers of 2 and of 5 needed to put the
	1550	* decimal and binary numbers over a common
	1551	* denominator. */
	1552	double significand; /* Sigificand of the binary number */
	1553	int binExponent; /* Exponent of the binary number */
	1554	int msb; /* Most significant bit position of an
	1555	* intermediate result */
	1556	int nDigits; /* Number of mp_digit's in an intermediate
	1557	* result */
	1558	mp_int twoMv; /* Approx binary value expressed as an exact
	1559	* integer scaled by the multiplier 2M */
	1560	mp_int twoMd; /* Exact decimal value expressed as an exact
	1561	* integer scaled by the multiplier 2M */
	1562	int scale; /* Scale factor for M */
	1563	int multiplier; /* Power of two to scale M */
	1564	double num, den; /* Numerator and denominator of the correction
	1565	* term */
	1566	double quot; /* Correction term */
	1567	double minincr; /* Lower bound on the absolute value of the
	1568	* correction term. */
	1569	int i;
	1570
	1571	/*
	1572	* The first approximation is always low. If we find that it's HUGE_VAL,
	1573	* we're done.
	1574	*/
	1575
	1576	if (approxResult == HUGE_VAL) {
	1577	return approxResult;
	1578	}
	1579
	1580	/*
	1581	* Find a common denominator for the decimal and binary fractions. The
	1582	* common denominator will be 2M2 + 5M5.
	1583	*/
	1584
	1585	significand = frexp(approxResult, &binExponent);
	1586	i = mantBits - binExponent;
	1587	if (i < 0) {
	1588	M2 = 0;
	1589	} else {
	1590	M2 = i;
	1591	}
	1592	if (exponent > 0) {
	1593	M5 = 0;
	1594	} else {
	1595	M5 = -exponent;
	1596	if ((M5-1) > M2) {
	1597	M2 = M5-1;
	1598	}
	1599	}
	1600
	1601	/*
	1602	* The floating point number is significand2*binExponent. Compute the
	1603	* large integer significand2(binExponent+M2+1). The 2*-1 bit of the
	1604	* significand (the most significant) corresponds to the
	1605	* 2*(binExponent+M2 + 1) bit of 2M2*v. Allocate enough digits to hold
	1606	* that quantity, then convert the significand to a large integer, scaled
	1607	* appropriately. Then multiply by the appropriate power of 5.
	1608	*/
	1609
	1610	msb = binExponent + M2; /* 1008 */
	1611	nDigits = msb / DIGIT_BIT + 1;
	1612	mp_init_size(&twoMv, nDigits);
	1613	i = (msb % DIGIT_BIT + 1);
	1614	twoMv.used = nDigits;
	1615	significand *= SafeLdExp(1.0, i);
	1616	while (--nDigits >= 0) {
	1617	twoMv.dp[nDigits] = (mp_digit) significand;
	1618	significand -= (mp_digit) significand;
	1619	significand = SafeLdExp(significand, DIGIT_BIT);
	1620	}
	1621	for (i = 0; i <= 8; ++i) {
	1622	if (M5 & (1 << i)) {
	1623	mp_mul(&twoMv, pow5+i, &twoMv);
	1624	}
	1625	}
	1626
	1627	/*
	1628	* Collect the decimal significand as a high precision integer. The least
	1629	* significant bit corresponds to bit M2+exponent+1 so it will need to be
	1630	* shifted left by that many bits after being multiplied by
	1631	* 5**(M5+exponent).
	1632	*/
	1633
	1634	mp_init_copy(&twoMd, exactSignificand);
	1635	for (i=0; i<=8; ++i) {
	1636	if ((M5+exponent) & (1 << i)) {
	1637	mp_mul(&twoMd, pow5+i, &twoMd);
	1638	}
	1639	}
	1640	mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
	1641	mp_sub(&twoMd, &twoMv, &twoMd);
	1642
	1643	/*
	1644	* The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
	1645	* term. Because 2M may well overflow a double, we need to scale the
	1646	* denominator by a factor of 2**binExponent-mantBits
	1647	*/
	1648
	1649	scale = binExponent - mantBits - 1;
	1650
	1651	mp_set(&twoMv, 1);
	1652	for (i=0; i<=8; ++i) {
	1653	if (M5 & (1 << i)) {
	1654	mp_mul(&twoMv, pow5+i, &twoMv);
	1655	}
	1656	}
	1657	multiplier = M2 + scale + 1;
	1658	if (multiplier > 0) {
	1659	mp_mul_2d(&twoMv, multiplier, &twoMv);
	1660	} else if (multiplier < 0) {
	1661	mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
	1662	}
	1663
	1664	/*
	1665	* If the result is less than unity, the error is less than 1/2 unit in
	1666	* the last place, so there's no correction to make.
	1667	*/
	1668
	1669	if (mp_cmp_mag(&twoMd, &twoMv) == MP_LT) {
	1670	mp_clear(&twoMd);
	1671	mp_clear(&twoMv);
	1672	return approxResult;
	1673	}
	1674
	1675	/*
	1676	* Convert the numerator and denominator of the corrector term accurately
	1677	* to floating point numbers.
	1678	*/
	1679
	1680	num = TclBignumToDouble(&twoMd);
	1681	den = TclBignumToDouble(&twoMv);
	1682
	1683	quot = SafeLdExp(num/den, scale);
	1684	minincr = SafeLdExp(1.0, binExponent-mantBits);
	1685
	1686	if (quot<0. && quot>-minincr) {
	1687	quot = -minincr;
	1688	} else if (quot>0. && quot<minincr) {
	1689	quot = minincr;
	1690	}
	1691
	1692	mp_clear(&twoMd);
	1693	mp_clear(&twoMv);
	1694
	1695	return approxResult + quot;
	1696	}
	1697
	1698	/*
	1699	*----------------------------------------------------------------------
	1700	*
	1701	* TclDoubleDigits --
	1702	*
	1703	* Converts a double to a string of digits.
	1704	*
	1705	* Results:
	1706	* Returns the position of the character in the string after which the
	1707	* decimal point should appear. Since the string contains only
	1708	* significant digits, the position may be less than zero or greater than
	1709	* the length of the string.
	1710	*
	1711	* Side effects:
	1712	* Stores the digits in the given buffer and sets 'signum' according to
	1713	* the sign of the number.
	1714	*
	1715	*----------------------------------------------------------------------
	1716
	1717	*/
	1718
	1719	int
	1720	TclDoubleDigits(
	1721	char buffer, / Buffer in which to store the result, must
	1722	* have at least 18 chars */
	1723	double v, /* Number to convert. Must be finite, and not
	1724	* NaN */
	1725	int signum) / Output: 1 if the number is negative.
	1726	* Should handle -0 correctly on the IEEE
	1727	* architecture. */
	1728	{
	1729	int e; /* Power of FLT_RADIX that satisfies
	1730	* v = f * FLT_RADIX*e /
	1731	int lowOK, highOK;
	1732	mp_int r; /* Scaled significand. */
	1733	mp_int s; /* Divisor such that v = r / s */
	1734	int smallestSig; /* Flag == 1 iff v's significand is the
	1735	* smallest that can be represented. */
	1736	mp_int mplus; /* Scaled epsilon: (r + 2* mplus) == v(+)
	1737	* where v(+) is the floating point successor
	1738	* of v. */
	1739	mp_int mminus; /* Scaled epsilon: (r - 2*mminus) == v(-)
	1740	* where v(-) is the floating point
	1741	* predecessor of v. */
	1742	mp_int temp;
	1743	int rfac2 = 0; /* Powers of 2 and 5 by which large */
	1744	int rfac5 = 0; /* integers should be scaled. */
	1745	int sfac2 = 0;
	1746	int sfac5 = 0;
	1747	int mplusfac2 = 0;
	1748	int mminusfac2 = 0;
	1749	char c;
	1750	int i, k, n;
	1751
	1752	/*
	1753	* Split the number into absolute value and signum.
	1754	*/
	1755
	1756	v = AbsoluteValue(v, signum);
	1757
	1758	/*
	1759	* Handle zero specially.
	1760	*/
	1761
	1762	if (v == 0.0) {
	1763	*buffer++ = '0';
	1764	*buffer++ = '\0';
	1765	return 1;
	1766	}
	1767
	1768	/*
	1769	* Find a large integer r, and integer e, such that
	1770	* v = r * FLT_RADIX**e
	1771	* and r is as small as possible. Also determine whether the significand
	1772	* is the smallest possible.
	1773	*/
	1774
	1775	smallestSig = GetIntegerTimesPower(v, &r, &e);
	1776
	1777	lowOK = highOK = (mp_iseven(&r));
	1778
	1779	/*
	1780	* We are going to want to develop integers r, s, mplus, and mminus such
	1781	* that v = r / s, v(+)-v / 2 = mplus / s; v-v(-) / 2 = mminus / s and
	1782	* then scale either s or r, mplus, mminus by an appropriate power of ten.
	1783	*
	1784	* We actually do this by keeping track of the powers of 2 and 5 by which
	1785	* f is multiplied to yield v and by which 1 is multiplied to yield s,
	1786	* mplus, and mminus.
	1787	*/
	1788
	1789	if (e >= 0) {
	1790	int bits = e * log2FLT_RADIX;
	1791
	1792	if (!smallestSig) {
	1793	/*
	1794	* Normal case, m+ and m- are both FLT_RADIX**e
	1795	*/
	1796
	1797	rfac2 = bits + 1;
	1798	sfac2 = 1;
	1799	mplusfac2 = bits;
	1800	mminusfac2 = bits;
	1801	} else {
	1802	/*
	1803	* If f is equal to the smallest significand, then we need another
	1804	* factor of FLT_RADIX in s to cope with stepping to the next
	1805	* smaller exponent when going to e's predecessor.
	1806	*/
	1807
	1808	rfac2 = bits + log2FLT_RADIX + 1;
	1809	sfac2 = 1 + log2FLT_RADIX;
	1810	mplusfac2 = bits + log2FLT_RADIX;
	1811	mminusfac2 = bits;
	1812	}
	1813	} else {
	1814	/*
	1815	* v has digits after the binary point
	1816	*/
	1817
	1818	if (e <= DBL_MIN_EXP-DBL_MANT_DIG \|\| !smallestSig) {
	1819	/*
	1820	* Either f isn't the smallest significand or e is the smallest
	1821	* exponent. mplus and mminus will both be 1.
	1822	*/
	1823
	1824	rfac2 = 1;
	1825	sfac2 = 1 - e * log2FLT_RADIX;
	1826	mplusfac2 = 0;
	1827	mminusfac2 = 0;
	1828	} else {
	1829	/*
	1830	* f is the smallest significand, but e is not the smallest
	1831	* exponent. We need to scale by FLT_RADIX again to cope with the
	1832	* fact that v's predecessor has a smaller exponent.
	1833	*/
	1834
	1835	rfac2 = 1 + log2FLT_RADIX;
	1836	sfac2 = 1 + log2FLT_RADIX * (1 - e);
	1837	mplusfac2 = FLT_RADIX;
	1838	mminusfac2 = 0;
	1839	}
	1840	}
	1841
	1842	/*
	1843	* Estimate the highest power of ten that will be needed to hold the
	1844	* result.
	1845	*/
	1846
	1847	k = (int) ceil(log(v) / log(10.));
	1848	if (k >= 0) {
	1849	sfac2 += k;
	1850	sfac5 = k;
	1851	} else {
	1852	rfac2 -= k;
	1853	mplusfac2 -= k;
	1854	mminusfac2 -= k;
	1855	rfac5 = -k;
	1856	}
	1857
	1858	/*
	1859	* Scale r, s, mplus, mminus by the appropriate powers of 2 and 5.
	1860	*/
	1861
	1862	mp_init_set(&mplus, 1);
	1863	for (i=0 ; i<=8 ; ++i) {
	1864	if (rfac5 & (1 << i)) {
	1865	mp_mul(&mplus, pow5+i, &mplus);
	1866	}
	1867	}
	1868	mp_mul(&r, &mplus, &r);
	1869	mp_mul_2d(&r, rfac2, &r);
	1870	mp_init_copy(&mminus, &mplus);
	1871	mp_mul_2d(&mplus, mplusfac2, &mplus);
	1872	mp_mul_2d(&mminus, mminusfac2, &mminus);
	1873	mp_init_set(&s, 1);
	1874	for (i=0 ; i<=8 ; ++i) {
	1875	if (sfac5 & (1 << i)) {
	1876	mp_mul(&s, pow5+i, &s);
	1877	}
	1878	}
	1879	mp_mul_2d(&s, sfac2, &s);
	1880
	1881	/*
	1882	* It is possible for k to be off by one because we used an inexact
	1883	* logarithm.
	1884	*/
	1885
	1886	mp_init(&temp);
	1887	mp_add(&r, &mplus, &temp);
	1888	i = mp_cmp_mag(&temp, &s);
	1889	if (i>0 \|\| (highOK && i==0)) {
	1890	mp_mul_d(&s, 10, &s);
	1891	++k;
	1892	} else {
	1893	mp_mul_d(&temp, 10, &temp);
	1894	i = mp_cmp_mag(&temp, &s);
	1895	if (i<0 \|\| (highOK && i==0)) {
	1896	mp_mul_d(&r, 10, &r);
	1897	mp_mul_d(&mplus, 10, &mplus);
	1898	mp_mul_d(&mminus, 10, &mminus);
	1899	--k;
	1900	}
	1901	}
	1902
	1903	/*
	1904	* At this point, k contains the power of ten by which we're scaling the
	1905	* result. r/s is at least 1/10 and strictly less than ten, and v = r/s *
	1906	* 10**k. mplus and mminus give the rounding limits.
	1907	*/
	1908
	1909	for (;;) {
	1910	int tc1, tc2;
	1911
	1912	mp_mul_d(&r, 10, &r);
	1913	mp_div(&r, &s, &temp, &r); /* temp = 10r / s; r = 10r mod s */
	1914	i = temp.dp[0];
	1915	mp_mul_d(&mplus, 10, &mplus);
	1916	mp_mul_d(&mminus, 10, &mminus);
	1917	tc1 = mp_cmp_mag(&r, &mminus);
	1918	if (lowOK) {
	1919	tc1 = (tc1 <= 0);
	1920	} else {
	1921	tc1 = (tc1 < 0);
	1922	}
	1923	mp_add(&r, &mplus, &temp);
	1924	tc2 = mp_cmp_mag(&temp, &s);
	1925	if (highOK) {
	1926	tc2 = (tc2 >= 0);
	1927	} else {
	1928	tc2= (tc2 > 0);
	1929	}
	1930	if (!tc1) {
	1931	if (!tc2) {
	1932	*buffer++ = '0' + i;
	1933	} else {
	1934	c = (char) (i + '1');
	1935	break;
	1936	}
	1937	} else {
	1938	if (!tc2) {
	1939	c = (char) (i + '0');
	1940	} else {
	1941	mp_mul_2d(&r, 1, &r);
	1942	n = mp_cmp_mag(&r, &s);
	1943	if (n < 0) {
	1944	c = (char) (i + '0');
	1945	} else {
	1946	c = (char) (i + '1');
	1947	}
	1948	}
	1949	break;
	1950	}
	1951	};
	1952	*buffer++ = c;
	1953	*buffer++ = '\0';
	1954
	1955	/*
	1956	* Free memory, and return.
	1957	*/
	1958
	1959	mp_clear_multi(&r, &s, &mplus, &mminus, &temp, NULL);
	1960	return k;
	1961	}
	1962
	1963	/*
	1964	*----------------------------------------------------------------------
	1965	*
	1966	* AbsoluteValue --
	1967	*
	1968	* Splits a 'double' into its absolute value and sign.
	1969	*
	1970	* Results:
	1971	* Returns the absolute value.
	1972	*
	1973	* Side effects:
	1974	* Stores the signum in '*signum'.
	1975	*
	1976	*----------------------------------------------------------------------
	1977	*/
	1978
	1979	static double
	1980	AbsoluteValue(
	1981	double v, /* Number to split */
	1982	int signum) / (Output) Sign of the number 1=-, 0=+ */
	1983	{
	1984	/*
	1985	* Take the absolute value of the number, and report the number's sign.
	1986	* Take special steps to preserve signed zeroes in IEEE floating point.
	1987	* (We can't use fpclassify, because that's a C9x feature and we still
	1988	* have to build on C89 compilers.)
	1989	*/
	1990
	1991	#ifndef IEEE_FLOATING_POINT
	1992	if (v >= 0.0) {
	1993	*signum = 0;
	1994	} else {
	1995	*signum = 1;
	1996	v = -v;
	1997	}
	1998	#else
	1999	union {
	2000	Tcl_WideUInt iv;
	2001	double dv;
	2002	} bitwhack;
	2003	bitwhack.dv = v;
	2004	if (n770_fp) {
	2005	bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
	2006	}
	2007	if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
	2008	*signum = 1;
	2009	bitwhack.iv &= ~((Tcl_WideUInt) 1 << 63);
	2010	if (n770_fp) {
	2011	bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
	2012	}
	2013	v = bitwhack.dv;
	2014	} else {
	2015	*signum = 0;
	2016	}
	2017	#endif
	2018	return v;
	2019	}
	2020
	2021	/*
	2022	*----------------------------------------------------------------------
	2023	*
	2024	* GetIntegerTimesPower --
	2025	*
	2026	* Converts a floating point number to an exact integer times a power of
	2027	* the floating point radix.
	2028	*
	2029	* Results:
	2030	* Returns 1 if it converted the smallest significand, 0 otherwise.
	2031	*
	2032	* Side effects:
	2033	* Initializes the integer value (does not just assign it), and stores
	2034	* the exponent.
	2035	*
	2036	*----------------------------------------------------------------------
	2037	*/
	2038
	2039	static int
	2040	GetIntegerTimesPower(
	2041	double v, /* Value to convert */
	2042	mp_int rPtr, / (Output) Integer value */
	2043	int ePtr) / (Output) Power of FLT_RADIX by which r must
	2044	* be multiplied to yield v*/
	2045	{
	2046	double a, f;
	2047	int e, i, n;
	2048
	2049	/*
	2050	* Develop f and e such that v = f * FLT_RADIX**e, with
	2051	* 1.0/FLT_RADIX <= f < 1.
	2052	*/
	2053
	2054	f = frexp(v, &e);
	2055	#if FLT_RADIX > 2
	2056	n = e % log2FLT_RADIX;
	2057	if (n > 0) {
	2058	n -= log2FLT_RADIX;
	2059	e += 1;
	2060	f *= ldexp(1.0, n);
	2061	}
	2062	e = (e - n) / log2FLT_RADIX;
	2063	#endif
	2064	if (f == 1.0) {
	2065	f = 1.0 / FLT_RADIX;
	2066	e += 1;
	2067	}
	2068
	2069	/*
	2070	* If the original number was denormalized, adjust e and f to be denormal
	2071	* as well.
	2072	*/
	2073
	2074	if (e < DBL_MIN_EXP) {
	2075	n = mantBits + (e - DBL_MIN_EXP)*log2FLT_RADIX;
	2076	f = ldexp(f, (e - DBL_MIN_EXP)*log2FLT_RADIX);
	2077	e = DBL_MIN_EXP;
	2078	n = (n + DIGIT_BIT - 1) / DIGIT_BIT;
	2079	} else {
	2080	n = mantDIGIT;
	2081	}
	2082
	2083	/*
	2084	* Now extract the base-2**DIGIT_BIT digits of f into a multi-precision
	2085	* integer r. Preserve the invariant v = r * 2*rfac2 FLT_RADIX**e by
	2086	* adjusting e.
	2087	*/
	2088
	2089	a = f;
	2090	n = mantDIGIT;
	2091	mp_init_size(rPtr, n);
	2092	rPtr->used = n;
	2093	rPtr->sign = MP_ZPOS;
	2094	i = (mantBits % DIGIT_BIT);
	2095	if (i == 0) {
	2096	i = DIGIT_BIT;
	2097	}
	2098	while (n > 0) {
	2099	a *= ldexp(1.0, i);
	2100	i = DIGIT_BIT;
	2101	rPtr->dp[--n] = (mp_digit) a;
	2102	a -= (mp_digit) a;
	2103	}
	2104	*ePtr = e - DBL_MANT_DIG;
	2105	return (f == 1.0 / FLT_RADIX);
	2106	}
	2107
	2108	/*
	2109	*----------------------------------------------------------------------
	2110	*
	2111	* TclInitDoubleConversion --
	2112	*
	2113	* Initializes constants that are needed for conversions to and from
	2114	* 'double'
	2115	*
	2116	* Results:
	2117	* None.
	2118	*
	2119	* Side effects:
	2120	* The log base 2 of the floating point radix, the number of bits in a
	2121	* double mantissa, and a table of the powers of five and ten are
	2122	* computed and stored.
	2123	*
	2124	*----------------------------------------------------------------------
	2125	*/
	2126
	2127	void
	2128	TclInitDoubleConversion(void)
	2129	{
	2130	int i;
	2131	int x;
	2132	Tcl_WideUInt u;
	2133	double d;
	2134
	2135	#ifdef IEEE_FLOATING_POINT
	2136	union {
	2137	double dv;
	2138	Tcl_WideUInt iv;
	2139	} bitwhack;
	2140	#endif
	2141
	2142	/*
	2143	* Initialize table of powers of 10 expressed as wide integers.
	2144	*/
	2145
	2146	maxpow10_wide = (int)
	2147	floor(sizeof(Tcl_WideUInt) * CHAR_BIT * log(2.) / log(10.));
	2148	pow10_wide = (Tcl_WideUInt *)
	2149	ckalloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
	2150	u = 1;
	2151	for (i = 0; i < maxpow10_wide; ++i) {
	2152	pow10_wide[i] = u;
	2153	u *= 10;
	2154	}
	2155	pow10_wide[i] = u;
	2156
	2157	/*
	2158	* Determine how many bits of precision a double has, and how many
	2159	* decimal digits that represents.
	2160	*/
	2161
	2162	if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
	2163	Tcl_Panic("This code doesn't work on a decimal machine!");
	2164	}
	2165	--log2FLT_RADIX;
	2166	mantBits = DBL_MANT_DIG * log2FLT_RADIX;
	2167	d = 1.0;
	2168
	2169	/*
	2170	* Initialize a table of powers of ten that can be exactly represented
	2171	* in a double.
	2172	*/
	2173
	2174	x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
	2175	if (x < MAXPOW) {
	2176	mmaxpow = x;
	2177	} else {
	2178	mmaxpow = MAXPOW;
	2179	}
	2180	for (i=0 ; i<=mmaxpow ; ++i) {
	2181	pow10vals[i] = d;
	2182	d *= 10.0;
	2183	}
	2184
	2185	/*
	2186	* Initialize a table of large powers of five.
	2187	*/
	2188
	2189	for (i=0; i<9; ++i) {
	2190	mp_init(pow5 + i);
	2191	}
	2192	mp_set(pow5, 5);
	2193	for (i=0; i<8; ++i) {
	2194	mp_sqr(pow5+i, pow5+i+1);
	2195	}
	2196
	2197	/*
	2198	* Determine the number of decimal digits to the left and right of the
	2199	* decimal point in the largest and smallest double, the smallest double
	2200	* that differs from zero, and the number of mp_digits needed to represent
	2201	* the significand of a double.
	2202	*/
	2203
	2204	tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
	2205	maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX)
	2206	+ 0.5 * log(10.)) / log(10.));
	2207	minDigits = (int) floor((DBL_MIN_EXP - DBL_MANT_DIG)
	2208	* log((double) FLT_RADIX) / log(10.));
	2209	mantDIGIT = (mantBits + DIGIT_BIT-1) / DIGIT_BIT;
	2210	log10_DIGIT_MAX = (int) floor(DIGIT_BIT * log(2.) / log(10.));
	2211
	2212	/*
	2213	* Nokia 770's software-emulated floating point is "middle endian": the
	2214	* bytes within a 32-bit word are little-endian (like the native
	2215	* integers), but the two words of a 'double' are presented most
	2216	* significant word first.
	2217	*/
	2218
	2219	#ifdef IEEE_FLOATING_POINT
	2220	bitwhack.dv = 1.000000238418579;
	2221	/* 3ff0 0000 4000 0000 */
	2222	if ((bitwhack.iv >> 32) == 0x3ff00000) {
	2223	n770_fp = 0;
	2224	} else if ((bitwhack.iv & 0xffffffff) == 0x3ff00000) {
	2225	n770_fp = 1;
	2226	} else {
	2227	Tcl_Panic("unknown floating point word order on this machine");
	2228	}
	2229	#endif
	2230	}
	2231
	2232	/*
	2233	*----------------------------------------------------------------------
	2234	*
	2235	* TclFinalizeDoubleConversion --
	2236	*
	2237	* Cleans up this file on exit.
	2238	*
	2239	* Results:
	2240	* None
	2241	*
	2242	* Side effects:
	2243	* Memory allocated by TclInitDoubleConversion is freed.
	2244	*
	2245	*----------------------------------------------------------------------
	2246	*/
	2247
	2248	void
	2249	TclFinalizeDoubleConversion(void)
	2250	{
	2251	int i;
	2252
	2253	Tcl_Free((char *) pow10_wide);
	2254	for (i=0; i<9; ++i) {
	2255	mp_clear(pow5 + i);
	2256	}
	2257	}
	2258
	2259	/*
	2260	*----------------------------------------------------------------------
	2261	*
	2262	* Tcl_InitBignumFromDouble --
	2263	*
	2264	* Extracts the integer part of a double and converts it to an arbitrary
	2265	* precision integer.
	2266	*
	2267	* Results:
	2268	* None.
	2269	*
	2270	* Side effects:
	2271	* Initializes the bignum supplied, and stores the converted number in
	2272	* it.
	2273	*
	2274	*----------------------------------------------------------------------
	2275	*/
	2276
	2277	int
	2278	Tcl_InitBignumFromDouble(
	2279	Tcl_Interp interp, / For error message */
	2280	double d, /* Number to convert */
	2281	mp_int b) / Place to store the result */
	2282	{
	2283	double fract;
	2284	int expt;
	2285
	2286	/*
	2287	* Infinite values can't convert to bignum.
	2288	*/
	2289
	2290	if (TclIsInfinite(d)) {
	2291	if (interp != NULL) {
	2292	const char *s = "integer value too large to represent";
	2293
	2294	Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1));
	2295	Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL);
	2296	}
	2297	return TCL_ERROR;
	2298	}
	2299
	2300	fract = frexp(d,&expt);
	2301	if (expt <= 0) {
	2302	mp_init(b);
	2303	mp_zero(b);
	2304	} else {
	2305	Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
	2306	int shift = expt - mantBits;
	2307
	2308	TclBNInitBignumFromWideInt(b, w);
	2309	if (shift < 0) {
	2310	mp_div_2d(b, -shift, b, NULL);
	2311	} else if (shift > 0) {
	2312	mp_mul_2d(b, shift, b);
	2313	}
	2314	}
	2315	return TCL_OK;
	2316	}
	2317
	2318	/*
	2319	*----------------------------------------------------------------------
	2320	*
	2321	* TclBignumToDouble --
	2322	*
	2323	* Convert an arbitrary-precision integer to a native floating point
	2324	* number.
	2325	*
	2326	* Results:
	2327	* Returns the converted number. Sets errno to ERANGE if the number is
	2328	* too large to convert.
	2329	*
	2330	*----------------------------------------------------------------------
	2331	*/
	2332
	2333	double
	2334	TclBignumToDouble(
	2335	mp_int a) / Integer to convert. */
	2336	{
	2337	mp_int b;
	2338	int bits, shift, i;
	2339	double r;
	2340
	2341	/*
	2342	* Determine how many bits we need, and extract that many from the input.
	2343	* Round to nearest unit in the last place.
	2344	*/
	2345
	2346	bits = mp_count_bits(a);
	2347	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	2348	errno = ERANGE;
	2349	if (a->sign == MP_ZPOS) {
	2350	return HUGE_VAL;
	2351	} else {
	2352	return -HUGE_VAL;
	2353	}
	2354	}
	2355	shift = mantBits + 1 - bits;
	2356	mp_init(&b);
	2357	if (shift > 0) {
	2358	mp_mul_2d(a, shift, &b);
	2359	} else if (shift < 0) {
	2360	mp_div_2d(a, -shift, &b, NULL);
	2361	} else {
	2362	mp_copy(a, &b);
	2363	}
	2364	mp_add_d(&b, 1, &b);
	2365	mp_div_2d(&b, 1, &b, NULL);
	2366
	2367	/*
	2368	* Accumulate the result, one mp_digit at a time.
	2369	*/
	2370
	2371	r = 0.0;
	2372	for (i=b.used-1 ; i>=0 ; --i) {
	2373	r = ldexp(r, DIGIT_BIT) + b.dp[i];
	2374	}
	2375	mp_clear(&b);
	2376
	2377	/*
	2378	* Scale the result to the correct number of bits.
	2379	*/
	2380
	2381	r = ldexp(r, bits - mantBits);
	2382
	2383	/*
	2384	* Return the result with the appropriate sign.
	2385	*/
	2386
	2387	if (a->sign == MP_ZPOS) {
	2388	return r;
	2389	} else {
	2390	return -r;
	2391	}
	2392	}
	2393
	2394	double
	2395	TclCeil(
	2396	mp_int a) / Integer to convert. */
	2397	{
	2398	double r = 0.0;
	2399	mp_int b;
	2400
	2401	mp_init(&b);
	2402	if (mp_cmp_d(a, 0) == MP_LT) {
	2403	mp_neg(a, &b);
	2404	r = -TclFloor(&b);
	2405	} else {
	2406	int bits = mp_count_bits(a);
	2407
	2408	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	2409	r = HUGE_VAL;
	2410	} else {
	2411	int i, exact = 1, shift = mantBits - bits;
	2412
	2413	if (shift > 0) {
	2414	mp_mul_2d(a, shift, &b);
	2415	} else if (shift < 0) {
	2416	mp_int d;
	2417	mp_init(&d);
	2418	mp_div_2d(a, -shift, &b, &d);
	2419	exact = mp_iszero(&d);
	2420	mp_clear(&d);
	2421	} else {
	2422	mp_copy(a, &b);
	2423	}
	2424	if (!exact) {
	2425	mp_add_d(&b, 1, &b);
	2426	}
	2427	for (i=b.used-1 ; i>=0 ; --i) {
	2428	r = ldexp(r, DIGIT_BIT) + b.dp[i];
	2429	}
	2430	r = ldexp(r, bits - mantBits);
	2431	}
	2432	}
	2433	mp_clear(&b);
	2434	return r;
	2435	}
	2436
	2437	double
	2438	TclFloor(
	2439	mp_int a) / Integer to convert. */
	2440	{
	2441	double r = 0.0;
	2442	mp_int b;
	2443
	2444	mp_init(&b);
	2445	if (mp_cmp_d(a, 0) == MP_LT) {
	2446	mp_neg(a, &b);
	2447	r = -TclCeil(&b);
	2448	} else {
	2449	int bits = mp_count_bits(a);
	2450
	2451	if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
	2452	r = DBL_MAX;
	2453	} else {
	2454	int i, shift = mantBits - bits;
	2455
	2456	if (shift > 0) {
	2457	mp_mul_2d(a, shift, &b);
	2458	} else if (shift < 0) {
	2459	mp_div_2d(a, -shift, &b, NULL);
	2460	} else {
	2461	mp_copy(a, &b);
	2462	}
	2463	for (i=b.used-1 ; i>=0 ; --i) {
	2464	r = ldexp(r, DIGIT_BIT) + b.dp[i];
	2465	}
	2466	r = ldexp(r, bits - mantBits);
	2467	}
	2468	}
	2469	mp_clear(&b);
	2470	return r;
	2471	}
	2472
	2473	/*
	2474	*----------------------------------------------------------------------
	2475	*
	2476	* BignumToBiasedFrExp --
	2477	*
	2478	* Convert an arbitrary-precision integer to a native floating point
	2479	* number in the range [0.5,1) times a power of two. NOTE: Intentionally
	2480	* converts to a number that's a few ulp too small, so that
	2481	* RefineApproximation will not overflow near the high end of the
	2482	* machine's arithmetic range.
	2483	*
	2484	* Results:
	2485	* Returns the converted number.
	2486	*
	2487	* Side effects:
	2488	* Stores the exponent of two in 'machexp'.
	2489	*
	2490	*----------------------------------------------------------------------
	2491	*/
	2492
	2493	static double
	2494	BignumToBiasedFrExp(
	2495	mp_int a, / Integer to convert */
	2496	int machexp) / Power of two */
	2497	{
	2498	mp_int b;
	2499	int bits;
	2500	int shift;
	2501	int i;
	2502	double r;
	2503
	2504	/*
	2505	* Determine how many bits we need, and extract that many from the input.
	2506	* Round to nearest unit in the last place.
	2507	*/
	2508
	2509	bits = mp_count_bits(a);
	2510	shift = mantBits - 2 - bits;
	2511	mp_init(&b);
	2512	if (shift > 0) {
	2513	mp_mul_2d(a, shift, &b);
	2514	} else if (shift < 0) {
	2515	mp_div_2d(a, -shift, &b, NULL);
	2516	} else {
	2517	mp_copy(a, &b);
	2518	}
	2519
	2520	/*
	2521	* Accumulate the result, one mp_digit at a time.
	2522	*/
	2523
	2524	r = 0.0;
	2525	for (i=b.used-1; i>=0; --i) {
	2526	r = ldexp(r, DIGIT_BIT) + b.dp[i];
	2527	}
	2528	mp_clear(&b);
	2529
	2530	/*
	2531	* Return the result with the appropriate sign.
	2532	*/
	2533
	2534	*machexp = bits - mantBits + 2;
	2535	return ((a->sign == MP_ZPOS) ? r : -r);
	2536	}
	2537
	2538	/*
	2539	*----------------------------------------------------------------------
	2540	*
	2541	* Pow10TimesFrExp --
	2542	*
	2543	* Multiply a power of ten by a number expressed as fraction and
	2544	* exponent.
	2545	*
	2546	* Results:
	2547	* Returns the significand of the result.
	2548	*
	2549	* Side effects:
	2550	* Overwrites the 'machexp' parameter with the exponent of the result.
	2551	*
	2552	* Assumes that 'exponent' is such that 10**exponent would be a double, even
	2553	* though 'fraction10*(machexp+exponent)' might overflow.
	2554	*
	2555	*----------------------------------------------------------------------
	2556	*/
	2557
	2558	static double
	2559	Pow10TimesFrExp(
	2560	int exponent, /* Power of 10 to multiply by */
	2561	double fraction, /* Significand of multiplicand */
	2562	int machexp) / On input, exponent of multiplicand. On
	2563	* output, exponent of result. */
	2564	{
	2565	int i, j;
	2566	int expt = *machexp;
	2567	double retval = fraction;
	2568
	2569	if (exponent > 0) {
	2570	/*
	2571	* Multiply by 10**exponent
	2572	*/
	2573
	2574	retval = frexp(retval * pow10vals[exponent&0xf], &j);
	2575	expt += j;
	2576	for (i=4; i<9; ++i) {
	2577	if (exponent & (1<<i)) {
	2578	retval = frexp(retval * pow_10_2_n[i], &j);
	2579	expt += j;
	2580	}
	2581	}
	2582	} else if (exponent < 0) {
	2583	/*
	2584	* Divide by 10**-exponent
	2585	*/
	2586
	2587	retval = frexp(retval / pow10vals[(-exponent) & 0xf], &j);
	2588	expt += j;
	2589	for (i=4; i<9; ++i) {
	2590	if ((-exponent) & (1<<i)) {
	2591	retval = frexp(retval / pow_10_2_n[i], &j);
	2592	expt += j;
	2593	}
	2594	}
	2595	}
	2596
	2597	*machexp = expt;
	2598	return retval;
	2599	}
	2600
	2601	/*
	2602	*----------------------------------------------------------------------
	2603	*
	2604	* SafeLdExp --
	2605	*
	2606	* Do an 'ldexp' operation, but handle denormals gracefully.
	2607	*
	2608	* Results:
	2609	* Returns the appropriately scaled value.
	2610	*
	2611	* On some platforms, 'ldexp' fails when presented with a number too
	2612	* small to represent as a normalized double. This routine does 'ldexp'
	2613	* in two steps for those numbers, to return correctly denormalized
	2614	* values.
	2615	*
	2616	*----------------------------------------------------------------------
	2617	*/
	2618
	2619	static double
	2620	SafeLdExp(
	2621	double fract,
	2622	int expt)
	2623	{
	2624	int minexpt = DBL_MIN_EXP * log2FLT_RADIX;
	2625	volatile double a, b, retval;
	2626
	2627	if (expt < minexpt) {
	2628	a = ldexp(fract, expt - mantBits - minexpt);
	2629	b = ldexp(1.0, mantBits + minexpt);
	2630	retval = a * b;
	2631	} else {
	2632	retval = ldexp(fract, expt);
	2633	}
	2634	return retval;
	2635	}
	2636
	2637	/*
	2638	*----------------------------------------------------------------------
	2639	*
	2640	* TclFormatNaN --
	2641	*
	2642	* Makes the string representation of a "Not a Number"
	2643	*
	2644	* Results:
	2645	* None.
	2646	*
	2647	* Side effects:
	2648	* Stores the string representation in the supplied buffer, which must be
	2649	* at least TCL_DOUBLE_SPACE characters.
	2650	*
	2651	*----------------------------------------------------------------------
	2652	*/
	2653
	2654	void
	2655	TclFormatNaN(
	2656	double value, /* The Not-a-Number to format. */
	2657	char buffer) / String representation. */
	2658	{
	2659	#ifndef IEEE_FLOATING_POINT
	2660	strcpy(buffer, "NaN");
	2661	return;
	2662	#else
	2663	union {
	2664	double dv;
	2665	Tcl_WideUInt iv;
	2666	} bitwhack;
	2667
	2668	bitwhack.dv = value;
	2669	if (n770_fp) {
	2670	bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
	2671	}
	2672	if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
	2673	bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63);
	2674	*buffer++ = '-';
	2675	}
	2676	*buffer++ = 'N';
	2677	*buffer++ = 'a';
	2678	*buffer++ = 'N';
	2679	bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
	2680	if (bitwhack.iv != 0) {
	2681	sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv);
	2682	} else {
	2683	*buffer = '\0';
	2684	}
	2685	#endif /* IEEE_FLOATING_POINT */
	2686	}
	2687
	2688	/*
	2689	*----------------------------------------------------------------------
	2690	*
	2691	* Nokia770Twiddle --
	2692	*
	2693	* Transpose the two words of a number for Nokia 770 floating
	2694	* point handling.
	2695	*
	2696	*----------------------------------------------------------------------
	2697	*/
	2698
	2699	static Tcl_WideUInt
	2700	Nokia770Twiddle(
	2701	Tcl_WideUInt w) /* Number to transpose */
	2702	{
	2703	return (((w >> 32) & 0xffffffff) \| (w << 32));
	2704	}
	2705
	2706	/*
	2707	*----------------------------------------------------------------------
	2708	*
	2709	* TclNokia770Doubles --
	2710	*
	2711	* Transpose the two words of a number for Nokia 770 floating
	2712	* point handling.
	2713	*
	2714	*----------------------------------------------------------------------
	2715	*/
	2716
	2717	int
	2718	TclNokia770Doubles(void)
	2719	{
	2720	return n770_fp;
	2721	}
	2722
	2723	/*
	2724	* Local Variables:
	2725	* mode: c
	2726	* c-basic-offset: 4
	2727	* fill-column: 78
	2728	* End:
	2729	*/

Note: See TracBrowser for help on using the repository browser.

Download in other formats: