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