Planet

navi

home

PPS

about

screenshots

download

development

forum

Context Navigation

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

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

Line
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: