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 */
18
19 #include "tclInt.h"
20 #include "tommath.h"
21 #include <math.h>
22
23 /*
24 * Define KILL_OCTAL to suppress interpretation of numbers with leading zero
25 * as octal. (Ceterum censeo: numeros octonarios delendos esse.)
26 */
27
28 #undef KILL_OCTAL
29
30 /*
31 * This code supports (at least hypothetically), IBM, Cray, VAX and IEEE-754
32 * floating point; of these, only IEEE-754 can represent NaN. IEEE-754 can be
33 * uniquely determined by radix and by the widths of significand and exponent.
34 */
35
36 #if (FLT_RADIX == 2) && (DBL_MANT_DIG == 53) && (DBL_MAX_EXP == 1024)
37 # define IEEE_FLOATING_POINT
38 #endif
39
40 /*
41 * gcc on x86 needs access to rounding controls, because of a questionable
42 * feature where it retains intermediate results as IEEE 'long double' values
43 * somewhat unpredictably. It is tempting to include fpu_control.h, but that
44 * file exists only on Linux; it is missing on Cygwin and MinGW. Most gcc-isms
45 * and ix86-isms are factored out here.
46 */
47
48 #if defined(__GNUC__) && defined(__i386)
49 typedef unsigned int fpu_control_t __attribute__ ((__mode__ (__HI__)));
50 #define _FPU_GETCW(cw) __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw))
51 #define _FPU_SETCW(cw) __asm__ __volatile__ ("fldcw %0" : : "m" (*&cw))
52 # define FPU_IEEE_ROUNDING 0x027f
53 # define ADJUST_FPU_CONTROL_WORD
54 #endif
55
56 /* Sun ProC needs sunmath for rounding control on x86 like gcc above.
57 *
58 *
59 */
60 #if defined(__sun) && defined(__i386) && !defined(__GNUC__)
61 #include <sunmath.h>
62 #endif
63
64 /*
65 * MIPS floating-point units need special settings in control registers
66 * to use gradual underflow as we expect. This fix is for the MIPSpro
67 * compiler.
68 */
69 #if defined(__sgi) && defined(_COMPILER_VERSION)
70 #include <sys/fpu.h>
71 #endif
72 /*
73 * HP's PA_RISC architecture uses 7ff4000000000000 to represent a quiet NaN.
74 * Everyone else uses 7ff8000000000000. (Why, HP, why?)
75 */
76
77 #ifdef __hppa
78 # define NAN_START 0x7ff4
79 # define NAN_MASK (((Tcl_WideUInt) 1) << 50)
80 #else
81 # define NAN_START 0x7ff8
82 # define NAN_MASK (((Tcl_WideUInt) 1) << 51)
83 #endif
84
85 /*
86 * Constants used by this file (most of which are only ever calculated at
87 * runtime).
88 */
89
90 /* Magic constants */
91
92 #define LOG10_2 0.3010299956639812
93 #define TWO_OVER_3LOG10 0.28952965460216784
94 #define LOG10_3HALVES_PLUS_FUDGE 0.1760912590558
95
96 /* Definitions of the parts of an IEEE754-format floating point number */
97
98 #define SIGN_BIT 0x80000000
99 /* Mask for the sign bit in the first
100 * word of a double */
101 #define EXP_MASK 0x7ff00000
102 /* Mask for the exponent field in the
103 * first word of a double */
104 #define EXP_SHIFT 20
105 /* Shift count to make the exponent an
106 * integer */
107 #define HIDDEN_BIT (((Tcl_WideUInt) 0x00100000) << 32)
108 /* Hidden 1 bit for the significand */
109 #define HI_ORDER_SIG_MASK 0x000fffff
110 /* Mask for the high-order part of the
111 * significand in the first word of a
112 * double */
113 #define SIG_MASK (((Tcl_WideUInt) HI_ORDER_SIG_MASK << 32) \
114 | 0xffffffff)
115 /* Mask for the 52-bit significand. */
116 #define FP_PRECISION 53
117 /* Number of bits of significand plus the
118 * hidden bit */
119 #define EXPONENT_BIAS 0x3ff
120 /* Bias of the exponent 0 */
121
122 /* Derived quantities */
123
124 #define TEN_PMAX 22
125 /* floor(FP_PRECISION*log(2)/log(5)) */
126 #define QUICK_MAX 14
127 /* floor((FP_PRECISION-1)*log(2)/log(10)) - 1 */
128 #define BLETCH 0x10
129 /* Highest power of two that is greater than
130 * DBL_MAX_10_EXP, divided by 16 */
131 #define DIGIT_GROUP 8
132 /* floor(DIGIT_BIT*log(2)/log(10)) */
133
134 /* Union used to dismantle floating point numbers. */
135
136 typedef union Double {
137 struct {
138 #ifdef WORDS_BIGENDIAN
139 int word0;
140 int word1;
141 #else
142 int word1;
143 int word0;
144 #endif
145 } w;
146 double d;
147 Tcl_WideUInt q;
148 } Double;
149
150 static int maxpow10_wide; /* The powers of ten that can be represented
151 * exactly as wide integers. */
152 static Tcl_WideUInt *pow10_wide;
153 #define MAXPOW 22
154 static double pow10vals[MAXPOW+1];
155 /* The powers of ten that can be represented
156 * exactly as IEEE754 doubles. */
157 static int mmaxpow; /* Largest power of ten that can be
158 * represented exactly in a 'double'. */
159 static int log10_DIGIT_MAX; /* The number of decimal digits that fit in an
160 * mp_digit. */
161 static int log2FLT_RADIX; /* Logarithm of the floating point radix. */
162 static int mantBits; /* Number of bits in a double's significand */
163 static mp_int pow5[9]; /* Table of powers of 5**(2**n), up to
164 * 5**256 */
165 static double tiny = 0.0; /* The smallest representable double */
166 static int maxDigits; /* The maximum number of digits to the left of
167 * the decimal point of a double. */
168 static int minDigits; /* The maximum number of digits to the right
169 * of the decimal point in a double. */
170 static const double pow_10_2_n[] = { /* Inexact higher powers of ten. */
171 1.0,
172 100.0,
173 10000.0,
174 1.0e+8,
175 1.0e+16,
176 1.0e+32,
177 1.0e+64,
178 1.0e+128,
179 1.0e+256
180 };
181
182 static int n770_fp; /* Flag is 1 on Nokia N770 floating point.
183 * Nokia's floating point has the words
184 * reversed: if big-endian is 7654 3210,
185 * and little-endian is 0123 4567,
186 * then Nokia's FP is 4567 0123;
187 * little-endian within the 32-bit words
188 * but big-endian between them. */
189
190 /* Table of powers of 5 that are small enough to fit in an mp_digit. */
191
192 static const mp_digit dpow5[13] = {
193 1, 5, 25, 125,
194 625, 3125, 15625, 78125,
195 390625, 1953125, 9765625, 48828125,
196 244140625
197 };
198
199 /* Table of powers: pow5_13[n] = 5**(13*2**(n+1)) */
200 static mp_int pow5_13[5]; /* Table of powers: 5**13, 5**26, 5**52,
201 * 5**104, 5**208 */
202 static const double tens[] = {
203 1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08, 1e09,
204 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
205 1e20, 1e21, 1e22
206 };
207
208 static const int itens [] = {
209 1,
210 10,
211 100,
212 1000,
213 10000,
214 100000,
215 1000000,
216 10000000,
217 100000000
218 };
219
220 static const double bigtens[] = {
221 1e016, 1e032, 1e064, 1e128, 1e256
222 };
223 #define N_BIGTENS 5
224
225 static const int log2pow5[27] = {
226 01, 3, 5, 7, 10, 12, 14, 17, 19, 21,
227 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
228 47, 49, 52, 54, 56, 59, 61
229 };
230 #define N_LOG2POW5 27
231
232 static const Tcl_WideUInt wuipow5[27] = {
233 (Tcl_WideUInt) 1, /* 5**0 */
234 (Tcl_WideUInt) 5,
235 (Tcl_WideUInt) 25,
236 (Tcl_WideUInt) 125,
237 (Tcl_WideUInt) 625,
238 (Tcl_WideUInt) 3125, /* 5**5 */
239 (Tcl_WideUInt) 3125*5,
240 (Tcl_WideUInt) 3125*25,
241 (Tcl_WideUInt) 3125*125,
242 (Tcl_WideUInt) 3125*625,
243 (Tcl_WideUInt) 3125*3125, /* 5**10 */
244 (Tcl_WideUInt) 3125*3125*5,
245 (Tcl_WideUInt) 3125*3125*25,
246 (Tcl_WideUInt) 3125*3125*125,
247 (Tcl_WideUInt) 3125*3125*625,
248 (Tcl_WideUInt) 3125*3125*3125, /* 5**15 */
249 (Tcl_WideUInt) 3125*3125*3125*5,
250 (Tcl_WideUInt) 3125*3125*3125*25,
251 (Tcl_WideUInt) 3125*3125*3125*125,
252 (Tcl_WideUInt) 3125*3125*3125*625,
253 (Tcl_WideUInt) 3125*3125*3125*3125, /* 5**20 */
254 (Tcl_WideUInt) 3125*3125*3125*3125*5,
255 (Tcl_WideUInt) 3125*3125*3125*3125*25,
256 (Tcl_WideUInt) 3125*3125*3125*3125*125,
257 (Tcl_WideUInt) 3125*3125*3125*3125*625,
258 (Tcl_WideUInt) 3125*3125*3125*3125*3125, /* 5**25 */
259 (Tcl_WideUInt) 3125*3125*3125*3125*3125*5 /* 5**26 */
260 };
261
262 /*
263 * Static functions defined in this file.
264 */
265
266 static int AccumulateDecimalDigit(unsigned, int,
267 Tcl_WideUInt *, mp_int *, int);
268 static double MakeHighPrecisionDouble(int signum,
269 mp_int *significand, int nSigDigs, int exponent);
270 static double MakeLowPrecisionDouble(int signum,
271 Tcl_WideUInt significand, int nSigDigs,
272 int exponent);
273 static double MakeNaN(int signum, Tcl_WideUInt tag);
274 static double RefineApproximation(double approx,
275 mp_int *exactSignificand, int exponent);
276 static void MulPow5(mp_int*, unsigned, mp_int*);
277 static int NormalizeRightward(Tcl_WideUInt*);
278 static int RequiredPrecision(Tcl_WideUInt);
279 static void DoubleToExpAndSig(double, Tcl_WideUInt*, int*, int*);
280 static void TakeAbsoluteValue(Double*, int*);
281 static char* FormatInfAndNaN(Double*, int*, char**);
282 static char* FormatZero(int*, char**);
283 static int ApproximateLog10(Tcl_WideUInt, int, int);
284 static int BetterLog10(double, int, int*);
285 static void ComputeScale(int, int, int*, int*, int*, int*);
286 static void SetPrecisionLimits(int, int, int*, int*, int*, int*);
287 static char* BumpUp(char*, char*, int*);
288 static int AdjustRange(double*, int);
289 static char* ShorteningQuickFormat(double, int, int, double,
290 char*, int*);
291 static char* StrictQuickFormat(double, int, int, double,
292 char*, int*);
293 static char* QuickConversion(double, int, int, int, int, int, int,
294 int*, char**);
295 static void CastOutPowersOf2(int*, int*, int*);
296 static char* ShorteningInt64Conversion(Double*, int, Tcl_WideUInt,
297 int, int, int, int, int, int, int, int, int,
298 int, int, int*, char**);
299 static char* StrictInt64Conversion(Double*, int, Tcl_WideUInt,
300 int, int, int, int, int, int,
301 int, int, int*, char**);
302 static int ShouldBankerRoundUpPowD(mp_int*, int, int);
303 static int ShouldBankerRoundUpToNextPowD(mp_int*, mp_int*,
304 int, int, int, mp_int*);
305 static char* ShorteningBignumConversionPowD(Double* dPtr,
306 int convType, Tcl_WideUInt bw, int b2, int b5,
307 int m2plus, int m2minus, int m5,
308 int sd, int k, int len,
309 int ilim, int ilim1, int* decpt,
310 char** endPtr);
311 static char* StrictBignumConversionPowD(Double* dPtr, int convType,
312 Tcl_WideUInt bw, int b2, int b5,
313 int sd, int k, int len,
314 int ilim, int ilim1, int* decpt,
315 char** endPtr);
316 static int ShouldBankerRoundUp(mp_int*, mp_int*, int);
317 static int ShouldBankerRoundUpToNext(mp_int*, mp_int*, mp_int*,
318 int, int, mp_int*);
319 static char* ShorteningBignumConversion(Double* dPtr, int convType,
320 Tcl_WideUInt bw, int b2,
321 int m2plus, int m2minus,
322 int s2, int s5, int k, int len,
323 int ilim, int ilim1, int* decpt,
324 char** endPtr);
325 static char* StrictBignumConversion(Double* dPtr, int convType,
326 Tcl_WideUInt bw, int b2,
327 int s2, int s5, int k, int len,
328 int ilim, int ilim1, int* decpt,
329 char** endPtr);
330 static double BignumToBiasedFrExp(mp_int *big, int *machexp);
331 static double Pow10TimesFrExp(int exponent, double fraction,
332 int *machexp);
333 static double SafeLdExp(double fraction, int exponent);
334 static Tcl_WideUInt Nokia770Twiddle(Tcl_WideUInt w);
335
336 /*
337 *----------------------------------------------------------------------
338 *
339 * TclParseNumber --
340 *
341 * Scans bytes, interpreted as characters in Tcl's internal encoding, and
342 * parses the longest prefix that is the string representation of a
343 * number in a format recognized by Tcl.
344 *
345 * The arguments bytes, numBytes, and objPtr are the inputs which
346 * determine the string to be parsed. If bytes is non-NULL, it points to
347 * the first byte to be scanned. If bytes is NULL, then objPtr must be
348 * non-NULL, and the string representation of objPtr will be scanned
349 * (generated first, if necessary). The numBytes argument determines the
350 * number of bytes to be scanned. If numBytes is negative, the first NUL
351 * byte encountered will terminate the scan. If numBytes is non-negative,
352 * then no more than numBytes bytes will be scanned.
353 *
354 * The argument flags is an input that controls the numeric formats
355 * recognized by the parser. The flag bits are:
356 *
357 * - TCL_PARSE_INTEGER_ONLY: accept only integer values; reject
358 * strings that denote floating point values (or accept only the
359 * leading portion of them that are integer values).
360 * - TCL_PARSE_SCAN_PREFIXES: ignore the prefixes 0b and 0o that are
361 * not part of the [scan] command's vocabulary. Use only in
362 * combination with TCL_PARSE_INTEGER_ONLY.
363 * - TCL_PARSE_OCTAL_ONLY: parse only in the octal format, whether
364 * or not a prefix is present that would lead to octal parsing.
365 * Use only in combination with TCL_PARSE_INTEGER_ONLY.
366 * - TCL_PARSE_HEXADECIMAL_ONLY: parse only in the hexadecimal format,
367 * whether or not a prefix is present that would lead to
368 * hexadecimal parsing. Use only in combination with
369 * TCL_PARSE_INTEGER_ONLY.
370 * - TCL_PARSE_DECIMAL_ONLY: parse only in the decimal format, no
371 * matter whether a 0 prefix would normally force a different
372 * base.
373 * - TCL_PARSE_NO_WHITESPACE: reject any leading/trailing whitespace
374 *
375 * The arguments interp and expected are inputs that control error
376 * message generation. If interp is NULL, no error message will be
377 * generated. If interp is non-NULL, then expected must also be non-NULL.
378 * When TCL_ERROR is returned, an error message will be left in the
379 * result of interp, and the expected argument will appear in the error
380 * message as the thing TclParseNumber expected, but failed to find in
381 * the string.
382 *
383 * The arguments objPtr and endPtrPtr as well as the return code are the
384 * outputs.
385 *
386 * When the parser cannot find any prefix of the string that matches a
387 * format it is looking for, TCL_ERROR is returned and an error message
388 * may be generated and returned as described above. The contents of
389 * objPtr will not be changed. If endPtrPtr is non-NULL, a pointer to the
390 * character in the string that terminated the scan will be written to
391 * *endPtrPtr.
392 *
393 * When the parser determines that the entire string matches a format it
394 * is looking for, TCL_OK is returned, and if objPtr is non-NULL, then
395 * the internal rep and Tcl_ObjType of objPtr are set to the "canonical"
396 * numeric value that matches the scanned string. If endPtrPtr is not
397 * NULL, a pointer to the end of the string will be written to *endPtrPtr
398 * (that is, either bytes+numBytes or a pointer to a terminating NUL
399 * byte).
400 *
401 * When the parser determines that a partial string matches a format it
402 * is looking for, the value of endPtrPtr determines what happens:
403 *
404 * - If endPtrPtr is NULL, then TCL_ERROR is returned, with error message
405 * generation as above.
406 *
407 * - If endPtrPtr is non-NULL, then TCL_OK is returned and objPtr
408 * internals are set as above. Also, a pointer to the first
409 * character following the parsed numeric string is written to
410 * *endPtrPtr.
411 *
412 * In some cases where the string being scanned is the string rep of
413 * objPtr, this routine can leave objPtr in an inconsistent state where
414 * its string rep and its internal rep do not agree. In these cases the
415 * internal rep will be in agreement with only some substring of the
416 * string rep. This might happen if the caller passes in a non-NULL bytes
417 * value that points somewhere into the string rep. It might happen if
418 * the caller passes in a numBytes value that limits the scan to only a
419 * prefix of the string rep. Or it might happen if a non-NULL value of
420 * endPtrPtr permits a TCL_OK return from only a partial string match. It
421 * is the responsibility of the caller to detect and correct such
422 * inconsistencies when they can and do arise.
423 *
424 * Results:
425 * Returns a standard Tcl result.
426 *
427 * Side effects:
428 * The string representaton of objPtr may be generated.
429 *
430 * The internal representation and Tcl_ObjType of objPtr may be changed.
431 * This may involve allocation and/or freeing of memory.
432 *
433 *----------------------------------------------------------------------
434 */
435
436 int
TclParseNumber(Tcl_Interp * interp,Tcl_Obj * objPtr,const char * expected,const char * bytes,int numBytes,const char ** endPtrPtr,int flags)437 TclParseNumber(
438 Tcl_Interp *interp, /* Used for error reporting. May be NULL. */
439 Tcl_Obj *objPtr, /* Object to receive the internal rep. */
440 const char *expected, /* Description of the type of number the
441 * caller expects to be able to parse
442 * ("integer", "boolean value", etc.). */
443 const char *bytes, /* Pointer to the start of the string to
444 * scan. */
445 int numBytes, /* Maximum number of bytes to scan, see
446 * above. */
447 const char **endPtrPtr, /* Place to store pointer to the character
448 * that terminated the scan. */
449 int flags) /* Flags governing the parse. */
450 {
451 enum State {
452 INITIAL, SIGNUM, ZERO, ZERO_X,
453 ZERO_O, ZERO_B, BINARY,
454 HEXADECIMAL, OCTAL, BAD_OCTAL, DECIMAL,
455 LEADING_RADIX_POINT, FRACTION,
456 EXPONENT_START, EXPONENT_SIGNUM, EXPONENT,
457 sI, sIN, sINF, sINFI, sINFIN, sINFINI, sINFINIT, sINFINITY
458 #ifdef IEEE_FLOATING_POINT
459 , sN, sNA, sNAN, sNANPAREN, sNANHEX, sNANFINISH
460 #endif
461 } state = INITIAL;
462 enum State acceptState = INITIAL;
463
464 int signum = 0; /* Sign of the number being parsed */
465 Tcl_WideUInt significandWide = 0;
466 /* Significand of the number being parsed (if
467 * no overflow) */
468 mp_int significandBig; /* Significand of the number being parsed (if
469 * it overflows significandWide) */
470 int significandOverflow = 0;/* Flag==1 iff significandBig is used */
471 Tcl_WideUInt octalSignificandWide = 0;
472 /* Significand of an octal number; needed
473 * because we don't know whether a number with
474 * a leading zero is octal or decimal until
475 * we've scanned forward to a '.' or 'e' */
476 mp_int octalSignificandBig; /* Significand of octal number once
477 * octalSignificandWide overflows */
478 int octalSignificandOverflow = 0;
479 /* Flag==1 if octalSignificandBig is used */
480 int numSigDigs = 0; /* Number of significant digits in the decimal
481 * significand */
482 int numTrailZeros = 0; /* Number of trailing zeroes at the current
483 * point in the parse. */
484 int numDigitsAfterDp = 0; /* Number of digits scanned after the decimal
485 * point */
486 int exponentSignum = 0; /* Signum of the exponent of a floating point
487 * number */
488 long exponent = 0; /* Exponent of a floating point number */
489 const char *p; /* Pointer to next character to scan */
490 size_t len; /* Number of characters remaining after p */
491 const char *acceptPoint; /* Pointer to position after last character in
492 * an acceptable number */
493 size_t acceptLen; /* Number of characters following that
494 * point. */
495 int status = TCL_OK; /* Status to return to caller */
496 char d = 0; /* Last hexadecimal digit scanned; initialized
497 * to avoid a compiler warning. */
498 int shift = 0; /* Amount to shift when accumulating binary */
499 int explicitOctal = 0;
500
501 #define ALL_BITS (~(Tcl_WideUInt)0)
502 #define MOST_BITS (ALL_BITS >> 1)
503
504 /*
505 * Initialize bytes to start of the object's string rep if the caller
506 * didn't pass anything else.
507 */
508
509 if (bytes == NULL) {
510 bytes = TclGetString(objPtr);
511 }
512
513 p = bytes;
514 len = numBytes;
515 acceptPoint = p;
516 acceptLen = len;
517 while (1) {
518 char c = len ? *p : '\0';
519 switch (state) {
520
521 case INITIAL:
522 /*
523 * Initial state. Acceptable characters are +, -, digits, period,
524 * I, N, and whitespace.
525 */
526
527 if (TclIsSpaceProc(c)) {
528 if (flags & TCL_PARSE_NO_WHITESPACE) {
529 goto endgame;
530 }
531 break;
532 } else if (c == '+') {
533 state = SIGNUM;
534 break;
535 } else if (c == '-') {
536 signum = 1;
537 state = SIGNUM;
538 break;
539 }
540 /* FALLTHROUGH */
541
542 case SIGNUM:
543 /*
544 * Scanned a leading + or -. Acceptable characters are digits,
545 * period, I, and N.
546 */
547
548 if (c == '0') {
549 if (flags & TCL_PARSE_DECIMAL_ONLY) {
550 state = DECIMAL;
551 } else {
552 state = ZERO;
553 }
554 break;
555 } else if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
556 goto zerox;
557 } else if (flags & TCL_PARSE_OCTAL_ONLY) {
558 goto zeroo;
559 } else if (isdigit(UCHAR(c))) {
560 significandWide = c - '0';
561 numSigDigs = 1;
562 state = DECIMAL;
563 break;
564 } else if (flags & TCL_PARSE_INTEGER_ONLY) {
565 goto endgame;
566 } else if (c == '.') {
567 state = LEADING_RADIX_POINT;
568 break;
569 } else if (c == 'I' || c == 'i') {
570 state = sI;
571 break;
572 #ifdef IEEE_FLOATING_POINT
573 } else if (c == 'N' || c == 'n') {
574 state = sN;
575 break;
576 #endif
577 }
578 goto endgame;
579
580 case ZERO:
581 /*
582 * Scanned a leading zero (perhaps with a + or -). Acceptable
583 * inputs are digits, period, X, b, and E. If 8 or 9 is encountered,
584 * the number can't be octal. This state and the OCTAL state
585 * differ only in whether they recognize 'X' and 'b'.
586 */
587
588 acceptState = state;
589 acceptPoint = p;
590 acceptLen = len;
591 if (c == 'x' || c == 'X') {
592 state = ZERO_X;
593 break;
594 }
595 if (flags & TCL_PARSE_HEXADECIMAL_ONLY) {
596 goto zerox;
597 }
598 if (flags & TCL_PARSE_SCAN_PREFIXES) {
599 goto zeroo;
600 }
601 if (c == 'b' || c == 'B') {
602 state = ZERO_B;
603 break;
604 }
605 if (c == 'o' || c == 'O') {
606 explicitOctal = 1;
607 state = ZERO_O;
608 break;
609 }
610 #ifdef KILL_OCTAL
611 goto decimal;
612 #endif
613 /* FALLTHROUGH */
614
615 case OCTAL:
616 /*
617 * Scanned an optional + or -, followed by a string of octal
618 * digits. Acceptable inputs are more digits, period, or E. If 8
619 * or 9 is encountered, commit to floating point.
620 */
621
622 acceptState = state;
623 acceptPoint = p;
624 acceptLen = len;
625 /* FALLTHROUGH */
626 case ZERO_O:
627 zeroo:
628 if (c == '0') {
629 numTrailZeros++;
630 state = OCTAL;
631 break;
632 } else if (c >= '1' && c <= '7') {
633 if (objPtr != NULL) {
634 shift = 3 * (numTrailZeros + 1);
635 significandOverflow = AccumulateDecimalDigit(
636 (unsigned)(c-'0'), numTrailZeros,
637 &significandWide, &significandBig,
638 significandOverflow);
639
640 if (!octalSignificandOverflow) {
641 /*
642 * Shifting by more bits than are in the value being
643 * shifted is at least de facto nonportable. Check for
644 * too large shifts first.
645 */
646
647 if ((octalSignificandWide != 0)
648 && (((size_t)shift >=
649 CHAR_BIT*sizeof(Tcl_WideUInt))
650 || (octalSignificandWide >
651 (~(Tcl_WideUInt)0 >> shift)))) {
652 octalSignificandOverflow = 1;
653 TclBNInitBignumFromWideUInt(&octalSignificandBig,
654 octalSignificandWide);
655 }
656 }
657 if (!octalSignificandOverflow) {
658 octalSignificandWide =
659 (octalSignificandWide << shift) + (c - '0');
660 } else {
661 mp_mul_2d(&octalSignificandBig, shift,
662 &octalSignificandBig);
663 mp_add_d(&octalSignificandBig, (mp_digit)(c - '0'),
664 &octalSignificandBig);
665 }
666 }
667 if (numSigDigs != 0) {
668 numSigDigs += numTrailZeros+1;
669 } else {
670 numSigDigs = 1;
671 }
672 numTrailZeros = 0;
673 state = OCTAL;
674 break;
675 }
676 /* FALLTHROUGH */
677
678 case BAD_OCTAL:
679 if (explicitOctal) {
680 /*
681 * No forgiveness for bad digits in explicitly octal numbers.
682 */
683
684 goto endgame;
685 }
686 if (flags & TCL_PARSE_INTEGER_ONLY) {
687 /*
688 * No seeking floating point when parsing only integer.
689 */
690
691 goto endgame;
692 }
693 #ifndef KILL_OCTAL
694
695 /*
696 * Scanned a number with a leading zero that contains an 8, 9,
697 * radix point or E. This is an invalid octal number, but might
698 * still be floating point.
699 */
700
701 if (c == '0') {
702 numTrailZeros++;
703 state = BAD_OCTAL;
704 break;
705 } else if (isdigit(UCHAR(c))) {
706 if (objPtr != NULL) {
707 significandOverflow = AccumulateDecimalDigit(
708 (unsigned)(c-'0'), numTrailZeros,
709 &significandWide, &significandBig,
710 significandOverflow);
711 }
712 if (numSigDigs != 0) {
713 numSigDigs += (numTrailZeros + 1);
714 } else {
715 numSigDigs = 1;
716 }
717 numTrailZeros = 0;
718 state = BAD_OCTAL;
719 break;
720 } else if (c == '.') {
721 state = FRACTION;
722 break;
723 } else if (c == 'E' || c == 'e') {
724 state = EXPONENT_START;
725 break;
726 }
727 #endif
728 goto endgame;
729
730 /*
731 * Scanned 0x. If state is HEXADECIMAL, scanned at least one
732 * character following the 0x. The only acceptable inputs are
733 * hexadecimal digits.
734 */
735
736 case HEXADECIMAL:
737 acceptState = state;
738 acceptPoint = p;
739 acceptLen = len;
740 /* FALLTHROUGH */
741
742 case ZERO_X:
743 zerox:
744 if (c == '0') {
745 numTrailZeros++;
746 state = HEXADECIMAL;
747 break;
748 } else if (isdigit(UCHAR(c))) {
749 d = (c-'0');
750 } else if (c >= 'A' && c <= 'F') {
751 d = (c-'A'+10);
752 } else if (c >= 'a' && c <= 'f') {
753 d = (c-'a'+10);
754 } else {
755 goto endgame;
756 }
757 if (objPtr != NULL) {
758 shift = 4 * (numTrailZeros + 1);
759 if (!significandOverflow) {
760 /*
761 * Shifting by more bits than are in the value being
762 * shifted is at least de facto nonportable. Check for too
763 * large shifts first.
764 */
765
766 if (significandWide != 0 &&
767 ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
768 significandWide > (~(Tcl_WideUInt)0 >> shift))) {
769 significandOverflow = 1;
770 TclBNInitBignumFromWideUInt(&significandBig,
771 significandWide);
772 }
773 }
774 if (!significandOverflow) {
775 significandWide = (significandWide << shift) + d;
776 } else {
777 mp_mul_2d(&significandBig, shift, &significandBig);
778 mp_add_d(&significandBig, (mp_digit) d, &significandBig);
779 }
780 }
781 numTrailZeros = 0;
782 state = HEXADECIMAL;
783 break;
784
785 case BINARY:
786 acceptState = state;
787 acceptPoint = p;
788 acceptLen = len;
789 case ZERO_B:
790 if (c == '0') {
791 numTrailZeros++;
792 state = BINARY;
793 break;
794 } else if (c != '1') {
795 goto endgame;
796 }
797 if (objPtr != NULL) {
798 shift = numTrailZeros + 1;
799 if (!significandOverflow) {
800 /*
801 * Shifting by more bits than are in the value being
802 * shifted is at least de facto nonportable. Check for too
803 * large shifts first.
804 */
805
806 if (significandWide != 0 &&
807 ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
808 significandWide > (~(Tcl_WideUInt)0 >> shift))) {
809 significandOverflow = 1;
810 TclBNInitBignumFromWideUInt(&significandBig,
811 significandWide);
812 }
813 }
814 if (!significandOverflow) {
815 significandWide = (significandWide << shift) + 1;
816 } else {
817 mp_mul_2d(&significandBig, shift, &significandBig);
818 mp_add_d(&significandBig, (mp_digit) 1, &significandBig);
819 }
820 }
821 numTrailZeros = 0;
822 state = BINARY;
823 break;
824
825 case DECIMAL:
826 /*
827 * Scanned an optional + or - followed by a string of decimal
828 * digits.
829 */
830
831 #ifdef KILL_OCTAL
832 decimal:
833 #endif
834 acceptState = state;
835 acceptPoint = p;
836 acceptLen = len;
837 if (c == '0') {
838 numTrailZeros++;
839 state = DECIMAL;
840 break;
841 } else if (isdigit(UCHAR(c))) {
842 if (objPtr != NULL) {
843 significandOverflow = AccumulateDecimalDigit(
844 (unsigned)(c - '0'), numTrailZeros,
845 &significandWide, &significandBig,
846 significandOverflow);
847 }
848 numSigDigs += numTrailZeros+1;
849 numTrailZeros = 0;
850 state = DECIMAL;
851 break;
852 } else if (flags & TCL_PARSE_INTEGER_ONLY) {
853 goto endgame;
854 } else if (c == '.') {
855 state = FRACTION;
856 break;
857 } else if (c == 'E' || c == 'e') {
858 state = EXPONENT_START;
859 break;
860 }
861 goto endgame;
862
863 /*
864 * Found a decimal point. If no digits have yet been scanned, E is
865 * not allowed; otherwise, it introduces the exponent. If at least
866 * one digit has been found, we have a possible complete number.
867 */
868
869 case FRACTION:
870 acceptState = state;
871 acceptPoint = p;
872 acceptLen = len;
873 if (c == 'E' || c=='e') {
874 state = EXPONENT_START;
875 break;
876 }
877 /* FALLTHROUGH */
878
879 case LEADING_RADIX_POINT:
880 if (c == '0') {
881 numDigitsAfterDp++;
882 numTrailZeros++;
883 state = FRACTION;
884 break;
885 } else if (isdigit(UCHAR(c))) {
886 numDigitsAfterDp++;
887 if (objPtr != NULL) {
888 significandOverflow = AccumulateDecimalDigit(
889 (unsigned)(c-'0'), numTrailZeros,
890 &significandWide, &significandBig,
891 significandOverflow);
892 }
893 if (numSigDigs != 0) {
894 numSigDigs += numTrailZeros+1;
895 } else {
896 numSigDigs = 1;
897 }
898 numTrailZeros = 0;
899 state = FRACTION;
900 break;
901 }
902 goto endgame;
903
904 case EXPONENT_START:
905 /*
906 * Scanned the E at the start of an exponent. Make sure a legal
907 * character follows before using the C library strtol routine,
908 * which allows whitespace.
909 */
910
911 if (c == '+') {
912 state = EXPONENT_SIGNUM;
913 break;
914 } else if (c == '-') {
915 exponentSignum = 1;
916 state = EXPONENT_SIGNUM;
917 break;
918 }
919 /* FALLTHROUGH */
920
921 case EXPONENT_SIGNUM:
922 /*
923 * Found the E at the start of the exponent, followed by a sign
924 * character.
925 */
926
927 if (isdigit(UCHAR(c))) {
928 exponent = c - '0';
929 state = EXPONENT;
930 break;
931 }
932 goto endgame;
933
934 case EXPONENT:
935 /*
936 * Found an exponent with at least one digit. Accumulate it,
937 * making sure to hard-pin it to LONG_MAX on overflow.
938 */
939
940 acceptState = state;
941 acceptPoint = p;
942 acceptLen = len;
943 if (isdigit(UCHAR(c))) {
944 if (exponent < (LONG_MAX - 9) / 10) {
945 exponent = 10 * exponent + (c - '0');
946 } else {
947 exponent = LONG_MAX;
948 }
949 state = EXPONENT;
950 break;
951 }
952 goto endgame;
953
954 /*
955 * Parse out INFINITY by simply spelling it out. INF is accepted
956 * as an abbreviation; other prefices are not.
957 */
958
959 case sI:
960 if (c == 'n' || c == 'N') {
961 state = sIN;
962 break;
963 }
964 goto endgame;
965 case sIN:
966 if (c == 'f' || c == 'F') {
967 state = sINF;
968 break;
969 }
970 goto endgame;
971 case sINF:
972 acceptState = state;
973 acceptPoint = p;
974 acceptLen = len;
975 if (c == 'i' || c == 'I') {
976 state = sINFI;
977 break;
978 }
979 goto endgame;
980 case sINFI:
981 if (c == 'n' || c == 'N') {
982 state = sINFIN;
983 break;
984 }
985 goto endgame;
986 case sINFIN:
987 if (c == 'i' || c == 'I') {
988 state = sINFINI;
989 break;
990 }
991 goto endgame;
992 case sINFINI:
993 if (c == 't' || c == 'T') {
994 state = sINFINIT;
995 break;
996 }
997 goto endgame;
998 case sINFINIT:
999 if (c == 'y' || c == 'Y') {
1000 state = sINFINITY;
1001 break;
1002 }
1003 goto endgame;
1004
1005 /*
1006 * Parse NaN's.
1007 */
1008 #ifdef IEEE_FLOATING_POINT
1009 case sN:
1010 if (c == 'a' || c == 'A') {
1011 state = sNA;
1012 break;
1013 }
1014 goto endgame;
1015 case sNA:
1016 if (c == 'n' || c == 'N') {
1017 state = sNAN;
1018 break;
1019 }
1020 goto endgame;
1021 case sNAN:
1022 acceptState = state;
1023 acceptPoint = p;
1024 acceptLen = len;
1025 if (c == '(') {
1026 state = sNANPAREN;
1027 break;
1028 }
1029 goto endgame;
1030
1031 /*
1032 * Parse NaN(hexdigits)
1033 */
1034 case sNANHEX:
1035 if (c == ')') {
1036 state = sNANFINISH;
1037 break;
1038 }
1039 /* FALLTHROUGH */
1040 case sNANPAREN:
1041 if (TclIsSpaceProc(c)) {
1042 break;
1043 }
1044 if (numSigDigs < 13) {
1045 if (c >= '0' && c <= '9') {
1046 d = c - '0';
1047 } else if (c >= 'a' && c <= 'f') {
1048 d = 10 + c - 'a';
1049 } else if (c >= 'A' && c <= 'F') {
1050 d = 10 + c - 'A';
1051 } else {
1052 goto endgame;
1053 }
1054 numSigDigs++;
1055 significandWide = (significandWide << 4) + d;
1056 state = sNANHEX;
1057 break;
1058 }
1059 goto endgame;
1060 case sNANFINISH:
1061 #endif
1062
1063 case sINFINITY:
1064 acceptState = state;
1065 acceptPoint = p;
1066 acceptLen = len;
1067 goto endgame;
1068 }
1069 p++;
1070 len--;
1071 }
1072
1073 endgame:
1074 if (acceptState == INITIAL) {
1075 /*
1076 * No numeric string at all found.
1077 */
1078
1079 status = TCL_ERROR;
1080 if (endPtrPtr != NULL) {
1081 *endPtrPtr = p;
1082 }
1083 } else {
1084 /*
1085 * Back up to the last accepting state in the lexer.
1086 */
1087
1088 p = acceptPoint;
1089 len = acceptLen;
1090 if (!(flags & TCL_PARSE_NO_WHITESPACE)) {
1091 /*
1092 * Accept trailing whitespace.
1093 */
1094
1095 while (len != 0 && TclIsSpaceProc(*p)) {
1096 p++;
1097 len--;
1098 }
1099 }
1100 if (endPtrPtr == NULL) {
1101 if ((len != 0) && ((numBytes > 0) || (*p != '\0'))) {
1102 status = TCL_ERROR;
1103 }
1104 } else {
1105 *endPtrPtr = p;
1106 }
1107 }
1108
1109 /*
1110 * Generate and store the appropriate internal rep.
1111 */
1112
1113 if (status == TCL_OK && objPtr != NULL) {
1114 TclFreeIntRep(objPtr);
1115 switch (acceptState) {
1116 case SIGNUM:
1117 case BAD_OCTAL:
1118 case ZERO_X:
1119 case ZERO_O:
1120 case ZERO_B:
1121 case LEADING_RADIX_POINT:
1122 case EXPONENT_START:
1123 case EXPONENT_SIGNUM:
1124 case sI:
1125 case sIN:
1126 case sINFI:
1127 case sINFIN:
1128 case sINFINI:
1129 case sINFINIT:
1130 #ifdef IEEE_FLOATING_POINT
1131 case sN:
1132 case sNA:
1133 case sNANPAREN:
1134 case sNANHEX:
1135 Tcl_Panic("TclParseNumber: bad acceptState %d parsing '%s'",
1136 acceptState, bytes);
1137 #endif
1138 case BINARY:
1139 shift = numTrailZeros;
1140 if (!significandOverflow && significandWide != 0 &&
1141 ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
1142 significandWide > (MOST_BITS + signum) >> shift)) {
1143 significandOverflow = 1;
1144 TclBNInitBignumFromWideUInt(&significandBig, significandWide);
1145 }
1146 if (shift) {
1147 if (!significandOverflow) {
1148 significandWide <<= shift;
1149 } else {
1150 mp_mul_2d(&significandBig, shift, &significandBig);
1151 }
1152 }
1153 goto returnInteger;
1154
1155 case HEXADECIMAL:
1156 /*
1157 * Returning a hex integer. Final scaling step.
1158 */
1159
1160 shift = 4 * numTrailZeros;
1161 if (!significandOverflow && significandWide !=0 &&
1162 ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
1163 significandWide > (MOST_BITS + signum) >> shift)) {
1164 significandOverflow = 1;
1165 TclBNInitBignumFromWideUInt(&significandBig, significandWide);
1166 }
1167 if (shift) {
1168 if (!significandOverflow) {
1169 significandWide <<= shift;
1170 } else {
1171 mp_mul_2d(&significandBig, shift, &significandBig);
1172 }
1173 }
1174 goto returnInteger;
1175
1176 case OCTAL:
1177 /*
1178 * Returning an octal integer. Final scaling step
1179 */
1180
1181 shift = 3 * numTrailZeros;
1182 if (!octalSignificandOverflow && octalSignificandWide != 0 &&
1183 ((size_t)shift >= CHAR_BIT*sizeof(Tcl_WideUInt) ||
1184 octalSignificandWide > (MOST_BITS + signum) >> shift)) {
1185 octalSignificandOverflow = 1;
1186 TclBNInitBignumFromWideUInt(&octalSignificandBig,
1187 octalSignificandWide);
1188 }
1189 if (shift) {
1190 if (!octalSignificandOverflow) {
1191 octalSignificandWide <<= shift;
1192 } else {
1193 mp_mul_2d(&octalSignificandBig, shift,
1194 &octalSignificandBig);
1195 }
1196 }
1197 if (!octalSignificandOverflow) {
1198 if (octalSignificandWide >
1199 (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
1200 #ifndef NO_WIDE_TYPE
1201 if (octalSignificandWide <= (MOST_BITS + signum)) {
1202 objPtr->typePtr = &tclWideIntType;
1203 if (signum) {
1204 objPtr->internalRep.wideValue =
1205 - (Tcl_WideInt) octalSignificandWide;
1206 } else {
1207 objPtr->internalRep.wideValue =
1208 (Tcl_WideInt) octalSignificandWide;
1209 }
1210 break;
1211 }
1212 #endif
1213 TclBNInitBignumFromWideUInt(&octalSignificandBig,
1214 octalSignificandWide);
1215 octalSignificandOverflow = 1;
1216 } else {
1217 objPtr->typePtr = &tclIntType;
1218 if (signum) {
1219 objPtr->internalRep.longValue =
1220 - (long) octalSignificandWide;
1221 } else {
1222 objPtr->internalRep.longValue =
1223 (long) octalSignificandWide;
1224 }
1225 }
1226 }
1227 if (octalSignificandOverflow) {
1228 if (signum) {
1229 mp_neg(&octalSignificandBig, &octalSignificandBig);
1230 }
1231 TclSetBignumIntRep(objPtr, &octalSignificandBig);
1232 }
1233 break;
1234
1235 case ZERO:
1236 case DECIMAL:
1237 significandOverflow = AccumulateDecimalDigit(0, numTrailZeros-1,
1238 &significandWide, &significandBig, significandOverflow);
1239 if (!significandOverflow && (significandWide > MOST_BITS+signum)) {
1240 significandOverflow = 1;
1241 TclBNInitBignumFromWideUInt(&significandBig, significandWide);
1242 }
1243 returnInteger:
1244 if (!significandOverflow) {
1245 if (significandWide >
1246 (Tcl_WideUInt)(((~(unsigned long)0) >> 1) + signum)) {
1247 #ifndef NO_WIDE_TYPE
1248 if (significandWide <= MOST_BITS+signum) {
1249 objPtr->typePtr = &tclWideIntType;
1250 if (signum) {
1251 objPtr->internalRep.wideValue =
1252 - (Tcl_WideInt) significandWide;
1253 } else {
1254 objPtr->internalRep.wideValue =
1255 (Tcl_WideInt) significandWide;
1256 }
1257 break;
1258 }
1259 #endif
1260 TclBNInitBignumFromWideUInt(&significandBig,
1261 significandWide);
1262 significandOverflow = 1;
1263 } else {
1264 objPtr->typePtr = &tclIntType;
1265 if (signum) {
1266 objPtr->internalRep.longValue =
1267 - (long) significandWide;
1268 } else {
1269 objPtr->internalRep.longValue =
1270 (long) significandWide;
1271 }
1272 }
1273 }
1274 if (significandOverflow) {
1275 if (signum) {
1276 mp_neg(&significandBig, &significandBig);
1277 }
1278 TclSetBignumIntRep(objPtr, &significandBig);
1279 }
1280 break;
1281
1282 case FRACTION:
1283 case EXPONENT:
1284
1285 /*
1286 * Here, we're parsing a floating-point number. 'significandWide'
1287 * or 'significandBig' contains the exact significand, according
1288 * to whether 'significandOverflow' is set. The desired floating
1289 * point value is significand * 10**k, where
1290 * k = numTrailZeros+exponent-numDigitsAfterDp.
1291 */
1292
1293 objPtr->typePtr = &tclDoubleType;
1294 if (exponentSignum) {
1295 exponent = - exponent;
1296 }
1297 if (!significandOverflow) {
1298 objPtr->internalRep.doubleValue = MakeLowPrecisionDouble(
1299 signum, significandWide, numSigDigs,
1300 (numTrailZeros + exponent - numDigitsAfterDp));
1301 } else {
1302 objPtr->internalRep.doubleValue = MakeHighPrecisionDouble(
1303 signum, &significandBig, numSigDigs,
1304 (numTrailZeros + exponent - numDigitsAfterDp));
1305 }
1306 break;
1307
1308 case sINF:
1309 case sINFINITY:
1310 if (signum) {
1311 objPtr->internalRep.doubleValue = -HUGE_VAL;
1312 } else {
1313 objPtr->internalRep.doubleValue = HUGE_VAL;
1314 }
1315 objPtr->typePtr = &tclDoubleType;
1316 break;
1317
1318 #ifdef IEEE_FLOATING_POINT
1319 case sNAN:
1320 case sNANFINISH:
1321 objPtr->internalRep.doubleValue = MakeNaN(signum, significandWide);
1322 objPtr->typePtr = &tclDoubleType;
1323 break;
1324 #endif
1325 case INITIAL:
1326 /* This case only to silence compiler warning */
1327 Tcl_Panic("TclParseNumber: state INITIAL can't happen here");
1328 }
1329 }
1330
1331 /*
1332 * Format an error message when an invalid number is encountered.
1333 */
1334
1335 if (status != TCL_OK) {
1336 if (interp != NULL) {
1337 Tcl_Obj *msg;
1338
1339 TclNewLiteralStringObj(msg, "expected ");
1340 Tcl_AppendToObj(msg, expected, -1);
1341 Tcl_AppendToObj(msg, " but got \"", -1);
1342 Tcl_AppendLimitedToObj(msg, bytes, numBytes, 50, "");
1343 Tcl_AppendToObj(msg, "\"", -1);
1344 if (state == BAD_OCTAL) {
1345 Tcl_AppendToObj(msg, " (looks like invalid octal number)", -1);
1346 }
1347 Tcl_SetObjResult(interp, msg);
1348 Tcl_SetErrorCode(interp, "TCL", "VALUE", "NUMBER", NULL);
1349 }
1350 }
1351
1352 /*
1353 * Free memory.
1354 */
1355
1356 if (octalSignificandOverflow) {
1357 mp_clear(&octalSignificandBig);
1358 }
1359 if (significandOverflow) {
1360 mp_clear(&significandBig);
1361 }
1362 return status;
1363 }
1364
1365 /*
1366 *----------------------------------------------------------------------
1367 *
1368 * AccumulateDecimalDigit --
1369 *
1370 * Consume a decimal digit in a number being scanned.
1371 *
1372 * Results:
1373 * Returns 1 if the number has overflowed to a bignum, 0 if it still fits
1374 * in a wide integer.
1375 *
1376 * Side effects:
1377 * Updates either the wide or bignum representation.
1378 *
1379 *----------------------------------------------------------------------
1380 */
1381
1382 static int
AccumulateDecimalDigit(unsigned digit,int numZeros,Tcl_WideUInt * wideRepPtr,mp_int * bignumRepPtr,int bignumFlag)1383 AccumulateDecimalDigit(
1384 unsigned digit, /* Digit being scanned. */
1385 int numZeros, /* Count of zero digits preceding the digit
1386 * being scanned. */
1387 Tcl_WideUInt *wideRepPtr, /* Representation of the partial number as a
1388 * wide integer. */
1389 mp_int *bignumRepPtr, /* Representation of the partial number as a
1390 * bignum. */
1391 int bignumFlag) /* Flag == 1 if the number overflowed previous
1392 * to this digit. */
1393 {
1394 int i, n;
1395 Tcl_WideUInt w;
1396
1397 /*
1398 * Try wide multiplication first
1399 */
1400
1401 if (!bignumFlag) {
1402 w = *wideRepPtr;
1403 if (w == 0) {
1404 /*
1405 * There's no need to multiply if the multiplicand is zero.
1406 */
1407
1408 *wideRepPtr = digit;
1409 return 0;
1410 } else if (numZeros >= maxpow10_wide
1411 || w > ((~(Tcl_WideUInt)0)-digit)/pow10_wide[numZeros+1]) {
1412 /*
1413 * Wide multiplication will overflow. Expand the
1414 * number to a bignum and fall through into the bignum case.
1415 */
1416
1417 TclBNInitBignumFromWideUInt(bignumRepPtr, w);
1418 } else {
1419 /*
1420 * Wide multiplication.
1421 */
1422 *wideRepPtr = w * pow10_wide[numZeros+1] + digit;
1423 return 0;
1424 }
1425 }
1426
1427 /*
1428 * Bignum multiplication.
1429 */
1430
1431 if (numZeros < log10_DIGIT_MAX) {
1432 /*
1433 * Up to about 8 zeros - single digit multiplication.
1434 */
1435
1436 mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[numZeros+1],
1437 bignumRepPtr);
1438 mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
1439 } else {
1440 /*
1441 * More than single digit multiplication. Multiply by the appropriate
1442 * small powers of 5, and then shift. Large strings of zeroes are
1443 * eaten 256 at a time; this is less efficient than it could be, but
1444 * seems implausible. We presume that DIGIT_BIT is at least 27. The
1445 * first multiplication, by up to 10**7, is done with a one-DIGIT
1446 * multiply (this presumes that DIGIT_BIT >= 24).
1447 */
1448
1449 n = numZeros + 1;
1450 mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
1451 for (i=3; i<=7; ++i) {
1452 if (n & (1 << i)) {
1453 mp_mul(bignumRepPtr, pow5+i, bignumRepPtr);
1454 }
1455 }
1456 while (n >= 256) {
1457 mp_mul(bignumRepPtr, pow5+8, bignumRepPtr);
1458 n -= 256;
1459 }
1460 mp_mul_2d(bignumRepPtr, (int)(numZeros+1)&~0x7, bignumRepPtr);
1461 mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
1462 }
1463
1464 return 1;
1465 }
1466
1467 /*
1468 *----------------------------------------------------------------------
1469 *
1470 * MakeLowPrecisionDouble --
1471 *
1472 * Makes the double precision number, signum*significand*10**exponent.
1473 *
1474 * Results:
1475 * Returns the constructed number.
1476 *
1477 * Common cases, where there are few enough digits that the number can be
1478 * represented with at most roundoff, are handled specially here. If the
1479 * number requires more than one rounded operation to compute, the code
1480 * promotes the significand to a bignum and calls MakeHighPrecisionDouble
1481 * to do it instead.
1482 *
1483 *----------------------------------------------------------------------
1484 */
1485
1486 static double
MakeLowPrecisionDouble(int signum,Tcl_WideUInt significand,int numSigDigs,int exponent)1487 MakeLowPrecisionDouble(
1488 int signum, /* 1 if the number is negative, 0 otherwise */
1489 Tcl_WideUInt significand, /* Significand of the number */
1490 int numSigDigs, /* Number of digits in the significand */
1491 int exponent) /* Power of ten */
1492 {
1493 double retval; /* Value of the number */
1494 mp_int significandBig; /* Significand expressed as a bignum */
1495
1496 /*
1497 * With gcc on x86, the floating point rounding mode is double-extended.
1498 * This causes the result of double-precision calculations to be rounded
1499 * twice: once to the precision of double-extended and then again to the
1500 * precision of double. Double-rounding introduces gratuitous errors of 1
1501 * ulp, so we need to change rounding mode to 53-bits.
1502 */
1503
1504 #if defined(__GNUC__) && defined(__i386)
1505 fpu_control_t roundTo53Bits = 0x027f;
1506 fpu_control_t oldRoundingMode;
1507 _FPU_GETCW(oldRoundingMode);
1508 _FPU_SETCW(roundTo53Bits);
1509 #endif
1510 #if defined(__sun) && defined(__i386) && !defined(__GNUC__)
1511 ieee_flags("set","precision","double",NULL);
1512 #endif
1513
1514 /*
1515 * Test for the easy cases.
1516 */
1517
1518 if (numSigDigs <= QUICK_MAX) {
1519 if (exponent >= 0) {
1520 if (exponent <= mmaxpow) {
1521 /*
1522 * The significand is an exact integer, and so is
1523 * 10**exponent. The product will be correct to within 1/2 ulp
1524 * without special handling.
1525 */
1526
1527 retval = (double)(Tcl_WideInt)significand * pow10vals[exponent];
1528 goto returnValue;
1529 } else {
1530 int diff = QUICK_MAX - numSigDigs;
1531 if (exponent-diff <= mmaxpow) {
1532 /*
1533 * 10**exponent is not an exact integer, but
1534 * 10**(exponent-diff) is exact, and so is
1535 * significand*10**diff, so we can still compute the value
1536 * with only one roundoff.
1537 */
1538
1539 volatile double factor =
1540 (double)(Tcl_WideInt)significand * pow10vals[diff];
1541 retval = factor * pow10vals[exponent-diff];
1542 goto returnValue;
1543 }
1544 }
1545 } else {
1546 if (exponent >= -mmaxpow) {
1547 /*
1548 * 10**-exponent is an exact integer, and so is the
1549 * significand. Compute the result by one division, again with
1550 * only one rounding.
1551 */
1552
1553 retval = (double)(Tcl_WideInt)significand / pow10vals[-exponent];
1554 goto returnValue;
1555 }
1556 }
1557 }
1558
1559 /*
1560 * All the easy cases have failed. Promote ths significand to bignum and
1561 * call MakeHighPrecisionDouble to do it the hard way.
1562 */
1563
1564 TclBNInitBignumFromWideUInt(&significandBig, significand);
1565 retval = MakeHighPrecisionDouble(0, &significandBig, numSigDigs,
1566 exponent);
1567 mp_clear(&significandBig);
1568
1569 /*
1570 * Come here to return the computed value.
1571 */
1572
1573 returnValue:
1574 if (signum) {
1575 retval = -retval;
1576 }
1577
1578 /*
1579 * On gcc on x86, restore the floating point mode word.
1580 */
1581
1582 #if defined(__GNUC__) && defined(__i386)
1583 _FPU_SETCW(oldRoundingMode);
1584 #endif
1585 #if defined(__sun) && defined(__i386) && !defined(__GNUC__)
1586 ieee_flags("clear","precision",NULL,NULL);
1587 #endif
1588
1589 return retval;
1590 }
1591
1592 /*
1593 *----------------------------------------------------------------------
1594 *
1595 * MakeHighPrecisionDouble --
1596 *
1597 * Makes the double precision number, signum*significand*10**exponent.
1598 *
1599 * Results:
1600 * Returns the constructed number.
1601 *
1602 * MakeHighPrecisionDouble is used when arbitrary-precision arithmetic is
1603 * needed to ensure correct rounding. It begins by calculating a
1604 * low-precision approximation to the desired number, and then refines
1605 * the answer in high precision.
1606 *
1607 *----------------------------------------------------------------------
1608 */
1609
1610 static double
MakeHighPrecisionDouble(int signum,mp_int * significand,int numSigDigs,int exponent)1611 MakeHighPrecisionDouble(
1612 int signum, /* 1=negative, 0=nonnegative */
1613 mp_int *significand, /* Exact significand of the number */
1614 int numSigDigs, /* Number of significant digits */
1615 int exponent) /* Power of 10 by which to multiply */
1616 {
1617 double retval;
1618 int machexp; /* Machine exponent of a power of 10 */
1619
1620 /*
1621 * With gcc on x86, the floating point rounding mode is double-extended.
1622 * This causes the result of double-precision calculations to be rounded
1623 * twice: once to the precision of double-extended and then again to the
1624 * precision of double. Double-rounding introduces gratuitous errors of 1
1625 * ulp, so we need to change rounding mode to 53-bits.
1626 */
1627
1628 #if defined(__GNUC__) && defined(__i386)
1629 fpu_control_t roundTo53Bits = 0x027f;
1630 fpu_control_t oldRoundingMode;
1631 _FPU_GETCW(oldRoundingMode);
1632 _FPU_SETCW(roundTo53Bits);
1633 #endif
1634 #if defined(__sun) && defined(__i386) && !defined(__GNUC__)
1635 ieee_flags("set","precision","double",NULL);
1636 #endif
1637
1638 /*
1639 * Quick checks for over/underflow.
1640 */
1641
1642 if (numSigDigs+exponent-1 > maxDigits) {
1643 retval = HUGE_VAL;
1644 goto returnValue;
1645 }
1646 if (numSigDigs+exponent-1 < minDigits) {
1647 retval = 0;
1648 goto returnValue;
1649 }
1650
1651 /*
1652 * Develop a first approximation to the significand. It is tempting simply
1653 * to force bignum to double, but that will overflow on input numbers like
1654 * 1.[string repeat 0 1000]1; while this is a not terribly likely
1655 * scenario, we still have to deal with it. Use fraction and exponent
1656 * instead. Once we have the significand, multiply by 10**exponent. Test
1657 * for overflow. Convert back to a double, and test for underflow.
1658 */
1659
1660 retval = BignumToBiasedFrExp(significand, &machexp);
1661 retval = Pow10TimesFrExp(exponent, retval, &machexp);
1662 if (machexp > DBL_MAX_EXP*log2FLT_RADIX) {
1663 retval = HUGE_VAL;
1664 goto returnValue;
1665 }
1666 retval = SafeLdExp(retval, machexp);
1667 if (tiny == 0.0) {
1668 tiny = SafeLdExp(1.0, DBL_MIN_EXP * log2FLT_RADIX - mantBits);
1669 }
1670 if (retval < tiny) {
1671 retval = tiny;
1672 }
1673
1674 /*
1675 * Refine the result twice. (The second refinement should be necessary
1676 * only if the best approximation is a power of 2 minus 1/2 ulp).
1677 */
1678
1679 retval = RefineApproximation(retval, significand, exponent);
1680 retval = RefineApproximation(retval, significand, exponent);
1681
1682 /*
1683 * Come here to return the computed value.
1684 */
1685
1686 returnValue:
1687 if (signum) {
1688 retval = -retval;
1689 }
1690
1691 /*
1692 * On gcc on x86, restore the floating point mode word.
1693 */
1694
1695 #if defined(__GNUC__) && defined(__i386)
1696 _FPU_SETCW(oldRoundingMode);
1697 #endif
1698 #if defined(__sun) && defined(__i386) && !defined(__GNUC__)
1699 ieee_flags("clear","precision",NULL,NULL);
1700 #endif
1701 return retval;
1702 }
1703
1704 /*
1705 *----------------------------------------------------------------------
1706 *
1707 * MakeNaN --
1708 *
1709 * Makes a "Not a Number" given a set of bits to put in the tag bits
1710 *
1711 * Note that a signalling NaN is never returned.
1712 *
1713 *----------------------------------------------------------------------
1714 */
1715
1716 #ifdef IEEE_FLOATING_POINT
1717 static double
MakeNaN(int signum,Tcl_WideUInt tags)1718 MakeNaN(
1719 int signum, /* Sign bit (1=negative, 0=nonnegative */
1720 Tcl_WideUInt tags) /* Tag bits to put in the NaN */
1721 {
1722 union {
1723 Tcl_WideUInt iv;
1724 double dv;
1725 } theNaN;
1726
1727 theNaN.iv = tags;
1728 theNaN.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
1729 if (signum) {
1730 theNaN.iv |= ((Tcl_WideUInt) (0x8000 | NAN_START)) << 48;
1731 } else {
1732 theNaN.iv |= ((Tcl_WideUInt) NAN_START) << 48;
1733 }
1734 if (n770_fp) {
1735 theNaN.iv = Nokia770Twiddle(theNaN.iv);
1736 }
1737 return theNaN.dv;
1738 }
1739 #endif
1740
1741 /*
1742 *----------------------------------------------------------------------
1743 *
1744 * RefineApproximation --
1745 *
1746 * Given a poor approximation to a floating point number, returns a
1747 * better one. (The better approximation is correct to within 1 ulp, and
1748 * is entirely correct if the poor approximation is correct to 1 ulp.)
1749 *
1750 * Results:
1751 * Returns the improved result.
1752 *
1753 *----------------------------------------------------------------------
1754 */
1755
1756 static double
RefineApproximation(double approxResult,mp_int * exactSignificand,int exponent)1757 RefineApproximation(
1758 double approxResult, /* Approximate result of conversion */
1759 mp_int *exactSignificand, /* Integer significand */
1760 int exponent) /* Power of 10 to multiply by significand */
1761 {
1762 int M2, M5; /* Powers of 2 and of 5 needed to put the
1763 * decimal and binary numbers over a common
1764 * denominator. */
1765 double significand; /* Sigificand of the binary number */
1766 int binExponent; /* Exponent of the binary number */
1767 int msb; /* Most significant bit position of an
1768 * intermediate result */
1769 int nDigits; /* Number of mp_digit's in an intermediate
1770 * result */
1771 mp_int twoMv; /* Approx binary value expressed as an exact
1772 * integer scaled by the multiplier 2M */
1773 mp_int twoMd; /* Exact decimal value expressed as an exact
1774 * integer scaled by the multiplier 2M */
1775 int scale; /* Scale factor for M */
1776 int multiplier; /* Power of two to scale M */
1777 double num, den; /* Numerator and denominator of the correction
1778 * term */
1779 double quot; /* Correction term */
1780 double minincr; /* Lower bound on the absolute value of the
1781 * correction term. */
1782 int roundToEven = 0; /* Flag == TRUE if we need to invoke
1783 * "round to even" functionality */
1784 double rteSignificand; /* Significand of the round-to-even result */
1785 int rteExponent; /* Exponent of the round-to-even result */
1786 Tcl_WideInt rteSigWide; /* Wide integer version of the significand
1787 * for testing evenness */
1788 int i;
1789
1790 /*
1791 * The first approximation is always low. If we find that it's HUGE_VAL,
1792 * we're done.
1793 */
1794
1795 if (approxResult == HUGE_VAL) {
1796 return approxResult;
1797 }
1798
1799 /*
1800 * Find a common denominator for the decimal and binary fractions. The
1801 * common denominator will be 2**M2 + 5**M5.
1802 */
1803
1804 significand = frexp(approxResult, &binExponent);
1805 i = mantBits - binExponent;
1806 if (i < 0) {
1807 M2 = 0;
1808 } else {
1809 M2 = i;
1810 }
1811 if (exponent > 0) {
1812 M5 = 0;
1813 } else {
1814 M5 = -exponent;
1815 if ((M5-1) > M2) {
1816 M2 = M5-1;
1817 }
1818 }
1819
1820 /*
1821 * The floating point number is significand*2**binExponent. Compute the
1822 * large integer significand*2**(binExponent+M2+1). The 2**-1 bit of the
1823 * significand (the most significant) corresponds to the
1824 * 2**(binExponent+M2 + 1) bit of 2*M2*v. Allocate enough digits to hold
1825 * that quantity, then convert the significand to a large integer, scaled
1826 * appropriately. Then multiply by the appropriate power of 5.
1827 */
1828
1829 msb = binExponent + M2; /* 1008 */
1830 nDigits = msb / DIGIT_BIT + 1;
1831 mp_init_size(&twoMv, nDigits);
1832 i = (msb % DIGIT_BIT + 1);
1833 twoMv.used = nDigits;
1834 significand *= SafeLdExp(1.0, i);
1835 while (--nDigits >= 0) {
1836 twoMv.dp[nDigits] = (mp_digit) significand;
1837 significand -= (mp_digit) significand;
1838 significand = SafeLdExp(significand, DIGIT_BIT);
1839 }
1840 for (i = 0; i <= 8; ++i) {
1841 if (M5 & (1 << i)) {
1842 mp_mul(&twoMv, pow5+i, &twoMv);
1843 }
1844 }
1845
1846 /*
1847 * Collect the decimal significand as a high precision integer. The least
1848 * significant bit corresponds to bit M2+exponent+1 so it will need to be
1849 * shifted left by that many bits after being multiplied by
1850 * 5**(M5+exponent).
1851 */
1852
1853 mp_init_copy(&twoMd, exactSignificand);
1854 for (i=0; i<=8; ++i) {
1855 if ((M5+exponent) & (1 << i)) {
1856 mp_mul(&twoMd, pow5+i, &twoMd);
1857 }
1858 }
1859 mp_mul_2d(&twoMd, M2+exponent+1, &twoMd);
1860 mp_sub(&twoMd, &twoMv, &twoMd);
1861
1862 /*
1863 * The result, 2Mv-2Md, needs to be divided by 2M to yield a correction
1864 * term. Because 2M may well overflow a double, we need to scale the
1865 * denominator by a factor of 2**binExponent-mantBits
1866 */
1867
1868 scale = binExponent - mantBits - 1;
1869
1870 mp_set(&twoMv, 1);
1871 for (i=0; i<=8; ++i) {
1872 if (M5 & (1 << i)) {
1873 mp_mul(&twoMv, pow5+i, &twoMv);
1874 }
1875 }
1876 multiplier = M2 + scale + 1;
1877 if (multiplier > 0) {
1878 mp_mul_2d(&twoMv, multiplier, &twoMv);
1879 } else if (multiplier < 0) {
1880 mp_div_2d(&twoMv, -multiplier, &twoMv, NULL);
1881 }
1882
1883 switch (mp_cmp_mag(&twoMd, &twoMv)) {
1884 case MP_LT:
1885 /*
1886 * If the result is less than unity, the error is less than 1/2 unit in
1887 * the last place, so there's no correction to make.
1888 */
1889 mp_clear(&twoMd);
1890 mp_clear(&twoMv);
1891 return approxResult;
1892 case MP_EQ:
1893 /*
1894 * If the result is exactly unity, we need to round to even.
1895 */
1896 roundToEven = 1;
1897 break;
1898 case MP_GT:
1899 break;
1900 }
1901
1902 if (roundToEven) {
1903 rteSignificand = frexp(approxResult, &rteExponent);
1904 rteSigWide = (Tcl_WideInt) ldexp(rteSignificand, FP_PRECISION);
1905 if ((rteSigWide & 1) == 0) {
1906 mp_clear(&twoMd);
1907 mp_clear(&twoMv);
1908 return approxResult;
1909 }
1910 }
1911
1912 /*
1913 * Convert the numerator and denominator of the corrector term accurately
1914 * to floating point numbers.
1915 */
1916
1917 num = TclBignumToDouble(&twoMd);
1918 den = TclBignumToDouble(&twoMv);
1919
1920 quot = SafeLdExp(num/den, scale);
1921 minincr = SafeLdExp(1.0, binExponent-mantBits);
1922
1923 if (quot<0. && quot>-minincr) {
1924 quot = -minincr;
1925 } else if (quot>0. && quot<minincr) {
1926 quot = minincr;
1927 }
1928
1929 mp_clear(&twoMd);
1930 mp_clear(&twoMv);
1931
1932 return approxResult + quot;
1933 }
1934
1935 /*
1936 *-----------------------------------------------------------------------------
1937 *
1938 * MultPow5 --
1939 *
1940 * Multiply a bignum by a power of 5.
1941 *
1942 * Side effects:
1943 * Stores base*5**n in result
1944 *
1945 *-----------------------------------------------------------------------------
1946 */
1947
1948 inline static void
MulPow5(mp_int * base,unsigned n,mp_int * result)1949 MulPow5(mp_int* base, /* Number to multiply */
1950 unsigned n, /* Power of 5 to multiply by */
1951 mp_int* result) /* Place to store the result */
1952 {
1953 mp_int* p = base;
1954 int n13 = n / 13;
1955 int r = n % 13;
1956 if (r != 0) {
1957 mp_mul_d(p, dpow5[r], result);
1958 p = result;
1959 }
1960 r = 0;
1961 while (n13 != 0) {
1962 if (n13 & 1) {
1963 mp_mul(p, pow5_13+r, result);
1964 p = result;
1965 }
1966 n13 >>= 1;
1967 ++r;
1968 }
1969 if (p != result) {
1970 mp_copy(p, result);
1971 }
1972 }
1973
1974 /*
1975 *-----------------------------------------------------------------------------
1976 *
1977 * NormalizeRightward --
1978 *
1979 * Shifts a number rightward until it is odd (that is, until the
1980 * least significant bit is nonzero.
1981 *
1982 * Results:
1983 * Returns the number of bit positions by which the number was shifted.
1984 *
1985 * Side effects:
1986 * Shifts the number in place; *wPtr is replaced by the shifted number.
1987 *
1988 *-----------------------------------------------------------------------------
1989 */
1990
1991 inline static int
NormalizeRightward(Tcl_WideUInt * wPtr)1992 NormalizeRightward(Tcl_WideUInt* wPtr)
1993 /* INOUT: Number to shift */
1994 {
1995 int rv = 0;
1996 Tcl_WideUInt w = *wPtr;
1997 if (!(w & (Tcl_WideUInt) 0xffffffff)) {
1998 w >>= 32; rv += 32;
1999 }
2000 if (!(w & (Tcl_WideUInt) 0xffff)) {
2001 w >>= 16; rv += 16;
2002 }
2003 if (!(w & (Tcl_WideUInt) 0xff)) {
2004 w >>= 8; rv += 8;
2005 }
2006 if (!(w & (Tcl_WideUInt) 0xf)) {
2007 w >>= 4; rv += 4;
2008 }
2009 if (!(w & 0x3)) {
2010 w >>= 2; rv += 2;
2011 }
2012 if (!(w & 0x1)) {
2013 w >>= 1; ++rv;
2014 }
2015 *wPtr = w;
2016 return rv;
2017 }
2018
2019 /*
2020 *-----------------------------------------------------------------------------0
2021 *
2022 * RequiredPrecision --
2023 *
2024 * Determines the number of bits needed to hold an intger.
2025 *
2026 * Results:
2027 * Returns the position of the most significant bit (0 - 63).
2028 * Returns 0 if the number is zero.
2029 *
2030 *----------------------------------------------------------------------------
2031 */
2032
2033 static int
RequiredPrecision(Tcl_WideUInt w)2034 RequiredPrecision(Tcl_WideUInt w)
2035 /* Number to interrogate */
2036 {
2037 int rv;
2038 unsigned long wi;
2039 if (w & ((Tcl_WideUInt) 0xffffffff << 32)) {
2040 wi = (unsigned long) (w >> 32); rv = 32;
2041 } else {
2042 wi = (unsigned long) w; rv = 0;
2043 }
2044 if (wi & 0xffff0000) {
2045 wi >>= 16; rv += 16;
2046 }
2047 if (wi & 0xff00) {
2048 wi >>= 8; rv += 8;
2049 }
2050 if (wi & 0xf0) {
2051 wi >>= 4; rv += 4;
2052 }
2053 if (wi & 0xc) {
2054 wi >>= 2; rv += 2;
2055 }
2056 if (wi & 0x2) {
2057 wi >>= 1; ++rv;
2058 }
2059 if (wi & 0x1) {
2060 ++rv;
2061 }
2062 return rv;
2063 }
2064
2065 /*
2066 *-----------------------------------------------------------------------------
2067 *
2068 * DoubleToExpAndSig --
2069 *
2070 * Separates a 'double' into exponent and significand.
2071 *
2072 * Side effects:
2073 * Stores the significand in '*significand' and the exponent in
2074 * '*expon' so that dv == significand * 2.0**expon, and significand
2075 * is odd. Also stores the position of the leftmost 1-bit in 'significand'
2076 * in 'bits'.
2077 *
2078 *-----------------------------------------------------------------------------
2079 */
2080
2081 inline static void
DoubleToExpAndSig(double dv,Tcl_WideUInt * significand,int * expon,int * bits)2082 DoubleToExpAndSig(double dv, /* Number to convert */
2083 Tcl_WideUInt* significand,
2084 /* OUTPUT: Significand of the number */
2085 int* expon, /* OUTPUT: Exponent to multiply the number by */
2086 int* bits) /* OUTPUT: Number of significant bits */
2087 {
2088 Double d; /* Number being converted */
2089 Tcl_WideUInt z; /* Significand under construction */
2090 int de; /* Exponent of the number */
2091 int k; /* Bit count */
2092
2093 d.d = dv;
2094
2095 /* Extract exponent and significand */
2096
2097 de = (d.w.word0 & EXP_MASK) >> EXP_SHIFT;
2098 z = d.q & SIG_MASK;
2099 if (de != 0) {
2100 z |= HIDDEN_BIT;
2101 k = NormalizeRightward(&z);
2102 *bits = FP_PRECISION - k;
2103 *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1);
2104 } else {
2105 k = NormalizeRightward(&z);
2106 *expon = k + (de - EXPONENT_BIAS) - (FP_PRECISION-1) + 1;
2107 *bits = RequiredPrecision(z);
2108 }
2109 *significand = z;
2110 }
2111
2112 /*
2113 *-----------------------------------------------------------------------------
2114 *
2115 * TakeAbsoluteValue --
2116 *
2117 * Takes the absolute value of a 'double' including 0, Inf and NaN
2118 *
2119 * Side effects:
2120 * The 'double' in *d is replaced with its absolute value. The
2121 * signum is stored in 'sign': 1 for negative, 0 for nonnegative.
2122 *
2123 *-----------------------------------------------------------------------------
2124 */
2125
2126 inline static void
TakeAbsoluteValue(Double * d,int * sign)2127 TakeAbsoluteValue(Double* d, /* Number to replace with absolute value */
2128 int* sign) /* Place to put the signum */
2129 {
2130 if (d->w.word0 & SIGN_BIT) {
2131 *sign = 1;
2132 d->w.word0 &= ~SIGN_BIT;
2133 } else {
2134 *sign = 0;
2135 }
2136 }
2137
2138 /*
2139 *-----------------------------------------------------------------------------
2140 *
2141 * FormatInfAndNaN --
2142 *
2143 * Bailout for formatting infinities and Not-A-Number.
2144 *
2145 * Results:
2146 * Returns one of the strings 'Infinity' and 'NaN'.
2147 *
2148 * Side effects:
2149 * Stores 9999 in *decpt, and sets '*endPtr' to designate the
2150 * terminating NUL byte of the string if 'endPtr' is not NULL.
2151 *
2152 * The string returned must be freed by the caller using 'ckfree'.
2153 *
2154 *-----------------------------------------------------------------------------
2155 */
2156
2157 inline static char*
FormatInfAndNaN(Double * d,int * decpt,char ** endPtr)2158 FormatInfAndNaN(Double* d, /* Exceptional number to format */
2159 int* decpt, /* Decimal point to set to a bogus value */
2160 char** endPtr) /* Pointer to the end of the formatted data */
2161 {
2162 char* retval;
2163 *decpt = 9999;
2164 if (!(d->w.word1) && !(d->w.word0 & HI_ORDER_SIG_MASK)) {
2165 retval = ckalloc(9);
2166 strcpy(retval, "Infinity");
2167 if (endPtr) {
2168 *endPtr = retval + 8;
2169 }
2170 } else {
2171 retval = ckalloc(4);
2172 strcpy(retval, "NaN");
2173 if (endPtr) {
2174 *endPtr = retval + 3;
2175 }
2176 }
2177 return retval;
2178 }
2179
2180 /*
2181 *-----------------------------------------------------------------------------
2182 *
2183 * FormatZero --
2184 *
2185 * Bailout to format a zero floating-point number.
2186 *
2187 * Results:
2188 * Returns the constant string "0"
2189 *
2190 * Side effects:
2191 * Stores 1 in '*decpt' and puts a pointer to the NUL byte terminating
2192 * the string in '*endPtr' if 'endPtr' is not NULL.
2193 *
2194 *-----------------------------------------------------------------------------
2195 */
2196
2197 inline static char*
FormatZero(int * decpt,char ** endPtr)2198 FormatZero(int* decpt, /* Location of the decimal point */
2199 char** endPtr) /* Pointer to the end of the formatted data */
2200 {
2201 char* retval = ckalloc(2);
2202 strcpy(retval, "0");
2203 if (endPtr) {
2204 *endPtr = retval+1;
2205 }
2206 *decpt = 0;
2207 return retval;
2208 }
2209
2210 /*
2211 *-----------------------------------------------------------------------------
2212 *
2213 * ApproximateLog10 --
2214 *
2215 * Computes a two-term Taylor series approximation to the common
2216 * log of a number, and computes the number's binary log.
2217 *
2218 * Results:
2219 * Return an approximation to floor(log10(bw*2**be)) that is either
2220 * exact or 1 too high.
2221 *
2222 *-----------------------------------------------------------------------------
2223 */
2224
2225 inline static int
ApproximateLog10(Tcl_WideUInt bw,int be,int bbits)2226 ApproximateLog10(Tcl_WideUInt bw,
2227 /* Integer significand of the number */
2228 int be, /* Power of two to scale bw */
2229 int bbits) /* Number of bits of precision in bw */
2230 {
2231 int i; /* Log base 2 of the number */
2232 int k; /* Floor(Log base 10 of the number) */
2233 double ds; /* Mantissa of the number */
2234 Double d2;
2235
2236 /*
2237 * Compute i and d2 such that d = d2*2**i, and 1 < d2 < 2.
2238 * Compute an approximation to log10(d),
2239 * log10(d) ~ log10(2) * i + log10(1.5)
2240 * + (significand-1.5)/(1.5 * log(10))
2241 */
2242
2243 d2.q = bw << (FP_PRECISION - bbits) & SIG_MASK;
2244 d2.w.word0 |= (EXPONENT_BIAS) << EXP_SHIFT;
2245 i = be + bbits - 1;
2246 ds = (d2.d - 1.5) * TWO_OVER_3LOG10
2247 + LOG10_3HALVES_PLUS_FUDGE
2248 + LOG10_2 * i;
2249 k = (int) ds;
2250 if (k > ds) {
2251 --k;
2252 }
2253 return k;
2254 }
2255
2256 /*
2257 *-----------------------------------------------------------------------------
2258 *
2259 * BetterLog10 --
2260 *
2261 * Improves the result of ApproximateLog10 for numbers in the range
2262 * 1 .. 10**(TEN_PMAX)-1
2263 *
2264 * Side effects:
2265 * Sets k_check to 0 if the new result is known to be exact, and to
2266 * 1 if it may still be one too high.
2267 *
2268 * Results:
2269 * Returns the improved approximation to log10(d)
2270 *
2271 *-----------------------------------------------------------------------------
2272 */
2273
2274 inline static int
BetterLog10(double d,int k,int * k_check)2275 BetterLog10(double d, /* Original number to format */
2276 int k, /* Characteristic(Log base 10) of the number */
2277 int* k_check) /* Flag == 1 if k is inexact */
2278 {
2279 /*
2280 * Performance hack. If k is in the range 0..TEN_PMAX, then we can
2281 * use a powers-of-ten table to check it.
2282 */
2283 if (k >= 0 && k <= TEN_PMAX) {
2284 if (d < tens[k]) {
2285 k--;
2286 }
2287 *k_check = 0;
2288 } else {
2289 *k_check = 1;
2290 }
2291 return k;
2292 }
2293
2294 /*
2295 *-----------------------------------------------------------------------------
2296 *
2297 * ComputeScale --
2298 *
2299 * Prepares to format a floating-point number as decimal.
2300 *
2301 * Parameters:
2302 * floor(log10*x) is k (or possibly k-1). floor(log2(x) is i.
2303 * The significand of x requires bbits bits to represent.
2304 *
2305 * Results:
2306 * Determines integers b2, b5, s2, s5 so that sig*2**b2*5**b5/2**s2*2**s5
2307 * exactly represents the value of the x/10**k. This value will lie
2308 * in the range [1 .. 10), and allows for computing successive digits
2309 * by multiplying sig%10 by 10.
2310 *
2311 *-----------------------------------------------------------------------------
2312 */
2313
2314 inline static void
ComputeScale(int be,int k,int * b2,int * b5,int * s2,int * s5)2315 ComputeScale(int be, /* Exponent part of number: d = bw * 2**be */
2316 int k, /* Characteristic of log10(number) */
2317 int* b2, /* OUTPUT: Power of 2 in the numerator */
2318 int* b5, /* OUTPUT: Power of 5 in the numerator */
2319 int* s2, /* OUTPUT: Power of 2 in the denominator */
2320 int* s5) /* OUTPUT: Power of 5 in the denominator */
2321 {
2322
2323 /*
2324 * Scale numerator and denominator powers of 2 so that the
2325 * input binary number is the ratio of integers
2326 */
2327 if (be <= 0) {
2328 *b2 = 0;
2329 *s2 = -be;
2330 } else {
2331 *b2 = be;
2332 *s2 = 0;
2333 }
2334
2335 /*
2336 * Scale numerator and denominator so that the output decimal number
2337 * is the ratio of integers
2338 */
2339 if (k >= 0) {
2340 *b5 = 0;
2341 *s5 = k;
2342 *s2 += k;
2343 } else {
2344 *b2 -= k;
2345 *b5 = -k;
2346 *s5 = 0;
2347 }
2348 }
2349
2350 /*
2351 *-----------------------------------------------------------------------------
2352 *
2353 * SetPrecisionLimits --
2354 *
2355 * Determines how many digits of significance should be computed
2356 * (and, hence, how much memory need be allocated) for formatting a
2357 * floating point number.
2358 *
2359 * Given that 'k' is floor(log10(x)):
2360 * if 'shortest' format is used, there will be at most 18 digits in the result.
2361 * if 'F' format is used, there will be at most 'ndigits' + k + 1 digits
2362 * if 'E' format is used, there will be exactly 'ndigits' digits.
2363 *
2364 * Side effects:
2365 * Adjusts '*ndigitsPtr' to have a valid value.
2366 * Stores the maximum memory allocation needed in *iPtr.
2367 * Sets '*iLimPtr' to the limiting number of digits to convert if k
2368 * has been guessed correctly, and '*iLim1Ptr' to the limiting number
2369 * of digits to convert if k has been guessed to be one too high.
2370 *
2371 *-----------------------------------------------------------------------------
2372 */
2373
2374 inline static void
SetPrecisionLimits(int convType,int k,int * ndigitsPtr,int * iPtr,int * iLimPtr,int * iLim1Ptr)2375 SetPrecisionLimits(int convType,
2376 /* Type of conversion:
2377 * TCL_DD_SHORTEST
2378 * TCL_DD_STEELE0
2379 * TCL_DD_E_FMT
2380 * TCL_DD_F_FMT */
2381 int k, /* Floor(log10(number to convert)) */
2382 int* ndigitsPtr,
2383 /* IN/OUT: Number of digits requested
2384 * (Will be adjusted if needed) */
2385 int* iPtr, /* OUT: Maximum number of digits
2386 * to return */
2387 int *iLimPtr,/* OUT: Number of digits of significance
2388 * if the bignum method is used.*/
2389 int *iLim1Ptr)
2390 /* OUT: Number of digits of significance
2391 * if the quick method is used. */
2392 {
2393 switch(convType) {
2394 case TCL_DD_SHORTEST0:
2395 case TCL_DD_STEELE0:
2396 *iLimPtr = *iLim1Ptr = -1;
2397 *iPtr = 18;
2398 *ndigitsPtr = 0;
2399 break;
2400 case TCL_DD_E_FORMAT:
2401 if (*ndigitsPtr <= 0) {
2402 *ndigitsPtr = 1;
2403 }
2404 *iLimPtr = *iLim1Ptr = *iPtr = *ndigitsPtr;
2405 break;
2406 case TCL_DD_F_FORMAT:
2407 *iPtr = *ndigitsPtr + k + 1;
2408 *iLimPtr = *iPtr;
2409 *iLim1Ptr = *iPtr - 1;
2410 if (*iPtr <= 0) {
2411 *iPtr = 1;
2412 }
2413 break;
2414 default:
2415 *iPtr = -1;
2416 *iLimPtr = -1;
2417 *iLim1Ptr = -1;
2418 Tcl_Panic("impossible conversion type in TclDoubleDigits");
2419 }
2420 }
2421
2422 /*
2423 *-----------------------------------------------------------------------------
2424 *
2425 * BumpUp --
2426 *
2427 * Increases a string of digits ending in a series of nines to
2428 * designate the next higher number. xxxxb9999... -> xxxx(b+1)0000...
2429 *
2430 * Results:
2431 * Returns a pointer to the end of the adjusted string.
2432 *
2433 * Side effects:
2434 * In the case that the string consists solely of '999999', sets it
2435 * to "1" and moves the decimal point (*kPtr) one place to the right.
2436 *
2437 *-----------------------------------------------------------------------------
2438 */
2439
2440
2441 inline static char*
BumpUp(char * s,char * retval,int * kPtr)2442 BumpUp(char* s, /* Cursor pointing one past the end of the
2443 * string */
2444 char* retval, /* Start of the string of digits */
2445 int* kPtr) /* Position of the decimal point */
2446 {
2447 while (*--s == '9') {
2448 if (s == retval) {
2449 ++(*kPtr);
2450 *s = '1';
2451 return s+1;
2452 }
2453 }
2454 ++*s;
2455 ++s;
2456 return s;
2457 }
2458
2459 /*
2460 *-----------------------------------------------------------------------------
2461 *
2462 * AdjustRange --
2463 *
2464 * Rescales a 'double' in preparation for formatting it using the
2465 * 'quick' double-to-string method.
2466 *
2467 * Results:
2468 * Returns the precision that has been lost in the prescaling as
2469 * a count of units in the least significant place.
2470 *
2471 *-----------------------------------------------------------------------------
2472 */
2473
2474 inline static int
AdjustRange(double * dPtr,int k)2475 AdjustRange(double* dPtr, /* INOUT: Number to adjust */
2476 int k) /* IN: floor(log10(d)) */
2477 {
2478 int ieps; /* Number of roundoff errors that have
2479 * accumulated */
2480 double d = *dPtr; /* Number to adjust */
2481 double ds;
2482 int i, j, j1;
2483
2484 ieps = 2;
2485
2486 if (k > 0) {
2487 /*
2488 * The number must be reduced to bring it into range.
2489 */
2490 ds = tens[k & 0xf];
2491 j = k >> 4;
2492 if (j & BLETCH) {
2493 j &= (BLETCH-1);
2494 d /= bigtens[N_BIGTENS - 1];
2495 ieps++;
2496 }
2497 i = 0;
2498 for (; j != 0; j>>=1) {
2499 if (j & 1) {
2500 ds *= bigtens[i];
2501 ++ieps;
2502 }
2503 ++i;
2504 }
2505 d /= ds;
2506 } else if ((j1 = -k) != 0) {
2507 /*
2508 * The number must be increased to bring it into range
2509 */
2510 d *= tens[j1 & 0xf];
2511 i = 0;
2512 for (j = j1>>4; j; j>>=1) {
2513 if (j & 1) {
2514 ieps++;
2515 d *= bigtens[i];
2516 }
2517 ++i;
2518 }
2519 }
2520
2521 *dPtr = d;
2522 return ieps;
2523 }
2524
2525 /*
2526 *-----------------------------------------------------------------------------
2527 *
2528 * ShorteningQuickFormat --
2529 *
2530 * Returns a 'quick' format of a double precision number to a string
2531 * of digits, preferring a shorter string of digits if the shorter
2532 * string is still within 1/2 ulp of the number.
2533 *
2534 * Results:
2535 * Returns the string of digits. Returns NULL if the 'quick' method
2536 * fails and the bignum method must be used.
2537 *
2538 * Side effects:
2539 * Stores the position of the decimal point at '*kPtr'.
2540 *
2541 *-----------------------------------------------------------------------------
2542 */
2543
2544 inline static char*
ShorteningQuickFormat(double d,int k,int ilim,double eps,char * retval,int * kPtr)2545 ShorteningQuickFormat(double d, /* Number to convert */
2546 int k, /* floor(log10(d)) */
2547 int ilim, /* Number of significant digits to return */
2548 double eps,
2549 /* Estimated roundoff error */
2550 char* retval,
2551 /* Buffer to receive the digit string */
2552 int* kPtr)
2553 /* Pointer to stash the position of
2554 * the decimal point */
2555 {
2556 char* s = retval; /* Cursor in the return value */
2557 int digit; /* Current digit */
2558 int i;
2559
2560 eps = 0.5 / tens[ilim-1] - eps;
2561 i = 0;
2562 for (;;) {
2563 /* Convert a digit */
2564
2565 digit = (int) d;
2566 d -= digit;
2567 *s++ = '0' + digit;
2568
2569 /*
2570 * Truncate the conversion if the string of digits is within
2571 * 1/2 ulp of the actual value.
2572 */
2573
2574 if (d < eps) {
2575 *kPtr = k;
2576 return s;
2577 }
2578 if ((1. - d) < eps) {
2579 *kPtr = k;
2580 return BumpUp(s, retval, kPtr);
2581 }
2582
2583 /*
2584 * Bail out if the conversion fails to converge to a sufficiently
2585 * precise value
2586 */
2587
2588 if (++i >= ilim) {
2589 return NULL;
2590 }
2591
2592 /*
2593 * Bring the next digit to the integer part.
2594 */
2595
2596 eps *= 10;
2597 d *= 10.0;
2598 }
2599 }
2600
2601 /*
2602 *-----------------------------------------------------------------------------
2603 *
2604 * StrictQuickFormat --
2605 *
2606 * Convert a double precision number of a string of a precise number
2607 * of digits, using the 'quick' double precision method.
2608 *
2609 * Results:
2610 * Returns the digit string, or NULL if the bignum method must be
2611 * used to do the formatting.
2612 *
2613 * Side effects:
2614 * Stores the position of the decimal point in '*kPtr'.
2615 *
2616 *-----------------------------------------------------------------------------
2617 */
2618
2619 inline static char*
StrictQuickFormat(double d,int k,int ilim,double eps,char * retval,int * kPtr)2620 StrictQuickFormat(double d, /* Number to convert */
2621 int k, /* floor(log10(d)) */
2622 int ilim, /* Number of significant digits to return */
2623 double eps, /* Estimated roundoff error */
2624 char* retval, /* Start of the digit string */
2625 int* kPtr) /* Pointer to stash the position of
2626 * the decimal point */
2627 {
2628 char* s = retval; /* Cursor in the return value */
2629 int digit; /* Current digit of the answer */
2630 int i;
2631
2632 eps *= tens[ilim-1];
2633 i = 1;
2634 for (;;) {
2635 /* Extract a digit */
2636 digit = (int) d;
2637 d -= digit;
2638 if (d == 0.0) {
2639 ilim = i;
2640 }
2641 *s++ = '0' + digit;
2642
2643 /*
2644 * When the given digit count is reached, handle trailing strings
2645 * of 0 and 9.
2646 */
2647 if (i == ilim) {
2648 if (d > 0.5 + eps) {
2649 *kPtr = k;
2650 return BumpUp(s, retval, kPtr);
2651 } else if (d < 0.5 - eps) {
2652 while (*--s == '0') {
2653 /* do nothing */
2654 }
2655 s++;
2656 *kPtr = k;
2657 return s;
2658 } else {
2659 return NULL;
2660 }
2661 }
2662
2663 /* Advance to the next digit */
2664 ++i;
2665 d *= 10.0;
2666 }
2667 }
2668
2669 /*
2670 *-----------------------------------------------------------------------------
2671 *
2672 * QuickConversion --
2673 *
2674 * Converts a floating point number the 'quick' way, when only a limited
2675 * number of digits is required and floating point arithmetic can
2676 * therefore be used for the intermediate results.
2677 *
2678 * Results:
2679 * Returns the converted string, or NULL if the bignum method must
2680 * be used.
2681 *
2682 *-----------------------------------------------------------------------------
2683 */
2684
2685 inline static char*
QuickConversion(double e,int k,int k_check,int flags,int len,int ilim,int ilim1,int * decpt,char ** endPtr)2686 QuickConversion(double e, /* Number to format */
2687 int k, /* floor(log10(d)), approximately */
2688 int k_check, /* 0 if k is exact, 1 if it may be too high */
2689 int flags, /* Flags passed to dtoa:
2690 * TCL_DD_SHORTEN_FLAG */
2691 int len, /* Length of the return value */
2692 int ilim, /* Number of digits to store */
2693 int ilim1, /* Number of digits to store if we
2694 * musguessed k */
2695 int* decpt, /* OUTPUT: Location of the decimal point */
2696 char** endPtr) /* OUTPUT: Pointer to the terminal null byte */
2697 {
2698 int ieps; /* Number of 1-ulp roundoff errors that have
2699 * accumulated in the calculation*/
2700 Double eps; /* Estimated roundoff error */
2701 char* retval; /* Returned string */
2702 char* end; /* Pointer to the terminal null byte in the
2703 * returned string */
2704 volatile double d; /* Workaround for a bug in mingw gcc 3.4.5 */
2705
2706 /*
2707 * Bring d into the range [1 .. 10)
2708 */
2709 ieps = AdjustRange(&e, k);
2710 d = e;
2711
2712 /*
2713 * If the guessed value of k didn't get d into range, adjust it
2714 * by one. If that leaves us outside the range in which quick format
2715 * is accurate, bail out.
2716 */
2717 if (k_check && d < 1. && ilim > 0) {
2718 if (ilim1 < 0) {
2719 return NULL;
2720 }
2721 ilim = ilim1;
2722 --k;
2723 d *= 10.0;
2724 ++ieps;
2725 }
2726
2727 /*
2728 * Compute estimated roundoff error
2729 */
2730 eps.d = ieps * d + 7.;
2731 eps.w.word0 -= (FP_PRECISION-1) << EXP_SHIFT;
2732
2733 /*
2734 * Handle the peculiar case where the result has no significant
2735 * digits.
2736 */
2737 retval = ckalloc(len + 1);
2738 if (ilim == 0) {
2739 d -= 5.;
2740 if (d > eps.d) {
2741 *retval = '1';
2742 *decpt = k;
2743 return retval;
2744 } else if (d < -eps.d) {
2745 *decpt = k;
2746 return retval;
2747 } else {
2748 ckfree(retval);
2749 return NULL;
2750 }
2751 }
2752
2753 /* Format the digit string */
2754
2755 if (flags & TCL_DD_SHORTEN_FLAG) {
2756 end = ShorteningQuickFormat(d, k, ilim, eps.d, retval, decpt);
2757 } else {
2758 end = StrictQuickFormat(d, k, ilim, eps.d, retval, decpt);
2759 }
2760 if (end == NULL) {
2761 ckfree(retval);
2762 return NULL;
2763 }
2764 *end = '\0';
2765 if (endPtr != NULL) {
2766 *endPtr = end;
2767 }
2768 return retval;
2769 }
2770
2771 /*
2772 *-----------------------------------------------------------------------------
2773 *
2774 * CastOutPowersOf2 --
2775 *
2776 * Adjust the factors 'b2', 'm2', and 's2' to cast out common powers
2777 * of 2 from numerator and denominator in preparation for the 'bignum'
2778 * method of floating point conversion.
2779 *
2780 *-----------------------------------------------------------------------------
2781 */
2782
2783 inline static void
CastOutPowersOf2(int * b2,int * m2,int * s2)2784 CastOutPowersOf2(int* b2, /* Power of 2 to multiply the significand */
2785 int* m2, /* Power of 2 to multiply 1/2 ulp */
2786 int* s2) /* Power of 2 to multiply the common
2787 * denominator */
2788 {
2789 int i;
2790 if (*m2 > 0 && *s2 > 0) { /* Find the smallest power of 2 in the
2791 * numerator */
2792 if (*m2 < *s2) { /* Find the lowest common denominatorr */
2793 i = *m2;
2794 } else {
2795 i = *s2;
2796 }
2797 *b2 -= i; /* Reduce to lowest terms */
2798 *m2 -= i;
2799 *s2 -= i;
2800 }
2801 }
2802
2803 /*
2804 *-----------------------------------------------------------------------------
2805 *
2806 * ShorteningInt64Conversion --
2807 *
2808 * Converts a double-precision number to the shortest string of
2809 * digits that reconverts exactly to the given number, or to
2810 * 'ilim' digits if that will yield a shorter result. The numerator and
2811 * denominator in David Gay's conversion algorithm are known to fit
2812 * in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
2813 *
2814 * Results:
2815 * Returns the string of significant decimal digits, in newly
2816 * allocated memory
2817 *
2818 * Side effects:
2819 * Stores the location of the decimal point in '*decpt' and the
2820 * location of the terminal null byte in '*endPtr'.
2821 *
2822 *-----------------------------------------------------------------------------
2823 */
2824
2825 inline static char*
ShorteningInt64Conversion(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int b5,int m2plus,int m2minus,int m5,int s2,int s5,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)2826 ShorteningInt64Conversion(Double* dPtr,
2827 /* Original number to convert */
2828 int convType,
2829 /* Type of conversion (shortest, Steele,
2830 E format, F format) */
2831 Tcl_WideUInt bw,
2832 /* Integer significand */
2833 int b2, int b5,
2834 /* Scale factor for the significand
2835 * in the numerator */
2836 int m2plus, int m2minus, int m5,
2837 /* Scale factors for 1/2 ulp in
2838 * the numerator (will be different if
2839 * bw == 1 */
2840 int s2, int s5,
2841 /* Scale factors for the denominator */
2842 int k,
2843 /* Number of output digits before the decimal
2844 * point */
2845 int len,
2846 /* Number of digits to allocate */
2847 int ilim,
2848 /* Number of digits to convert if b >= s */
2849 int ilim1,
2850 /* Number of digits to convert if b < s */
2851 int* decpt,
2852 /* OUTPUT: Position of the decimal point */
2853 char** endPtr)
2854 /* OUTPUT: Position of the terminal '\0'
2855 * at the end of the returned string */
2856 {
2857
2858 char* retval = ckalloc(len + 1);
2859 /* Output buffer */
2860 Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
2861 /* Numerator of the fraction being converted */
2862 Tcl_WideUInt S = wuipow5[s5] << s2;
2863 /* Denominator of the fraction being
2864 * converted */
2865 Tcl_WideUInt mplus, mminus; /* Ranges for testing whether the result
2866 * is within roundoff of being exact */
2867 int digit; /* Current output digit */
2868 char* s = retval; /* Cursor in the output buffer */
2869 int i; /* Current position in the output buffer */
2870
2871 /* Adjust if the logarithm was guessed wrong */
2872
2873 if (b < S) {
2874 b = 10 * b;
2875 ++m2plus; ++m2minus; ++m5;
2876 ilim = ilim1;
2877 --k;
2878 }
2879
2880 /* Compute roundoff ranges */
2881
2882 mplus = wuipow5[m5] << m2plus;
2883 mminus = wuipow5[m5] << m2minus;
2884
2885 /* Loop through the digits */
2886
2887 i = 1;
2888 for (;;) {
2889 digit = (int)(b / S);
2890 if (digit > 10) {
2891 Tcl_Panic("wrong digit!");
2892 }
2893 b = b % S;
2894
2895 /*
2896 * Does the current digit put us on the low side of the exact value
2897 * but within within roundoff of being exact?
2898 */
2899 if (b < mplus
2900 || (b == mplus
2901 && convType != TCL_DD_STEELE0
2902 && (dPtr->w.word1 & 1) == 0)) {
2903 /*
2904 * Make sure we shouldn't be rounding *up* instead,
2905 * in case the next number above is closer
2906 */
2907 if (2 * b > S
2908 || (2 * b == S
2909 && (digit & 1) != 0)) {
2910 ++digit;
2911 if (digit == 10) {
2912 *s++ = '9';
2913 s = BumpUp(s, retval, &k);
2914 break;
2915 }
2916 }
2917
2918 /* Stash the current digit */
2919
2920 *s++ = '0' + digit;
2921 break;
2922 }
2923
2924 /*
2925 * Does one plus the current digit put us within roundoff of the
2926 * number?
2927 */
2928 if (b > S - mminus
2929 || (b == S - mminus
2930 && convType != TCL_DD_STEELE0
2931 && (dPtr->w.word1 & 1) == 0)) {
2932 if (digit == 9) {
2933 *s++ = '9';
2934 s = BumpUp(s, retval, &k);
2935 break;
2936 }
2937 ++digit;
2938 *s++ = '0' + digit;
2939 break;
2940 }
2941
2942 /*
2943 * Have we converted all the requested digits?
2944 */
2945 *s++ = '0' + digit;
2946 if (i == ilim) {
2947 if (2*b > S
2948 || (2*b == S && (digit & 1) != 0)) {
2949 s = BumpUp(s, retval, &k);
2950 }
2951 break;
2952 }
2953
2954 /* Advance to the next digit */
2955
2956 b = 10 * b;
2957 mplus = 10 * mplus;
2958 mminus = 10 * mminus;
2959 ++i;
2960 }
2961
2962 /*
2963 * Endgame - store the location of the decimal point and the end of the
2964 * string.
2965 */
2966 *s = '\0';
2967 *decpt = k;
2968 if (endPtr) {
2969 *endPtr = s;
2970 }
2971 return retval;
2972 }
2973
2974 /*
2975 *-----------------------------------------------------------------------------
2976 *
2977 * StrictInt64Conversion --
2978 *
2979 * Converts a double-precision number to a fixed-length string of
2980 * 'ilim' digits that reconverts exactly to the given number.
2981 * ('ilim' should be replaced with 'ilim1' in the case where
2982 * log10(d) has been overestimated). The numerator and
2983 * denominator in David Gay's conversion algorithm are known to fit
2984 * in Tcl_WideUInt, giving considerably faster arithmetic than mp_int's.
2985 *
2986 * Results:
2987 * Returns the string of significant decimal digits, in newly
2988 * allocated memory
2989 *
2990 * Side effects:
2991 * Stores the location of the decimal point in '*decpt' and the
2992 * location of the terminal null byte in '*endPtr'.
2993 *
2994 *-----------------------------------------------------------------------------
2995 */
2996
2997 inline static char*
StrictInt64Conversion(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int b5,int s2,int s5,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)2998 StrictInt64Conversion(Double* dPtr,
2999 /* Original number to convert */
3000 int convType,
3001 /* Type of conversion (shortest, Steele,
3002 E format, F format) */
3003 Tcl_WideUInt bw,
3004 /* Integer significand */
3005 int b2, int b5,
3006 /* Scale factor for the significand
3007 * in the numerator */
3008 int s2, int s5,
3009 /* Scale factors for the denominator */
3010 int k,
3011 /* Number of output digits before the decimal
3012 * point */
3013 int len,
3014 /* Number of digits to allocate */
3015 int ilim,
3016 /* Number of digits to convert if b >= s */
3017 int ilim1,
3018 /* Number of digits to convert if b < s */
3019 int* decpt,
3020 /* OUTPUT: Position of the decimal point */
3021 char** endPtr)
3022 /* OUTPUT: Position of the terminal '\0'
3023 * at the end of the returned string */
3024 {
3025
3026 char* retval = ckalloc(len + 1);
3027 /* Output buffer */
3028 Tcl_WideUInt b = (bw * wuipow5[b5]) << b2;
3029 /* Numerator of the fraction being converted */
3030 Tcl_WideUInt S = wuipow5[s5] << s2;
3031 /* Denominator of the fraction being
3032 * converted */
3033 int digit; /* Current output digit */
3034 char* s = retval; /* Cursor in the output buffer */
3035 int i; /* Current position in the output buffer */
3036
3037 /* Adjust if the logarithm was guessed wrong */
3038
3039 if (b < S) {
3040 b = 10 * b;
3041 ilim = ilim1;
3042 --k;
3043 }
3044
3045 /* Loop through the digits */
3046
3047 i = 1;
3048 for (;;) {
3049 digit = (int)(b / S);
3050 if (digit > 10) {
3051 Tcl_Panic("wrong digit!");
3052 }
3053 b = b % S;
3054
3055 /*
3056 * Have we converted all the requested digits?
3057 */
3058 *s++ = '0' + digit;
3059 if (i == ilim) {
3060 if (2*b > S
3061 || (2*b == S && (digit & 1) != 0)) {
3062 s = BumpUp(s, retval, &k);
3063 } else {
3064 while (*--s == '0') {
3065 /* do nothing */
3066 }
3067 ++s;
3068 }
3069 break;
3070 }
3071
3072 /* Advance to the next digit */
3073
3074 b = 10 * b;
3075 ++i;
3076 }
3077
3078 /*
3079 * Endgame - store the location of the decimal point and the end of the
3080 * string.
3081 */
3082 *s = '\0';
3083 *decpt = k;
3084 if (endPtr) {
3085 *endPtr = s;
3086 }
3087 return retval;
3088 }
3089
3090 /*
3091 *-----------------------------------------------------------------------------
3092 *
3093 * ShouldBankerRoundUpPowD --
3094 *
3095 * Test whether bankers' rounding should round a digit up. Assumption
3096 * is made that the denominator of the fraction being tested is
3097 * a power of 2**DIGIT_BIT.
3098 *
3099 * Results:
3100 * Returns 1 iff the fraction is more than 1/2, or if the fraction
3101 * is exactly 1/2 and the digit is odd.
3102 *
3103 *-----------------------------------------------------------------------------
3104 */
3105
3106 inline static int
ShouldBankerRoundUpPowD(mp_int * b,int sd,int isodd)3107 ShouldBankerRoundUpPowD(mp_int* b,
3108 /* Numerator of the fraction */
3109 int sd, /* Denominator is 2**(sd*DIGIT_BIT) */
3110 int isodd)
3111 /* 1 if the digit is odd, 0 if even */
3112 {
3113 int i;
3114 static const mp_digit topbit = (1<<(DIGIT_BIT-1));
3115 if (b->used < sd || (b->dp[sd-1] & topbit) == 0) {
3116 return 0;
3117 }
3118 if (b->dp[sd-1] != topbit) {
3119 return 1;
3120 }
3121 for (i = sd-2; i >= 0; --i) {
3122 if (b->dp[i] != 0) {
3123 return 1;
3124 }
3125 }
3126 return isodd;
3127 }
3128
3129 /*
3130 *-----------------------------------------------------------------------------
3131 *
3132 * ShouldBankerRoundUpToNextPowD --
3133 *
3134 * Tests whether bankers' rounding will round down in the
3135 * "denominator is a power of 2**MP_DIGIT" case.
3136 *
3137 * Results:
3138 * Returns 1 if the rounding will be performed - which increases the
3139 * digit by one - and 0 otherwise.
3140 *
3141 *-----------------------------------------------------------------------------
3142 */
3143
3144 inline static int
ShouldBankerRoundUpToNextPowD(mp_int * b,mp_int * m,int sd,int convType,int isodd,mp_int * temp)3145 ShouldBankerRoundUpToNextPowD(mp_int* b,
3146 /* Numerator of the fraction */
3147 mp_int* m,
3148 /* Numerator of the rounding tolerance */
3149 int sd,
3150 /* Common denominator is 2**(sd*DIGIT_BIT) */
3151 int convType,
3152 /* Conversion type: STEELE defeats
3153 * round-to-even (Not sure why one wants to
3154 * do this; I copied it from Gay) FIXME */
3155 int isodd,
3156 /* 1 if the integer significand is odd */
3157 mp_int* temp)
3158 /* Work area for the calculation */
3159 {
3160 int i;
3161
3162 /*
3163 * Compare B and S-m -- which is the same as comparing B+m and S --
3164 * which we do by computing b+m and doing a bitwhack compare against
3165 * 2**(DIGIT_BIT*sd)
3166 */
3167 mp_add(b, m, temp);
3168 if (temp->used <= sd) { /* too few digits to be > S */
3169 return 0;
3170 }
3171 if (temp->used > sd+1 || temp->dp[sd] > 1) {
3172 /* >= 2s */
3173 return 1;
3174 }
3175 for (i = sd-1; i >= 0; --i) {
3176 /* check for ==s */
3177 if (temp->dp[i] != 0) { /* > s */
3178 return 1;
3179 }
3180 }
3181 if (convType == TCL_DD_STEELE0) {
3182 /* biased rounding */
3183 return 0;
3184 }
3185 return isodd;
3186 }
3187
3188 /*
3189 *-----------------------------------------------------------------------------
3190 *
3191 * ShorteningBignumConversionPowD --
3192 *
3193 * Converts a double-precision number to the shortest string of
3194 * digits that reconverts exactly to the given number, or to
3195 * 'ilim' digits if that will yield a shorter result. The denominator
3196 * in David Gay's conversion algorithm is known to be a power of
3197 * 2**DIGIT_BIT, and hence the division in the main loop may be replaced
3198 * by a digit shift and mask.
3199 *
3200 * Results:
3201 * Returns the string of significant decimal digits, in newly
3202 * allocated memory
3203 *
3204 * Side effects:
3205 * Stores the location of the decimal point in '*decpt' and the
3206 * location of the terminal null byte in '*endPtr'.
3207 *
3208 *-----------------------------------------------------------------------------
3209 */
3210
3211 inline static char*
ShorteningBignumConversionPowD(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int b5,int m2plus,int m2minus,int m5,int sd,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)3212 ShorteningBignumConversionPowD(Double* dPtr,
3213 /* Original number to convert */
3214 int convType,
3215 /* Type of conversion (shortest, Steele,
3216 E format, F format) */
3217 Tcl_WideUInt bw,
3218 /* Integer significand */
3219 int b2, int b5,
3220 /* Scale factor for the significand
3221 * in the numerator */
3222 int m2plus, int m2minus, int m5,
3223 /* Scale factors for 1/2 ulp in
3224 * the numerator (will be different if
3225 * bw == 1 */
3226 int sd,
3227 /* Scale factor for the denominator */
3228 int k,
3229 /* Number of output digits before the decimal
3230 * point */
3231 int len,
3232 /* Number of digits to allocate */
3233 int ilim,
3234 /* Number of digits to convert if b >= s */
3235 int ilim1,
3236 /* Number of digits to convert if b < s */
3237 int* decpt,
3238 /* OUTPUT: Position of the decimal point */
3239 char** endPtr)
3240 /* OUTPUT: Position of the terminal '\0'
3241 * at the end of the returned string */
3242 {
3243
3244 char* retval = ckalloc(len + 1);
3245 /* Output buffer */
3246 mp_int b; /* Numerator of the fraction being converted */
3247 mp_int mplus, mminus; /* Bounds for roundoff */
3248 mp_digit digit; /* Current output digit */
3249 char* s = retval; /* Cursor in the output buffer */
3250 int i; /* Index in the output buffer */
3251 mp_int temp;
3252 int r1;
3253
3254 /*
3255 * b = bw * 2**b2 * 5**b5
3256 * mminus = 5**m5
3257 */
3258
3259 TclBNInitBignumFromWideUInt(&b, bw);
3260 mp_init_set_int(&mminus, 1);
3261 MulPow5(&b, b5, &b);
3262 mp_mul_2d(&b, b2, &b);
3263
3264 /* Adjust if the logarithm was guessed wrong */
3265
3266 if (b.used <= sd) {
3267 mp_mul_d(&b, 10, &b);
3268 ++m2plus; ++m2minus; ++m5;
3269 ilim = ilim1;
3270 --k;
3271 }
3272
3273 /*
3274 * mminus = 5**m5 * 2**m2minus
3275 * mplus = 5**m5 * 2**m2plus
3276 */
3277
3278 mp_mul_2d(&mminus, m2minus, &mminus);
3279 MulPow5(&mminus, m5, &mminus);
3280 if (m2plus > m2minus) {
3281 mp_init_copy(&mplus, &mminus);
3282 mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
3283 }
3284 mp_init(&temp);
3285
3286 /* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
3287 * by mp_digit extraction */
3288
3289 i = 0;
3290 for (;;) {
3291 if (b.used <= sd) {
3292 digit = 0;
3293 } else {
3294 digit = b.dp[sd];
3295 if (b.used > sd+1 || digit >= 10) {
3296 Tcl_Panic("wrong digit!");
3297 }
3298 --b.used; mp_clamp(&b);
3299 }
3300
3301 /*
3302 * Does the current digit put us on the low side of the exact value
3303 * but within within roundoff of being exact?
3304 */
3305
3306 r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
3307 if (r1 == MP_LT
3308 || (r1 == MP_EQ
3309 && convType != TCL_DD_STEELE0
3310 && (dPtr->w.word1 & 1) == 0)) {
3311 /*
3312 * Make sure we shouldn't be rounding *up* instead,
3313 * in case the next number above is closer
3314 */
3315 if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
3316 ++digit;
3317 if (digit == 10) {
3318 *s++ = '9';
3319 s = BumpUp(s, retval, &k);
3320 break;
3321 }
3322 }
3323
3324 /* Stash the last digit */
3325
3326 *s++ = '0' + digit;
3327 break;
3328 }
3329
3330 /*
3331 * Does one plus the current digit put us within roundoff of the
3332 * number?
3333 */
3334
3335 if (ShouldBankerRoundUpToNextPowD(&b, &mminus, sd,
3336 convType, dPtr->w.word1 & 1,
3337 &temp)) {
3338 if (digit == 9) {
3339 *s++ = '9';
3340 s = BumpUp(s, retval, &k);
3341 break;
3342 }
3343 ++digit;
3344 *s++ = '0' + digit;
3345 break;
3346 }
3347
3348 /*
3349 * Have we converted all the requested digits?
3350 */
3351 *s++ = '0' + digit;
3352 if (i == ilim) {
3353 if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
3354 s = BumpUp(s, retval, &k);
3355 }
3356 break;
3357 }
3358
3359 /* Advance to the next digit */
3360
3361 mp_mul_d(&b, 10, &b);
3362 mp_mul_d(&mminus, 10, &mminus);
3363 if (m2plus > m2minus) {
3364 mp_mul_2d(&mminus, m2plus-m2minus, &mplus);
3365 }
3366 ++i;
3367 }
3368
3369 /*
3370 * Endgame - store the location of the decimal point and the end of the
3371 * string.
3372 */
3373 if (m2plus > m2minus) {
3374 mp_clear(&mplus);
3375 }
3376 mp_clear_multi(&b, &mminus, &temp, NULL);
3377 *s = '\0';
3378 *decpt = k;
3379 if (endPtr) {
3380 *endPtr = s;
3381 }
3382 return retval;
3383 }
3384
3385 /*
3386 *-----------------------------------------------------------------------------
3387 *
3388 * StrictBignumConversionPowD --
3389 *
3390 * Converts a double-precision number to a fixed-lengt string of
3391 * 'ilim' digits (or 'ilim1' if log10(d) has been overestimated.)
3392 * The denominator in David Gay's conversion algorithm is known to
3393 * be a power of 2**DIGIT_BIT, and hence the division in the main
3394 * loop may be replaced by a digit shift and mask.
3395 *
3396 * Results:
3397 * Returns the string of significant decimal digits, in newly
3398 * allocated memory.
3399 *
3400 * Side effects:
3401 * Stores the location of the decimal point in '*decpt' and the
3402 * location of the terminal null byte in '*endPtr'.
3403 *
3404 *-----------------------------------------------------------------------------
3405 */
3406
3407 inline static char*
StrictBignumConversionPowD(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int b5,int sd,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)3408 StrictBignumConversionPowD(Double* dPtr,
3409 /* Original number to convert */
3410 int convType,
3411 /* Type of conversion (shortest, Steele,
3412 E format, F format) */
3413 Tcl_WideUInt bw,
3414 /* Integer significand */
3415 int b2, int b5,
3416 /* Scale factor for the significand
3417 * in the numerator */
3418 int sd,
3419 /* Scale factor for the denominator */
3420 int k,
3421 /* Number of output digits before the decimal
3422 * point */
3423 int len,
3424 /* Number of digits to allocate */
3425 int ilim,
3426 /* Number of digits to convert if b >= s */
3427 int ilim1,
3428 /* Number of digits to convert if b < s */
3429 int* decpt,
3430 /* OUTPUT: Position of the decimal point */
3431 char** endPtr)
3432 /* OUTPUT: Position of the terminal '\0'
3433 * at the end of the returned string */
3434 {
3435
3436 char* retval = ckalloc(len + 1);
3437 /* Output buffer */
3438 mp_int b; /* Numerator of the fraction being converted */
3439 mp_digit digit; /* Current output digit */
3440 char* s = retval; /* Cursor in the output buffer */
3441 int i; /* Index in the output buffer */
3442 mp_int temp;
3443
3444 /*
3445 * b = bw * 2**b2 * 5**b5
3446 */
3447
3448 TclBNInitBignumFromWideUInt(&b, bw);
3449 MulPow5(&b, b5, &b);
3450 mp_mul_2d(&b, b2, &b);
3451
3452 /* Adjust if the logarithm was guessed wrong */
3453
3454 if (b.used <= sd) {
3455 mp_mul_d(&b, 10, &b);
3456 ilim = ilim1;
3457 --k;
3458 }
3459 mp_init(&temp);
3460
3461 /*
3462 * Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
3463 * by mp_digit extraction
3464 */
3465
3466 i = 1;
3467 for (;;) {
3468 if (b.used <= sd) {
3469 digit = 0;
3470 } else {
3471 digit = b.dp[sd];
3472 if (b.used > sd+1 || digit >= 10) {
3473 Tcl_Panic("wrong digit!");
3474 }
3475 --b.used; mp_clamp(&b);
3476 }
3477
3478 /*
3479 * Have we converted all the requested digits?
3480 */
3481 *s++ = '0' + digit;
3482 if (i == ilim) {
3483 if (ShouldBankerRoundUpPowD(&b, sd, digit&1)) {
3484 s = BumpUp(s, retval, &k);
3485 } else {
3486 while (*--s == '0') {
3487 /* do nothing */
3488 }
3489 ++s;
3490 }
3491 break;
3492 }
3493
3494 /* Advance to the next digit */
3495
3496 mp_mul_d(&b, 10, &b);
3497 ++i;
3498 }
3499
3500 /*
3501 * Endgame - store the location of the decimal point and the end of the
3502 * string.
3503 */
3504 mp_clear_multi(&b, &temp, NULL);
3505 *s = '\0';
3506 *decpt = k;
3507 if (endPtr) {
3508 *endPtr = s;
3509 }
3510 return retval;
3511 }
3512
3513 /*
3514 *-----------------------------------------------------------------------------
3515 *
3516 * ShouldBankerRoundUp --
3517 *
3518 * Tests whether a digit should be rounded up or down when finishing
3519 * bignum-based floating point conversion.
3520 *
3521 * Results:
3522 * Returns 1 if the number needs to be rounded up, 0 otherwise.
3523 *
3524 *-----------------------------------------------------------------------------
3525 */
3526
3527 inline static int
ShouldBankerRoundUp(mp_int * twor,mp_int * S,int isodd)3528 ShouldBankerRoundUp(mp_int* twor,
3529 /* 2x the remainder from thd division that
3530 * produced the last digit */
3531 mp_int* S, /* Denominator */
3532 int isodd) /* Flag == 1 if the last digit is odd */
3533 {
3534 int r = mp_cmp_mag(twor, S);
3535 switch (r) {
3536 case MP_LT:
3537 return 0;
3538 case MP_EQ:
3539 return isodd;
3540 case MP_GT:
3541 return 1;
3542 }
3543 Tcl_Panic("in ShouldBankerRoundUp, trichotomy fails!");
3544 return 0;
3545 }
3546
3547 /*
3548 *-----------------------------------------------------------------------------
3549 *
3550 * ShouldBankerRoundUpToNext --
3551 *
3552 * Tests whether the remainder is great enough to force rounding
3553 * to the next higher digit.
3554 *
3555 * Results:
3556 * Returns 1 if the number should be rounded up, 0 otherwise.
3557 *
3558 *-----------------------------------------------------------------------------
3559 */
3560
3561 inline static int
ShouldBankerRoundUpToNext(mp_int * b,mp_int * m,mp_int * S,int convType,int isodd,mp_int * temp)3562 ShouldBankerRoundUpToNext(mp_int* b,
3563 /* Remainder from the division that produced
3564 * the last digit. */
3565 mp_int* m,
3566 /* Numerator of the rounding tolerance */
3567 mp_int* S,
3568 /* Denominator */
3569 int convType,
3570 /* Conversion type: STEELE0 defeats
3571 * round-to-even. (Not sure why one would
3572 * want this; I coped it from Gay. FIXME */
3573 int isodd,
3574 /* 1 if the integer significand is odd */
3575 mp_int* temp)
3576 /* Work area needed for the calculation */
3577 {
3578 int r;
3579 /* Compare b and S-m: this is the same as comparing B+m and S. */
3580 mp_add(b, m, temp);
3581 r = mp_cmp_mag(temp, S);
3582 switch(r) {
3583 case MP_LT:
3584 return 0;
3585 case MP_EQ:
3586 if (convType == TCL_DD_STEELE0) {
3587 return 0;
3588 } else {
3589 return isodd;
3590 }
3591 case MP_GT:
3592 return 1;
3593 }
3594 Tcl_Panic("in ShouldBankerRoundUpToNext, trichotomy fails!");
3595 return 0;
3596 }
3597
3598 /*
3599 *-----------------------------------------------------------------------------
3600 *
3601 * ShorteningBignumConversion --
3602 *
3603 * Convert a floating point number to a variable-length digit string
3604 * using the multiprecision method.
3605 *
3606 * Results:
3607 * Returns the string of digits.
3608 *
3609 * Side effects:
3610 * Stores the position of the decimal point in *decpt.
3611 * Stores a pointer to the end of the number in *endPtr.
3612 *
3613 *-----------------------------------------------------------------------------
3614 */
3615
3616 inline static char*
ShorteningBignumConversion(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int m2plus,int m2minus,int s2,int s5,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)3617 ShorteningBignumConversion(Double* dPtr,
3618 /* Original number being converted */
3619 int convType,
3620 /* Conversion type */
3621 Tcl_WideUInt bw,
3622 /* Integer significand and exponent */
3623 int b2,
3624 /* Scale factor for the significand */
3625 int m2plus, int m2minus,
3626 /* Scale factors for 1/2 ulp in numerator */
3627 int s2, int s5,
3628 /* Scale factors for denominator */
3629 int k,
3630 /* Guessed position of the decimal point */
3631 int len,
3632 /* Size of the digit buffer to allocate */
3633 int ilim,
3634 /* Number of digits to convert if b >= s */
3635 int ilim1,
3636 /* Number of digits to convert if b < s */
3637 int* decpt,
3638 /* OUTPUT: Position of the decimal point */
3639 char** endPtr)
3640 /* OUTPUT: Pointer to the end of the number */
3641 {
3642 char* retval = ckalloc(len+1);
3643 /* Buffer of digits to return */
3644 char* s = retval; /* Cursor in the return value */
3645 mp_int b; /* Numerator of the result */
3646 mp_int mminus; /* 1/2 ulp below the result */
3647 mp_int mplus; /* 1/2 ulp above the result */
3648 mp_int S; /* Denominator of the result */
3649 mp_int dig; /* Current digit of the result */
3650 int digit; /* Current digit of the result */
3651 mp_int temp; /* Work area */
3652 int minit = 1; /* Fudge factor for when we misguess k */
3653 int i;
3654 int r1;
3655
3656 /*
3657 * b = bw * 2**b2 * 5**b5
3658 * S = 2**s2 * 5*s5
3659 */
3660
3661 TclBNInitBignumFromWideUInt(&b, bw);
3662 mp_mul_2d(&b, b2, &b);
3663 mp_init_set_int(&S, 1);
3664 MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
3665
3666 /*
3667 * Handle the case where we guess the position of the decimal point
3668 * wrong.
3669 */
3670
3671 if (mp_cmp_mag(&b, &S) == MP_LT) {
3672 mp_mul_d(&b, 10, &b);
3673 minit = 10;
3674 ilim =ilim1;
3675 --k;
3676 }
3677
3678 /* mminus = 2**m2minus * 5**m5 */
3679
3680 mp_init_set_int(&mminus, minit);
3681 mp_mul_2d(&mminus, m2minus, &mminus);
3682 if (m2plus > m2minus) {
3683 mp_init_copy(&mplus, &mminus);
3684 mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
3685 }
3686 mp_init(&temp);
3687
3688 /* Loop through the digits */
3689
3690 mp_init(&dig);
3691 i = 1;
3692 for (;;) {
3693 mp_div(&b, &S, &dig, &b);
3694 if (dig.used > 1 || dig.dp[0] >= 10) {
3695 Tcl_Panic("wrong digit!");
3696 }
3697 digit = dig.dp[0];
3698
3699 /*
3700 * Does the current digit leave us with a remainder small enough to
3701 * round to it?
3702 */
3703
3704 r1 = mp_cmp_mag(&b, (m2plus > m2minus)? &mplus : &mminus);
3705 if (r1 == MP_LT
3706 || (r1 == MP_EQ
3707 && convType != TCL_DD_STEELE0
3708 && (dPtr->w.word1 & 1) == 0)) {
3709 mp_mul_2d(&b, 1, &b);
3710 if (ShouldBankerRoundUp(&b, &S, digit&1)) {
3711 ++digit;
3712 if (digit == 10) {
3713 *s++ = '9';
3714 s = BumpUp(s, retval, &k);
3715 break;
3716 }
3717 }
3718 *s++ = '0' + digit;
3719 break;
3720 }
3721
3722 /*
3723 * Does the current digit leave us with a remainder large enough
3724 * to commit to rounding up to the next higher digit?
3725 */
3726
3727 if (ShouldBankerRoundUpToNext(&b, &mminus, &S, convType,
3728 dPtr->w.word1 & 1, &temp)) {
3729 ++digit;
3730 if (digit == 10) {
3731 *s++ = '9';
3732 s = BumpUp(s, retval, &k);
3733 break;
3734 }
3735 *s++ = '0' + digit;
3736 break;
3737 }
3738
3739 /* Have we converted all the requested digits? */
3740
3741 *s++ = '0' + digit;
3742 if (i == ilim) {
3743 mp_mul_2d(&b, 1, &b);
3744 if (ShouldBankerRoundUp(&b, &S, digit&1)) {
3745 s = BumpUp(s, retval, &k);
3746 }
3747 break;
3748 }
3749
3750 /* Advance to the next digit */
3751
3752 if (s5 > 0) {
3753
3754 /* Can possibly shorten the denominator */
3755 mp_mul_2d(&b, 1, &b);
3756 mp_mul_2d(&mminus, 1, &mminus);
3757 if (m2plus > m2minus) {
3758 mp_mul_2d(&mplus, 1, &mplus);
3759 }
3760 mp_div_d(&S, 5, &S, NULL);
3761 --s5;
3762 /*
3763 * IDEA: It might possibly be a win to fall back to
3764 * int64 arithmetic here if S < 2**64/10. But it's
3765 * a win only for a fairly narrow range of magnitudes
3766 * so perhaps not worth bothering. We already know that
3767 * we shorten the denominator by at least 1 mp_digit, perhaps
3768 * 2. as we do the conversion for 17 digits of significance.
3769 * Possible savings:
3770 * 10**26 1 trip through loop before fallback possible
3771 * 10**27 1 trip
3772 * 10**28 2 trips
3773 * 10**29 3 trips
3774 * 10**30 4 trips
3775 * 10**31 5 trips
3776 * 10**32 6 trips
3777 * 10**33 7 trips
3778 * 10**34 8 trips
3779 * 10**35 9 trips
3780 * 10**36 10 trips
3781 * 10**37 11 trips
3782 * 10**38 12 trips
3783 * 10**39 13 trips
3784 * 10**40 14 trips
3785 * 10**41 15 trips
3786 * 10**42 16 trips
3787 * thereafter no gain.
3788 */
3789 } else {
3790 mp_mul_d(&b, 10, &b);
3791 mp_mul_d(&mminus, 10, &mminus);
3792 if (m2plus > m2minus) {
3793 mp_mul_2d(&mplus, 10, &mplus);
3794 }
3795 }
3796
3797 ++i;
3798 }
3799
3800
3801 /*
3802 * Endgame - store the location of the decimal point and the end of the
3803 * string.
3804 */
3805 if (m2plus > m2minus) {
3806 mp_clear(&mplus);
3807 }
3808 mp_clear_multi(&b, &mminus, &temp, &dig, &S, NULL);
3809 *s = '\0';
3810 *decpt = k;
3811 if (endPtr) {
3812 *endPtr = s;
3813 }
3814 return retval;
3815
3816 }
3817
3818 /*
3819 *-----------------------------------------------------------------------------
3820 *
3821 * StrictBignumConversion --
3822 *
3823 * Convert a floating point number to a fixed-length digit string
3824 * using the multiprecision method.
3825 *
3826 * Results:
3827 * Returns the string of digits.
3828 *
3829 * Side effects:
3830 * Stores the position of the decimal point in *decpt.
3831 * Stores a pointer to the end of the number in *endPtr.
3832 *
3833 *-----------------------------------------------------------------------------
3834 */
3835
3836 inline static char*
StrictBignumConversion(Double * dPtr,int convType,Tcl_WideUInt bw,int b2,int s2,int s5,int k,int len,int ilim,int ilim1,int * decpt,char ** endPtr)3837 StrictBignumConversion(Double* dPtr,
3838 /* Original number being converted */
3839 int convType,
3840 /* Conversion type */
3841 Tcl_WideUInt bw,
3842 /* Integer significand and exponent */
3843 int b2, /* Scale factor for the significand */
3844 int s2, int s5,
3845 /* Scale factors for denominator */
3846 int k, /* Guessed position of the decimal point */
3847 int len, /* Size of the digit buffer to allocate */
3848 int ilim,
3849 /* Number of digits to convert if b >= s */
3850 int ilim1,
3851 /* Number of digits to convert if b < s */
3852 int* decpt,
3853 /* OUTPUT: Position of the decimal point */
3854 char** endPtr)
3855 /* OUTPUT: Pointer to the end of the number */
3856 {
3857 char* retval = ckalloc(len+1);
3858 /* Buffer of digits to return */
3859 char* s = retval; /* Cursor in the return value */
3860 mp_int b; /* Numerator of the result */
3861 mp_int S; /* Denominator of the result */
3862 mp_int dig; /* Current digit of the result */
3863 int digit; /* Current digit of the result */
3864 mp_int temp; /* Work area */
3865 int g; /* Size of the current digit groun */
3866 int i, j;
3867
3868 /*
3869 * b = bw * 2**b2 * 5**b5
3870 * S = 2**s2 * 5*s5
3871 */
3872
3873 mp_init_multi(&temp, &dig, NULL);
3874 TclBNInitBignumFromWideUInt(&b, bw);
3875 mp_mul_2d(&b, b2, &b);
3876 mp_init_set_int(&S, 1);
3877 MulPow5(&S, s5, &S); mp_mul_2d(&S, s2, &S);
3878
3879 /*
3880 * Handle the case where we guess the position of the decimal point
3881 * wrong.
3882 */
3883
3884 if (mp_cmp_mag(&b, &S) == MP_LT) {
3885 mp_mul_d(&b, 10, &b);
3886 ilim =ilim1;
3887 --k;
3888 }
3889
3890 /* Convert the leading digit */
3891
3892 i = 0;
3893 mp_div(&b, &S, &dig, &b);
3894 if (dig.used > 1 || dig.dp[0] >= 10) {
3895 Tcl_Panic("wrong digit!");
3896 }
3897 digit = dig.dp[0];
3898
3899 /* Is a single digit all that was requested? */
3900
3901 *s++ = '0' + digit;
3902 if (++i >= ilim) {
3903 mp_mul_2d(&b, 1, &b);
3904 if (ShouldBankerRoundUp(&b, &S, digit&1)) {
3905 s = BumpUp(s, retval, &k);
3906 }
3907 } else {
3908
3909 for (;;) {
3910
3911 /* Shift by a group of digits. */
3912
3913 g = ilim - i;
3914 if (g > DIGIT_GROUP) {
3915 g = DIGIT_GROUP;
3916 }
3917 if (s5 >= g) {
3918 mp_div_d(&S, dpow5[g], &S, NULL);
3919 s5 -= g;
3920 } else if (s5 > 0) {
3921 mp_div_d(&S, dpow5[s5], &S, NULL);
3922 mp_mul_d(&b, dpow5[g - s5], &b);
3923 s5 = 0;
3924 } else {
3925 mp_mul_d(&b, dpow5[g], &b);
3926 }
3927 mp_mul_2d(&b, g, &b);
3928
3929 /*
3930 * As with the shortening bignum conversion, it's possible at
3931 * this point that we will have reduced the denominator to
3932 * less than 2**64/10, at which point it would be possible to
3933 * fall back to to int64 arithmetic. But the potential payoff
3934 * is tremendously less - unless we're working in F format -
3935 * because we know that three groups of digits will always
3936 * suffice for %#.17e, the longest format that doesn't introduce
3937 * empty precision.
3938 */
3939
3940 /* Extract the next group of digits */
3941
3942 mp_div(&b, &S, &dig, &b);
3943 if (dig.used > 1) {
3944 Tcl_Panic("wrong digit!");
3945 }
3946 digit = dig.dp[0];
3947 for (j = g-1; j >= 0; --j) {
3948 int t = itens[j];
3949 *s++ = digit / t + '0';
3950 digit %= t;
3951 }
3952 i += g;
3953
3954 /* Have we converted all the requested digits? */
3955
3956 if (i == ilim) {
3957 mp_mul_2d(&b, 1, &b);
3958 if (ShouldBankerRoundUp(&b, &S, digit&1)) {
3959 s = BumpUp(s, retval, &k);
3960 } else {
3961 while (*--s == '0') {
3962 /* do nothing */
3963 }
3964 ++s;
3965 }
3966 break;
3967 }
3968 }
3969 }
3970 /*
3971 * Endgame - store the location of the decimal point and the end of the
3972 * string.
3973 */
3974 mp_clear_multi(&b, &S, &temp, &dig, NULL);
3975 *s = '\0';
3976 *decpt = k;
3977 if (endPtr) {
3978 *endPtr = s;
3979 }
3980 return retval;
3981
3982 }
3983
3984 /*
3985 *-----------------------------------------------------------------------------
3986 *
3987 * TclDoubleDigits --
3988 *
3989 * Core of Tcl's conversion of double-precision floating point numbers
3990 * to decimal.
3991 *
3992 * Results:
3993 * Returns a newly-allocated string of digits.
3994 *
3995 * Side effects:
3996 * Sets *decpt to the index of the character in the string before the
3997 * place that the decimal point should go. If 'endPtr' is not NULL,
3998 * sets endPtr to point to the terminating '\0' byte of the string.
3999 * Sets *sign to 1 if a minus sign should be printed with the number,
4000 * or 0 if a plus sign (or no sign) should appear.
4001 *
4002 * This function is a service routine that produces the string of digits
4003 * for floating-point-to-decimal conversion. It can do a number of things
4004 * according to the 'flags' argument. Valid values for 'flags' include:
4005 * TCL_DD_SHORTEST - This is the default for floating point conversion
4006 * if ::tcl_precision is 0. It constructs the shortest string
4007 * of digits that will reconvert to the given number when scanned.
4008 * For floating point numbers that are exactly between two
4009 * decimal numbers, it resolves using the 'round to even' rule.
4010 * With this value, the 'ndigits' parameter is ignored.
4011 * TCL_DD_STEELE - This value is not recommended and may be removed
4012 * in the future. It follows the conversion algorithm outlined
4013 * in "How to Print Floating-Point Numbers Accurately" by
4014 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90,
4015 * pp. 112-126]. This rule has the effect of rendering 1e23
4016 * as 9.9999999999999999e22 - which is a 'better' approximation
4017 * in the sense that it will reconvert correctly even if
4018 * a subsequent input conversion is 'round up' or 'round down'
4019 * rather than 'round to nearest', but is surprising otherwise.
4020 * TCL_DD_E_FORMAT - This value is used to prepare numbers for %e
4021 * format conversion (or for default floating->string if
4022 * tcl_precision is not 0). It constructs a string of at most
4023 * 'ndigits' digits, choosing the one that is closest to the
4024 * given number (and resolving ties with 'round to even').
4025 * It is allowed to return fewer than 'ndigits' if the number
4026 * converts exactly; if the TCL_DD_E_FORMAT|TCL_DD_SHORTEN_FLAG
4027 * is supplied instead, it also returns fewer digits if the
4028 * shorter string will still reconvert to the given input number.
4029 * In any case, strings of trailing zeroes are suppressed.
4030 * TCL_DD_F_FORMAT - This value is used to prepare numbers for %f
4031 * format conversion. It requests that conversion proceed until
4032 * 'ndigits' digits after the decimal point have been converted.
4033 * It is possible for this format to result in a zero-length
4034 * string if the number is sufficiently small. Again, it
4035 * is permissible for TCL_DD_F_FORMAT to return fewer digits
4036 * for a number that converts exactly, and changing the
4037 * argument to TCL_DD_F_FORMAT|TCL_DD_SHORTEN_FLAG will allow
4038 * the routine also to return fewer digits if the shorter string
4039 * will still reconvert without loss to the given input number.
4040 * Strings of trailing zeroes are suppressed.
4041 *
4042 * To any of these flags may be OR'ed TCL_DD_NO_QUICK; this flag
4043 * requires all calculations to be done in exact arithmetic. Normally,
4044 * E and F format with fewer than about 14 digits will be done with
4045 * a quick floating point approximation and fall back on the exact
4046 * arithmetic only if the input number is close enough to the
4047 * midpoint between two decimal strings that more precision is needed
4048 * to resolve which string is correct.
4049 *
4050 * The value stored in the 'decpt' argument on return may be negative
4051 * (indicating that the decimal point falls to the left of the string)
4052 * or greater than the length of the string. In addition, the value -9999
4053 * is used as a sentinel to indicate that the string is one of the special
4054 * values "Infinity" and "NaN", and that no decimal point should be inserted.
4055 *
4056 *-----------------------------------------------------------------------------
4057 */
4058 char*
TclDoubleDigits(double dv,int ndigits,int flags,int * decpt,int * sign,char ** endPtr)4059 TclDoubleDigits(double dv, /* Number to convert */
4060 int ndigits, /* Number of digits requested */
4061 int flags, /* Conversion flags */
4062 int* decpt, /* OUTPUT: Position of the decimal point */
4063 int* sign, /* OUTPUT: 1 if the result is negative */
4064 char** endPtr) /* OUTPUT: If not NULL, receives a pointer
4065 * to one character beyond the end
4066 * of the returned string */
4067 {
4068 int convType = (flags & TCL_DD_CONVERSION_TYPE_MASK);
4069 /* Type of conversion being performed
4070 * TCL_DD_SHORTEST0
4071 * TCL_DD_STEELE0
4072 * TCL_DD_E_FORMAT
4073 * TCL_DD_F_FORMAT */
4074 Double d; /* Union for deconstructing doubles */
4075 Tcl_WideUInt bw; /* Integer significand */
4076 int be; /* Power of 2 by which b must be multiplied */
4077 int bbits; /* Number of bits needed to represent b */
4078 int denorm; /* Flag == 1 iff the input number was
4079 * denormalized */
4080 int k; /* Estimate of floor(log10(d)) */
4081 int k_check; /* Flag == 1 if d is near enough to a
4082 * power of ten that k must be checked */
4083 int b2, b5, s2, s5; /* Powers of 2 and 5 in the numerator and
4084 * denominator of intermediate results */
4085 int ilim = -1, ilim1 = -1; /* Number of digits to convert, and number
4086 * to convert if log10(d) has been
4087 * overestimated */
4088 char* retval; /* Return value from this function */
4089 int i = -1;
4090
4091 /* Put the input number into a union for bit-whacking */
4092
4093 d.d = dv;
4094
4095 /*
4096 * Handle the cases of negative numbers (by taking the absolute value:
4097 * this includes -Inf and -NaN!), infinity, Not a Number, and zero.
4098 */
4099
4100 TakeAbsoluteValue(&d, sign);
4101 if ((d.w.word0 & EXP_MASK) == EXP_MASK) {
4102 return FormatInfAndNaN(&d, decpt, endPtr);
4103 }
4104 if (d.d == 0.0) {
4105 return FormatZero(decpt, endPtr);
4106 }
4107
4108 /*
4109 * Unpack the floating point into a wide integer and an exponent.
4110 * Determine the number of bits that the big integer requires, and
4111 * compute a quick approximation (which may be one too high) of
4112 * ceil(log10(d.d)).
4113 */
4114 denorm = ((d.w.word0 & EXP_MASK) == 0);
4115 DoubleToExpAndSig(d.d, &bw, &be, &bbits);
4116 k = ApproximateLog10(bw, be, bbits);
4117 k = BetterLog10(d.d, k, &k_check);
4118
4119 /* At this point, we have:
4120 * d is the number to convert.
4121 * bw are significand and exponent: d == bw*2**be,
4122 * bbits is the length of bw: 2**bbits-1 <= bw < 2**bbits
4123 * k is either ceil(log10(d)) or ceil(log10(d))+1. k_check is 0
4124 * if we know that k is exactly ceil(log10(d)) and 1 if we need to
4125 * check.
4126 * We want a rational number
4127 * r = b * 10**(1-k) = bw * 2**b2 * 5**b5 / (2**s2 / 5**s5),
4128 * with b2, b5, s2, s5 >= 0. Note that the most significant decimal
4129 * digit is floor(r) and that successive digits can be obtained
4130 * by setting r <- 10*floor(r) (or b <= 10 * (b % S)).
4131 * Find appropriate b2, b5, s2, s5.
4132 */
4133
4134 ComputeScale(be, k, &b2, &b5, &s2, &s5);
4135
4136 /*
4137 * Correct an incorrect caller-supplied 'ndigits'.
4138 * Also determine:
4139 * i = The maximum number of decimal digits that will be returned in the
4140 * formatted string. This is k + 1 + ndigits for F format, 18 for
4141 * shortest and Steele, and ndigits for E format.
4142 * ilim = The number of significant digits to convert if
4143 * k has been guessed correctly. This is -1 for shortest and Steele
4144 * (which stop when all significance has been lost), 'ndigits'
4145 * for E format, and 'k + 1 + ndigits' for F format.
4146 * ilim1 = The minimum number of significant digits to convert if
4147 * k has been guessed 1 too high. This, too, is -1 for shortest
4148 * and Steele, and 'ndigits' for E format, but it's 'ndigits-1'
4149 * for F format.
4150 */
4151
4152 SetPrecisionLimits(convType, k, &ndigits, &i, &ilim, &ilim1);
4153
4154 /*
4155 * Try to do low-precision conversion in floating point rather
4156 * than resorting to expensive multiprecision arithmetic
4157 */
4158 if (ilim >= 0 && ilim <= QUICK_MAX && !(flags & TCL_DD_NO_QUICK)) {
4159 if ((retval = QuickConversion(d.d, k, k_check, flags,
4160 i, ilim, ilim1,
4161 decpt, endPtr)) != NULL) {
4162 return retval;
4163 }
4164 }
4165
4166 /*
4167 * For shortening conversions, determine the upper and lower bounds
4168 * for the remainder at which we can stop.
4169 * m+ = (2**m2plus * 5**m5) / (2**s2 * 5**s5) is the limit on the
4170 * high side, and
4171 * m- = (2**m2minus * 5**m5) / (2**s2 * 5**s5) is the limit on the
4172 * low side.
4173 * We may need to increase s2 to put m2plus, m2minus, b2 over a
4174 * common denominator.
4175 */
4176
4177 if (flags & TCL_DD_SHORTEN_FLAG) {
4178 int m2minus = b2;
4179 int m2plus;
4180 int m5 = b5;
4181 int len = i;
4182
4183 /*
4184 * Find the quantity i so that (2**i*5**b5)/(2**s2*5**s5)
4185 * is 1/2 unit in the least significant place of the floating
4186 * point number.
4187 */
4188 if (denorm) {
4189 i = be + EXPONENT_BIAS + (FP_PRECISION-1);
4190 } else {
4191 i = 1 + FP_PRECISION - bbits;
4192 }
4193 b2 += i;
4194 s2 += i;
4195
4196 /*
4197 * Reduce the fractions to lowest terms, since the above calculation
4198 * may have left excess powers of 2 in numerator and denominator
4199 */
4200 CastOutPowersOf2(&b2, &m2minus, &s2);
4201
4202 /*
4203 * In the special case where bw==1, the nearest floating point number
4204 * to it on the low side is 1/4 ulp below it. Adjust accordingly.
4205 */
4206 m2plus = m2minus;
4207 if (!denorm && bw == 1) {
4208 ++b2;
4209 ++s2;
4210 ++m2plus;
4211 }
4212
4213 if (s5+1 < N_LOG2POW5
4214 && s2+1 + log2pow5[s5+1] <= 64) {
4215 /*
4216 * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
4217 * word, then all our intermediate calculations can be done
4218 * using exact 64-bit arithmetic with no need for expensive
4219 * multiprecision operations. (This will be true for all numbers
4220 * in the range [1.0e-3 .. 1.0e+24]).
4221 */
4222
4223 return ShorteningInt64Conversion(&d, convType, bw, b2, b5,
4224 m2plus, m2minus, m5,
4225 s2, s5, k, len, ilim, ilim1,
4226 decpt, endPtr);
4227 } else if (s5 == 0) {
4228 /*
4229 * The denominator is a power of 2, so we can replace division
4230 * by digit shifts. First we round up s2 to a multiple of
4231 * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
4232 * into a version of the comparison that's specialized for
4233 * the 'power of mp_digit in the denominator' case.
4234 */
4235 if (s2 % DIGIT_BIT != 0) {
4236 int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
4237 b2 += delta;
4238 m2plus += delta;
4239 m2minus += delta;
4240 s2 += delta;
4241 }
4242 return ShorteningBignumConversionPowD(&d, convType, bw, b2, b5,
4243 m2plus, m2minus, m5,
4244 s2/DIGIT_BIT, k, len,
4245 ilim, ilim1, decpt, endPtr);
4246 } else {
4247
4248 /*
4249 * Alas, there's no helpful special case; use full-up
4250 * bignum arithmetic for the conversion
4251 */
4252
4253 return ShorteningBignumConversion(&d, convType, bw,
4254 b2, m2plus, m2minus,
4255 s2, s5, k, len,
4256 ilim, ilim1, decpt, endPtr);
4257
4258 }
4259
4260 } else {
4261
4262 /* Non-shortening conversion */
4263
4264 int len = i;
4265
4266 /* Reduce numerator and denominator to lowest terms */
4267
4268 if (b2 >= s2 && s2 > 0) {
4269 b2 -= s2; s2 = 0;
4270 } else if (s2 >= b2 && b2 > 0) {
4271 s2 -= b2; b2 = 0;
4272 }
4273
4274 if (s5+1 < N_LOG2POW5
4275 && s2+1 + log2pow5[s5+1] <= 64) {
4276 /*
4277 * If 10*2**s2*5**s5 == 2**(s2+1)+5**(s5+1) fits in a 64-bit
4278 * word, then all our intermediate calculations can be done
4279 * using exact 64-bit arithmetic with no need for expensive
4280 * multiprecision operations.
4281 */
4282
4283 return StrictInt64Conversion(&d, convType, bw, b2, b5,
4284 s2, s5, k, len, ilim, ilim1,
4285 decpt, endPtr);
4286
4287 } else if (s5 == 0) {
4288 /*
4289 * The denominator is a power of 2, so we can replace division
4290 * by digit shifts. First we round up s2 to a multiple of
4291 * DIGIT_BIT, and adjust m2 and b2 accordingly. Then we launch
4292 * into a version of the comparison that's specialized for
4293 * the 'power of mp_digit in the denominator' case.
4294 */
4295 if (s2 % DIGIT_BIT != 0) {
4296 int delta = DIGIT_BIT - (s2 % DIGIT_BIT);
4297 b2 += delta;
4298 s2 += delta;
4299 }
4300 return StrictBignumConversionPowD(&d, convType, bw, b2, b5,
4301 s2/DIGIT_BIT, k, len,
4302 ilim, ilim1, decpt, endPtr);
4303 } else {
4304 /*
4305 * There are no helpful special cases, but at least we know
4306 * in advance how many digits we will convert. We can run the
4307 * conversion in steps of DIGIT_GROUP digits, so as to
4308 * have many fewer mp_int divisions.
4309 */
4310 return StrictBignumConversion(&d, convType, bw, b2, s2, s5,
4311 k, len, ilim, ilim1, decpt, endPtr);
4312 }
4313 }
4314 }
4315
4316
4317 /*
4318 *----------------------------------------------------------------------
4319 *
4320 * TclInitDoubleConversion --
4321 *
4322 * Initializes constants that are needed for conversions to and from
4323 * 'double'
4324 *
4325 * Results:
4326 * None.
4327 *
4328 * Side effects:
4329 * The log base 2 of the floating point radix, the number of bits in a
4330 * double mantissa, and a table of the powers of five and ten are
4331 * computed and stored.
4332 *
4333 *----------------------------------------------------------------------
4334 */
4335
4336 void
TclInitDoubleConversion(void)4337 TclInitDoubleConversion(void)
4338 {
4339 int i;
4340 int x;
4341 Tcl_WideUInt u;
4342 double d;
4343
4344 #ifdef IEEE_FLOATING_POINT
4345 union {
4346 double dv;
4347 Tcl_WideUInt iv;
4348 } bitwhack;
4349 #endif
4350
4351 #if defined(__sgi) && defined(_COMPILER_VERSION)
4352 union fpc_csr mipsCR;
4353
4354 mipsCR.fc_word = get_fpc_csr();
4355 mipsCR.fc_struct.flush = 0;
4356 set_fpc_csr(mipsCR.fc_word);
4357 #endif
4358
4359 /*
4360 * Initialize table of powers of 10 expressed as wide integers.
4361 */
4362
4363 maxpow10_wide = (int)
4364 floor(sizeof(Tcl_WideUInt) * CHAR_BIT * log(2.) / log(10.));
4365 pow10_wide = (Tcl_WideUInt *)
4366 ckalloc((maxpow10_wide + 1) * sizeof(Tcl_WideUInt));
4367 u = 1;
4368 for (i = 0; i < maxpow10_wide; ++i) {
4369 pow10_wide[i] = u;
4370 u *= 10;
4371 }
4372 pow10_wide[i] = u;
4373
4374 /*
4375 * Determine how many bits of precision a double has, and how many
4376 * decimal digits that represents.
4377 */
4378
4379 if (frexp((double) FLT_RADIX, &log2FLT_RADIX) != 0.5) {
4380 Tcl_Panic("This code doesn't work on a decimal machine!");
4381 }
4382 log2FLT_RADIX--;
4383 mantBits = DBL_MANT_DIG * log2FLT_RADIX;
4384 d = 1.0;
4385
4386 /*
4387 * Initialize a table of powers of ten that can be exactly represented
4388 * in a double.
4389 */
4390
4391 x = (int) (DBL_MANT_DIG * log((double) FLT_RADIX) / log(5.0));
4392 if (x < MAXPOW) {
4393 mmaxpow = x;
4394 } else {
4395 mmaxpow = MAXPOW;
4396 }
4397 for (i=0 ; i<=mmaxpow ; ++i) {
4398 pow10vals[i] = d;
4399 d *= 10.0;
4400 }
4401
4402 /*
4403 * Initialize a table of large powers of five.
4404 */
4405
4406 for (i=0; i<9; ++i) {
4407 mp_init(pow5 + i);
4408 }
4409 mp_set(pow5, 5);
4410 for (i=0; i<8; ++i) {
4411 mp_sqr(pow5+i, pow5+i+1);
4412 }
4413 mp_init_set_int(pow5_13, 1220703125);
4414 for (i = 1; i < 5; ++i) {
4415 mp_init(pow5_13 + i);
4416 mp_sqr(pow5_13 + i - 1, pow5_13 + i);
4417 }
4418
4419 /*
4420 * Determine the number of decimal digits to the left and right of the
4421 * decimal point in the largest and smallest double, the smallest double
4422 * that differs from zero, and the number of mp_digits needed to represent
4423 * the significand of a double.
4424 */
4425
4426 maxDigits = (int) ((DBL_MAX_EXP * log((double) FLT_RADIX)
4427 + 0.5 * log(10.)) / log(10.));
4428 minDigits = (int) floor((DBL_MIN_EXP - DBL_MANT_DIG)
4429 * log((double) FLT_RADIX) / log(10.));
4430 log10_DIGIT_MAX = (int) floor(DIGIT_BIT * log(2.) / log(10.));
4431
4432 /*
4433 * Nokia 770's software-emulated floating point is "middle endian": the
4434 * bytes within a 32-bit word are little-endian (like the native
4435 * integers), but the two words of a 'double' are presented most
4436 * significant word first.
4437 */
4438
4439 #ifdef IEEE_FLOATING_POINT
4440 bitwhack.dv = 1.000000238418579;
4441 /* 3ff0 0000 4000 0000 */
4442 if ((bitwhack.iv >> 32) == 0x3ff00000) {
4443 n770_fp = 0;
4444 } else if ((bitwhack.iv & 0xffffffff) == 0x3ff00000) {
4445 n770_fp = 1;
4446 } else {
4447 Tcl_Panic("unknown floating point word order on this machine");
4448 }
4449 #endif
4450 }
4451
4452 /*
4453 *----------------------------------------------------------------------
4454 *
4455 * TclFinalizeDoubleConversion --
4456 *
4457 * Cleans up this file on exit.
4458 *
4459 * Results:
4460 * None
4461 *
4462 * Side effects:
4463 * Memory allocated by TclInitDoubleConversion is freed.
4464 *
4465 *----------------------------------------------------------------------
4466 */
4467
4468 void
TclFinalizeDoubleConversion(void)4469 TclFinalizeDoubleConversion(void)
4470 {
4471 int i;
4472
4473 ckfree((char *) pow10_wide);
4474 for (i=0; i<9; ++i) {
4475 mp_clear(pow5 + i);
4476 }
4477 }
4478
4479 /*
4480 *----------------------------------------------------------------------
4481 *
4482 * Tcl_InitBignumFromDouble --
4483 *
4484 * Extracts the integer part of a double and converts it to an arbitrary
4485 * precision integer.
4486 *
4487 * Results:
4488 * None.
4489 *
4490 * Side effects:
4491 * Initializes the bignum supplied, and stores the converted number in
4492 * it.
4493 *
4494 *----------------------------------------------------------------------
4495 */
4496
4497 int
Tcl_InitBignumFromDouble(Tcl_Interp * interp,double d,mp_int * b)4498 Tcl_InitBignumFromDouble(
4499 Tcl_Interp *interp, /* For error message */
4500 double d, /* Number to convert */
4501 mp_int *b) /* Place to store the result */
4502 {
4503 double fract;
4504 int expt;
4505
4506 /*
4507 * Infinite values can't convert to bignum.
4508 */
4509
4510 if (TclIsInfinite(d)) {
4511 if (interp != NULL) {
4512 const char *s = "integer value too large to represent";
4513
4514 Tcl_SetObjResult(interp, Tcl_NewStringObj(s, -1));
4515 Tcl_SetErrorCode(interp, "ARITH", "IOVERFLOW", s, NULL);
4516 }
4517 return TCL_ERROR;
4518 }
4519
4520 fract = frexp(d,&expt);
4521 if (expt <= 0) {
4522 mp_init(b);
4523 mp_zero(b);
4524 } else {
4525 Tcl_WideInt w = (Tcl_WideInt) ldexp(fract, mantBits);
4526 int shift = expt - mantBits;
4527
4528 TclBNInitBignumFromWideInt(b, w);
4529 if (shift < 0) {
4530 mp_div_2d(b, -shift, b, NULL);
4531 } else if (shift > 0) {
4532 mp_mul_2d(b, shift, b);
4533 }
4534 }
4535 return TCL_OK;
4536 }
4537
4538 /*
4539 *----------------------------------------------------------------------
4540 *
4541 * TclBignumToDouble --
4542 *
4543 * Convert an arbitrary-precision integer to a native floating point
4544 * number.
4545 *
4546 * Results:
4547 * Returns the converted number. Sets errno to ERANGE if the number is
4548 * too large to convert.
4549 *
4550 *----------------------------------------------------------------------
4551 */
4552
4553 double
TclBignumToDouble(mp_int * a)4554 TclBignumToDouble(
4555 mp_int *a) /* Integer to convert. */
4556 {
4557 mp_int b;
4558 int bits, shift, i, lsb;
4559 double r;
4560
4561
4562 /*
4563 * We need a 'mantBits'-bit significand. Determine what shift will
4564 * give us that.
4565 */
4566
4567 bits = mp_count_bits(a);
4568 if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
4569 errno = ERANGE;
4570 if (a->sign == MP_ZPOS) {
4571 return HUGE_VAL;
4572 } else {
4573 return -HUGE_VAL;
4574 }
4575 }
4576 shift = mantBits - bits;
4577
4578 /*
4579 * If shift > 0, shift the significand left by the requisite number of
4580 * bits. If shift == 0, the significand is already exactly 'mantBits'
4581 * in length. If shift < 0, we will need to shift the significand right
4582 * by the requisite number of bits, and round it. If the '1-shift'
4583 * least significant bits are 0, but the 'shift'th bit is nonzero,
4584 * then the significand lies exactly between two values and must be
4585 * 'rounded to even'.
4586 */
4587
4588 mp_init(&b);
4589 if (shift == 0) {
4590 mp_copy(a, &b);
4591 } else if (shift > 0) {
4592 mp_mul_2d(a, shift, &b);
4593 } else if (shift < 0) {
4594 lsb = mp_cnt_lsb(a);
4595 if (lsb == -1-shift) {
4596
4597 /*
4598 * Round to even
4599 */
4600
4601 mp_div_2d(a, -shift, &b, NULL);
4602 if (mp_isodd(&b)) {
4603 if (b.sign == MP_ZPOS) {
4604 mp_add_d(&b, 1, &b);
4605 } else {
4606 mp_sub_d(&b, 1, &b);
4607 }
4608 }
4609 } else {
4610
4611 /*
4612 * Ordinary rounding
4613 */
4614
4615 mp_div_2d(a, -1-shift, &b, NULL);
4616 if (b.sign == MP_ZPOS) {
4617 mp_add_d(&b, 1, &b);
4618 } else {
4619 mp_sub_d(&b, 1, &b);
4620 }
4621 mp_div_2d(&b, 1, &b, NULL);
4622 }
4623 }
4624
4625 /*
4626 * Accumulate the result, one mp_digit at a time.
4627 */
4628
4629 r = 0.0;
4630 for (i=b.used-1 ; i>=0 ; --i) {
4631 r = ldexp(r, DIGIT_BIT) + b.dp[i];
4632 }
4633 mp_clear(&b);
4634
4635 /*
4636 * Scale the result to the correct number of bits.
4637 */
4638
4639 r = ldexp(r, bits - mantBits);
4640
4641 /*
4642 * Return the result with the appropriate sign.
4643 */
4644
4645 if (a->sign == MP_ZPOS) {
4646 return r;
4647 } else {
4648 return -r;
4649 }
4650 }
4651
4652 /*
4653 *-----------------------------------------------------------------------------
4654 *
4655 * TclCeil --
4656 *
4657 * Computes the smallest floating point number that is at least the
4658 * mp_int argument.
4659 *
4660 * Results:
4661 * Returns the floating point number.
4662 *
4663 *-----------------------------------------------------------------------------
4664 */
4665
4666 double
TclCeil(mp_int * a)4667 TclCeil(
4668 mp_int *a) /* Integer to convert. */
4669 {
4670 double r = 0.0;
4671 mp_int b;
4672
4673 mp_init(&b);
4674 if (mp_cmp_d(a, 0) == MP_LT) {
4675 mp_neg(a, &b);
4676 r = -TclFloor(&b);
4677 } else {
4678 int bits = mp_count_bits(a);
4679
4680 if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
4681 r = HUGE_VAL;
4682 } else {
4683 int i, exact = 1, shift = mantBits - bits;
4684
4685 if (shift > 0) {
4686 mp_mul_2d(a, shift, &b);
4687 } else if (shift < 0) {
4688 mp_int d;
4689 mp_init(&d);
4690 mp_div_2d(a, -shift, &b, &d);
4691 exact = mp_iszero(&d);
4692 mp_clear(&d);
4693 } else {
4694 mp_copy(a, &b);
4695 }
4696 if (!exact) {
4697 mp_add_d(&b, 1, &b);
4698 }
4699 for (i=b.used-1 ; i>=0 ; --i) {
4700 r = ldexp(r, DIGIT_BIT) + b.dp[i];
4701 }
4702 r = ldexp(r, bits - mantBits);
4703 }
4704 }
4705 mp_clear(&b);
4706 return r;
4707 }
4708
4709 /*
4710 *-----------------------------------------------------------------------------
4711 *
4712 * TclFloor --
4713 *
4714 * Computes the largest floating point number less than or equal to
4715 * the mp_int argument.
4716 *
4717 * Results:
4718 * Returns the floating point value.
4719 *
4720 *-----------------------------------------------------------------------------
4721 */
4722
4723 double
TclFloor(mp_int * a)4724 TclFloor(
4725 mp_int *a) /* Integer to convert. */
4726 {
4727 double r = 0.0;
4728 mp_int b;
4729
4730 mp_init(&b);
4731 if (mp_cmp_d(a, 0) == MP_LT) {
4732 mp_neg(a, &b);
4733 r = -TclCeil(&b);
4734 } else {
4735 int bits = mp_count_bits(a);
4736
4737 if (bits > DBL_MAX_EXP*log2FLT_RADIX) {
4738 r = DBL_MAX;
4739 } else {
4740 int i, shift = mantBits - bits;
4741
4742 if (shift > 0) {
4743 mp_mul_2d(a, shift, &b);
4744 } else if (shift < 0) {
4745 mp_div_2d(a, -shift, &b, NULL);
4746 } else {
4747 mp_copy(a, &b);
4748 }
4749 for (i=b.used-1 ; i>=0 ; --i) {
4750 r = ldexp(r, DIGIT_BIT) + b.dp[i];
4751 }
4752 r = ldexp(r, bits - mantBits);
4753 }
4754 }
4755 mp_clear(&b);
4756 return r;
4757 }
4758
4759 /*
4760 *----------------------------------------------------------------------
4761 *
4762 * BignumToBiasedFrExp --
4763 *
4764 * Convert an arbitrary-precision integer to a native floating point
4765 * number in the range [0.5,1) times a power of two. NOTE: Intentionally
4766 * converts to a number that's a few ulp too small, so that
4767 * RefineApproximation will not overflow near the high end of the
4768 * machine's arithmetic range.
4769 *
4770 * Results:
4771 * Returns the converted number.
4772 *
4773 * Side effects:
4774 * Stores the exponent of two in 'machexp'.
4775 *
4776 *----------------------------------------------------------------------
4777 */
4778
4779 static double
BignumToBiasedFrExp(mp_int * a,int * machexp)4780 BignumToBiasedFrExp(
4781 mp_int *a, /* Integer to convert */
4782 int *machexp) /* Power of two */
4783 {
4784 mp_int b;
4785 int bits;
4786 int shift;
4787 int i;
4788 double r;
4789
4790 /*
4791 * Determine how many bits we need, and extract that many from the input.
4792 * Round to nearest unit in the last place.
4793 */
4794
4795 bits = mp_count_bits(a);
4796 shift = mantBits - 2 - bits;
4797 mp_init(&b);
4798 if (shift > 0) {
4799 mp_mul_2d(a, shift, &b);
4800 } else if (shift < 0) {
4801 mp_div_2d(a, -shift, &b, NULL);
4802 } else {
4803 mp_copy(a, &b);
4804 }
4805
4806 /*
4807 * Accumulate the result, one mp_digit at a time.
4808 */
4809
4810 r = 0.0;
4811 for (i=b.used-1; i>=0; --i) {
4812 r = ldexp(r, DIGIT_BIT) + b.dp[i];
4813 }
4814 mp_clear(&b);
4815
4816 /*
4817 * Return the result with the appropriate sign.
4818 */
4819
4820 *machexp = bits - mantBits + 2;
4821 return ((a->sign == MP_ZPOS) ? r : -r);
4822 }
4823
4824 /*
4825 *----------------------------------------------------------------------
4826 *
4827 * Pow10TimesFrExp --
4828 *
4829 * Multiply a power of ten by a number expressed as fraction and
4830 * exponent.
4831 *
4832 * Results:
4833 * Returns the significand of the result.
4834 *
4835 * Side effects:
4836 * Overwrites the 'machexp' parameter with the exponent of the result.
4837 *
4838 * Assumes that 'exponent' is such that 10**exponent would be a double, even
4839 * though 'fraction*10**(machexp+exponent)' might overflow.
4840 *
4841 *----------------------------------------------------------------------
4842 */
4843
4844 static double
Pow10TimesFrExp(int exponent,double fraction,int * machexp)4845 Pow10TimesFrExp(
4846 int exponent, /* Power of 10 to multiply by */
4847 double fraction, /* Significand of multiplicand */
4848 int *machexp) /* On input, exponent of multiplicand. On
4849 * output, exponent of result. */
4850 {
4851 int i, j;
4852 int expt = *machexp;
4853 double retval = fraction;
4854
4855 if (exponent > 0) {
4856 /*
4857 * Multiply by 10**exponent
4858 */
4859
4860 retval = frexp(retval * pow10vals[exponent&0xf], &j);
4861 expt += j;
4862 for (i=4; i<9; ++i) {
4863 if (exponent & (1<<i)) {
4864 retval = frexp(retval * pow_10_2_n[i], &j);
4865 expt += j;
4866 }
4867 }
4868 } else if (exponent < 0) {
4869 /*
4870 * Divide by 10**-exponent
4871 */
4872
4873 retval = frexp(retval / pow10vals[(-exponent) & 0xf], &j);
4874 expt += j;
4875 for (i=4; i<9; ++i) {
4876 if ((-exponent) & (1<<i)) {
4877 retval = frexp(retval / pow_10_2_n[i], &j);
4878 expt += j;
4879 }
4880 }
4881 }
4882
4883 *machexp = expt;
4884 return retval;
4885 }
4886
4887 /*
4888 *----------------------------------------------------------------------
4889 *
4890 * SafeLdExp --
4891 *
4892 * Do an 'ldexp' operation, but handle denormals gracefully.
4893 *
4894 * Results:
4895 * Returns the appropriately scaled value.
4896 *
4897 * On some platforms, 'ldexp' fails when presented with a number too
4898 * small to represent as a normalized double. This routine does 'ldexp'
4899 * in two steps for those numbers, to return correctly denormalized
4900 * values.
4901 *
4902 *----------------------------------------------------------------------
4903 */
4904
4905 static double
SafeLdExp(double fract,int expt)4906 SafeLdExp(
4907 double fract,
4908 int expt)
4909 {
4910 int minexpt = DBL_MIN_EXP * log2FLT_RADIX;
4911 volatile double a, b, retval;
4912
4913 if (expt < minexpt) {
4914 a = ldexp(fract, expt - mantBits - minexpt);
4915 b = ldexp(1.0, mantBits + minexpt);
4916 retval = a * b;
4917 } else {
4918 retval = ldexp(fract, expt);
4919 }
4920 return retval;
4921 }
4922
4923 /*
4924 *----------------------------------------------------------------------
4925 *
4926 * TclFormatNaN --
4927 *
4928 * Makes the string representation of a "Not a Number"
4929 *
4930 * Results:
4931 * None.
4932 *
4933 * Side effects:
4934 * Stores the string representation in the supplied buffer, which must be
4935 * at least TCL_DOUBLE_SPACE characters.
4936 *
4937 *----------------------------------------------------------------------
4938 */
4939
4940 void
TclFormatNaN(double value,char * buffer)4941 TclFormatNaN(
4942 double value, /* The Not-a-Number to format. */
4943 char *buffer) /* String representation. */
4944 {
4945 #ifndef IEEE_FLOATING_POINT
4946 strcpy(buffer, "NaN");
4947 return;
4948 #else
4949 union {
4950 double dv;
4951 Tcl_WideUInt iv;
4952 } bitwhack;
4953
4954 bitwhack.dv = value;
4955 if (n770_fp) {
4956 bitwhack.iv = Nokia770Twiddle(bitwhack.iv);
4957 }
4958 if (bitwhack.iv & ((Tcl_WideUInt) 1 << 63)) {
4959 bitwhack.iv &= ~ ((Tcl_WideUInt) 1 << 63);
4960 *buffer++ = '-';
4961 }
4962 *buffer++ = 'N';
4963 *buffer++ = 'a';
4964 *buffer++ = 'N';
4965 bitwhack.iv &= (((Tcl_WideUInt) 1) << 51) - 1;
4966 if (bitwhack.iv != 0) {
4967 sprintf(buffer, "(%" TCL_LL_MODIFIER "x)", bitwhack.iv);
4968 } else {
4969 *buffer = '\0';
4970 }
4971 #endif /* IEEE_FLOATING_POINT */
4972 }
4973
4974 /*
4975 *----------------------------------------------------------------------
4976 *
4977 * Nokia770Twiddle --
4978 *
4979 * Transpose the two words of a number for Nokia 770 floating
4980 * point handling.
4981 *
4982 *----------------------------------------------------------------------
4983 */
4984
4985 static Tcl_WideUInt
Nokia770Twiddle(Tcl_WideUInt w)4986 Nokia770Twiddle(
4987 Tcl_WideUInt w) /* Number to transpose */
4988 {
4989 return (((w >> 32) & 0xffffffff) | (w << 32));
4990 }
4991
4992 /*
4993 *----------------------------------------------------------------------
4994 *
4995 * TclNokia770Doubles --
4996 *
4997 * Transpose the two words of a number for Nokia 770 floating
4998 * point handling.
4999 *
5000 *----------------------------------------------------------------------
5001 */
5002
5003 int
TclNokia770Doubles(void)5004 TclNokia770Doubles(void)
5005 {
5006 return n770_fp;
5007 }
5008
5009 /*
5010 * Local Variables:
5011 * mode: c
5012 * c-basic-offset: 4
5013 * fill-column: 78
5014 * End:
5015 */
5016