src/libsphinxbase/util/dtoa.c

   1 /****************************************************************
   2  *
   3  * The author of this software is David M. Gay.
   4  *
   5  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
   6  *
   7  * Permission to use, copy, modify, and distribute this software for any
   8  * purpose without fee is hereby granted, provided that this entire notice
   9  * is included in all copies of any software which is or includes a copy
  10  * or modification of this software and in all copies of the supporting
  11  * documentation for such software.
  12  *
  13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
  15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17  *
  18  ***************************************************************/
  19
  20 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  21  * with " at " changed at "@" and " dot " changed to ".").      */
  22
  23 /* On a machine with IEEE extended-precision registers, it is
  24  * necessary to specify double-precision (53-bit) rounding precision
  25  * before invoking strtod or dtoa.  If the machine uses (the equivalent
  26  * of) Intel 80x87 arithmetic, the call
  27  *      _control87(PC_53, MCW_PC);
  28  * does this with many compilers.  Whether this or another call is
  29  * appropriate depends on the compiler; for this to work, it may be
  30  * necessary to #include "float.h" or another system-dependent header
  31  * file.
  32  */
  33
  34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  35  *
  36  * This strtod returns a nearest machine number to the input decimal
  37  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
  38  * broken by the IEEE round-even rule.  Otherwise ties are broken by
  39  * biased rounding (add half and chop).
  40  *
  41  * Inspired loosely by William D. Clinger's paper "How to Read Floating
  42  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  43  *
  44  * Modifications:
  45  *
  46  *      1. We only require IEEE, IBM, or VAX double-precision
  47  *              arithmetic (not IEEE double-extended).
  48  *      2. We get by with floating-point arithmetic in a case that
  49  *              Clinger missed -- when we're computing d * 10^n
  50  *              for a small integer d and the integer n is not too
  51  *              much larger than 22 (the maximum integer k for which
  52  *              we can represent 10^k exactly), we may be able to
  53  *              compute (d*10^k) * 10^(e-k) with just one roundoff.
  54  *      3. Rather than a bit-at-a-time adjustment of the binary
  55  *              result in the hard case, we use floating-point
  56  *              arithmetic to determine the adjustment to within
  57  *              one bit; only in really hard cases do we need to
  58  *              compute a second residual.
  59  *      4. Because of 3., we don't need a large table of powers of 10
  60  *              for ten-to-e (just some small tables, e.g. of 10^k
  61  *              for 0 <= k <= 22).
  62  */
  63
  64 /*
  65  * This file has been modified to remove dtoa() and all
  66  * non-reentrancy.  This makes it slower, but it also makes life a lot
  67  * easier on Windows and other platforms without static lock
  68  * initializers (grumble).
  69  */
  70
  71 /* Added by dhuggins@cs.cmu.edu to use autoconf results. */
  72 /* We do not care about the VAX. */
  73 #include "config.h"
  74 #ifdef WORDS_BIGENDIAN
  75 #define IEEE_MC68k
  76 #else
  77 #define IEEE_8087
  78 #endif
  79 #ifndef HAVE_LONG_LONG
  80 #define NO_LONG_LONG
  81 #endif
  82 #define Omit_Private_Memory
  83 #include "sphinxbase/ckd_alloc.h"
  84 #undef USE_LOCALE
  85
  86 /* Correct totally bogus typedefs in this code. */
  87 #include "sphinxbase/prim_type.h"
  88 #define Long int32   /* ZOMG */
  89 #define ULong uint32 /* WTF */
  90
  91 /*
  92  * #define IEEE_8087 for IEEE-arithmetic machines where the least
  93  *      significant byte has the lowest address.
  94  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  95  *      significant byte has the lowest address.
  96  * #define Long int on machines with 32-bit ints and 64-bit longs.
  97  * #define IBM for IBM mainframe-style floating-point arithmetic.
  98  * #define VAX for VAX-style floating-point arithmetic (D_floating).
  99  * #define No_leftright to omit left-right logic in fast floating-point
 100  *      computation of dtoa.
 101  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 102  *      and strtod and dtoa should round accordingly.
 103  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
 104  *      and Honor_FLT_ROUNDS is not #defined.
 105  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
 106  *      that use extended-precision instructions to compute rounded
 107  *      products and quotients) with IBM.
 108  * #define ROUND_BIASED for IEEE-format with biased rounding.
 109  * #define Inaccurate_Divide for IEEE-format with correctly rounded
 110  *      products but inaccurate quotients, e.g., for Intel i860.
 111  * #define NO_LONG_LONG on machines that do not have a "long long"
 112  *      integer type (of >= 64 bits).  On such machines, you can
 113  *      #define Just_16 to store 16 bits per 32-bit Long when doing
 114  *      high-precision integer arithmetic.  Whether this speeds things
 115  *      up or slows things down depends on the machine and the number
 116  *      being converted.  If long long is available and the name is
 117  *      something other than "long long", #define Llong to be the name,
 118  *      and if "unsigned Llong" does not work as an unsigned version of
 119  *      Llong, #define #ULLong to be the corresponding unsigned type.
 120  * #define KR_headers for old-style C function headers.
 121  * #define Bad_float_h if your system lacks a float.h or if it does not
 122  *      define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
 123  *      FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
 124  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
 125  *      if memory is available and otherwise does something you deem
 126  *      appropriate.  If MALLOC is undefined, malloc will be invoked
 127  *      directly -- and assumed always to succeed.
 128  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
 129  *      memory allocations from a private pool of memory when possible.
 130  *      When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
 131  *      unless #defined to be a different length.  This default length
 132  *      suffices to get rid of MALLOC calls except for unusual cases,
 133  *      such as decimal-to-binary conversion of a very long string of
 134  *      digits.  The longest string dtoa can return is about 751 bytes
 135  *      long.  For conversions by strtod of strings of 800 digits and
 136  *      all dtoa conversions in single-threaded executions with 8-byte
 137  *      pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
 138  *      pointers, PRIVATE_MEM >= 7112 appears adequate.
 139  * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
 140  *      #defined automatically on IEEE systems.  On such systems,
 141  *      when INFNAN_CHECK is #defined, strtod checks
 142  *      for Infinity and NaN (case insensitively).  On some systems
 143  *      (e.g., some HP systems), it may be necessary to #define NAN_WORD0
 144  *      appropriately -- to the most significant word of a quiet NaN.
 145  *      (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
 146  *      When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
 147  *      strtod also accepts (case insensitively) strings of the form
 148  *      NaN(x), where x is a string of hexadecimal digits and spaces;
 149  *      if there is only one string of hexadecimal digits, it is taken
 150  *      for the 52 fraction bits of the resulting NaN; if there are two
 151  *      or more strings of hex digits, the first is for the high 20 bits,
 152  *      the second and subsequent for the low 32 bits, with intervening
 153  *      white space ignored; but if this results in none of the 52
 154  *      fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
 155  *      and NAN_WORD1 are used instead.
 156  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
 157  *      multiple threads.  In this case, you must provide (or suitably
 158  *      #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
 159  *      by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
 160  *      in pow5mult, ensures lazy evaluation of only one copy of high
 161  *      powers of 5; omitting this lock would introduce a small
 162  *      probability of wasting memory, but would otherwise be harmless.)
 163  *      You must also invoke freedtoa(s) to free the value s returned by
 164  *      dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
 165  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
 166  *      avoids underflows on inputs whose result does not underflow.
 167  *      If you #define NO_IEEE_Scale on a machine that uses IEEE-format
 168  *      floating-point numbers and flushes underflows to zero rather
 169  *      than implementing gradual underflow, then you must also #define
 170  *      Sudden_Underflow.
 171  * #define YES_ALIAS to permit aliasing certain double values with
 172  *      arrays of ULongs.  This leads to slightly better code with
 173  *      some compilers and was always used prior to 19990916, but it
 174  *      is not strictly legal and can cause trouble with aggressively
 175  *      optimizing compilers (e.g., gcc 2.95.1 under -O2).
 176  * #define USE_LOCALE to use the current locale's decimal_point value.
 177  * #define SET_INEXACT if IEEE arithmetic is being used and extra
 178  *      computation should be done to set the inexact flag when the
 179  *      result is inexact and avoid setting inexact when the result
 180  *      is exact.  In this case, dtoa.c must be compiled in
 181  *      an environment, perhaps provided by #include "dtoa.c" in a
 182  *      suitable wrapper, that defines two functions,
 183  *              int get_inexact(void);
 184  *              void clear_inexact(void);
 185  *      such that get_inexact() returns a nonzero value if the
 186  *      inexact bit is already set, and clear_inexact() sets the
 187  *      inexact bit to 0.  When SET_INEXACT is #defined, strtod
 188  *      also does extra computations to set the underflow and overflow
 189  *      flags when appropriate (i.e., when the result is tiny and
 190  *      inexact or when it is a numeric value rounded to +-infinity).
 191  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
 192  *      the result overflows to +-Infinity or underflows to 0.
 193  */
 194
 195 #ifndef Long
 196 #define Long long
 197 #endif
 198 #ifndef ULong
 199 typedef unsigned Long ULong;
 200 #endif
 201
 202 #ifdef DEBUG
 203 #include "stdio.h"
 204 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
 205 #endif
 206
 207 #include "stdlib.h"
 208 #include "string.h"
 209
 210 #ifdef USE_LOCALE
 211 #include "locale.h"
 212 #endif
 213
 214 /* Private memory and other non-reentrant stuff removed. */
 215
 216 #undef IEEE_Arith
 217 #undef Avoid_Underflow
 218 #ifdef IEEE_MC68k
 219 #define IEEE_Arith
 220 #endif
 221 #ifdef IEEE_8087
 222 #define IEEE_Arith
 223 #endif
 224
 225 #ifdef IEEE_Arith
 226 #ifndef NO_INFNAN_CHECK
 227 #undef INFNAN_CHECK
 228 #define INFNAN_CHECK
 229 #endif
 230 #else
 231 #undef INFNAN_CHECK
 232 #endif
 233
 234 #include "errno.h"
 235
 236 #ifdef Bad_float_h
 237
 238 #ifdef IEEE_Arith
 239 #define DBL_DIG 15
 240 #define DBL_MAX_10_EXP 308
 241 #define DBL_MAX_EXP 1024
 242 #define FLT_RADIX 2
 243 #endif /*IEEE_Arith*/
 244
 245 #ifdef IBM
 246 #define DBL_DIG 16
 247 #define DBL_MAX_10_EXP 75
 248 #define DBL_MAX_EXP 63
 249 #define FLT_RADIX 16
 250 #define DBL_MAX 7.2370055773322621e+75
 251 #endif
 252
 253 #ifdef VAX
 254 #define DBL_DIG 16
 255 #define DBL_MAX_10_EXP 38
 256 #define DBL_MAX_EXP 127
 257 #define FLT_RADIX 2
 258 #define DBL_MAX 1.7014118346046923e+38
 259 #endif
 260
 261 #ifndef LONG_MAX
 262 #define LONG_MAX 2147483647
 263 #endif
 264
 265 #else /* ifndef Bad_float_h */
 266 #include "float.h"
 267 #endif /* Bad_float_h */
 268
 269 #ifndef __MATH_H__
 270 #include "math.h"
 271 #endif
 272
 273 #ifdef __cplusplus
 274 extern "C" {
 275 #endif
 276
 277 #ifndef CONST
 278 #ifdef KR_headers
 279 #define CONST /* blank */
 280 #else
 281 #define CONST const
 282 #endif
 283 #endif
 284
 285
 286 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
 287 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
 288 #endif
 289
 290 /** Union to extract the bytes of a double */
 291 typedef union { double d; ULong L[2]; } U;
 292
 293 #ifdef YES_ALIAS
 294 #define dval(x) x
 295 #ifdef IEEE_8087
 296 #define word0(x) ((ULong *)&x)[1]
 297 #define word1(x) ((ULong *)&x)[0]
 298 #else
 299 #define word0(x) ((ULong *)&x)[0]
 300 #define word1(x) ((ULong *)&x)[1]
 301 #endif
 302 #else
 303 #ifdef IEEE_8087
 304 #define word0(x) ((U*)&x)->L[1]
 305 #define word1(x) ((U*)&x)->L[0]
 306 #else
 307 #define word0(x) ((U*)&x)->L[0]
 308 #define word1(x) ((U*)&x)->L[1]
 309 #endif
 310 #define dval(x) ((U*)&x)->d
 311 #endif
 312
 313 /* The following definition of Storeinc is appropriate for MIPS processors.
 314  * An alternative that might be better on some machines is
 315  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
 316  */
 317 #if defined(IEEE_8087) + defined(VAX)
 318 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
 319 ((unsigned short *)a)[0] = (unsigned short)c, a++)
 320 #else
 321 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
 322 ((unsigned short *)a)[1] = (unsigned short)c, a++)
 323 #endif
 324
 325 /* #define P DBL_MANT_DIG */
 326 /* Ten_pmax = floor(P*log(2)/log(5)) */
 327 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
 328 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
 329 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
 330
 331 #ifdef IEEE_Arith
 332 #define Exp_shift  20
 333 #define Exp_shift1 20
 334 #define Exp_msk1    0x100000
 335 #define Exp_msk11   0x100000
 336 #define Exp_mask  0x7ff00000
 337 #define P 53
 338 #define Bias 1023
 339 #define Emin (-1022)
 340 #define Exp_1  0x3ff00000
 341 #define Exp_11 0x3ff00000
 342 #define Ebits 11
 343 #define Frac_mask  0xfffff
 344 #define Frac_mask1 0xfffff
 345 #define Ten_pmax 22
 346 #define Bletch 0x10
 347 #define Bndry_mask  0xfffff
 348 #define Bndry_mask1 0xfffff
 349 #define LSB 1
 350 #define Sign_bit 0x80000000
 351 #define Log2P 1
 352 #define Tiny0 0
 353 #define Tiny1 1
 354 #define Quick_max 14
 355 #define Int_max 14
 356 #ifndef NO_IEEE_Scale
 357 #define Avoid_Underflow
 358 #ifdef Flush_Denorm     /* debugging option */
 359 #undef Sudden_Underflow
 360 #endif
 361 #endif
 362
 363 #ifndef Flt_Rounds
 364 #ifdef FLT_ROUNDS
 365 #define Flt_Rounds FLT_ROUNDS
 366 #else
 367 #define Flt_Rounds 1
 368 #endif
 369 #endif /*Flt_Rounds*/
 370
 371 #ifdef Honor_FLT_ROUNDS
 372 #define Rounding rounding
 373 #undef Check_FLT_ROUNDS
 374 #define Check_FLT_ROUNDS
 375 #else
 376 #define Rounding Flt_Rounds
 377 #endif
 378
 379 #else /* ifndef IEEE_Arith */
 380 #undef Check_FLT_ROUNDS
 381 #undef Honor_FLT_ROUNDS
 382 #undef SET_INEXACT
 383 #undef  Sudden_Underflow
 384 #define Sudden_Underflow
 385 #ifdef IBM
 386 #undef Flt_Rounds
 387 #define Flt_Rounds 0
 388 #define Exp_shift  24
 389 #define Exp_shift1 24
 390 #define Exp_msk1   0x1000000
 391 #define Exp_msk11  0x1000000
 392 #define Exp_mask  0x7f000000
 393 #define P 14
 394 #define Bias 65
 395 #define Exp_1  0x41000000
 396 #define Exp_11 0x41000000
 397 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
 398 #define Frac_mask  0xffffff
 399 #define Frac_mask1 0xffffff
 400 #define Bletch 4
 401 #define Ten_pmax 22
 402 #define Bndry_mask  0xefffff
 403 #define Bndry_mask1 0xffffff
 404 #define LSB 1
 405 #define Sign_bit 0x80000000
 406 #define Log2P 4
 407 #define Tiny0 0x100000
 408 #define Tiny1 0
 409 #define Quick_max 14
 410 #define Int_max 15
 411 #else /* VAX */
 412 #undef Flt_Rounds
 413 #define Flt_Rounds 1
 414 #define Exp_shift  23
 415 #define Exp_shift1 7
 416 #define Exp_msk1    0x80
 417 #define Exp_msk11   0x800000
 418 #define Exp_mask  0x7f80
 419 #define P 56
 420 #define Bias 129
 421 #define Exp_1  0x40800000
 422 #define Exp_11 0x4080
 423 #define Ebits 8
 424 #define Frac_mask  0x7fffff
 425 #define Frac_mask1 0xffff007f
 426 #define Ten_pmax 24
 427 #define Bletch 2
 428 #define Bndry_mask  0xffff007f
 429 #define Bndry_mask1 0xffff007f
 430 #define LSB 0x10000
 431 #define Sign_bit 0x8000
 432 #define Log2P 1
 433 #define Tiny0 0x80
 434 #define Tiny1 0
 435 #define Quick_max 15
 436 #define Int_max 15
 437 #endif /* IBM, VAX */
 438 #endif /* IEEE_Arith */
 439
 440 #ifndef IEEE_Arith
 441 #define ROUND_BIASED
 442 #endif
 443
 444 #ifdef RND_PRODQUOT
 445 #define rounded_product(a,b) a = rnd_prod(a, b)
 446 #define rounded_quotient(a,b) a = rnd_quot(a, b)
 447 #ifdef KR_headers
 448 extern double rnd_prod(), rnd_quot();
 449 #else
 450 extern double rnd_prod(double, double), rnd_quot(double, double);
 451 #endif
 452 #else
 453 #define rounded_product(a,b) a *= b
 454 #define rounded_quotient(a,b) a /= b
 455 #endif
 456
 457 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
 458 #define Big1 0xffffffff
 459
 460 #ifndef Pack_32
 461 #define Pack_32
 462 #endif
 463
 464 #ifdef KR_headers
 465 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
 466 #else
 467 #define FFFFFFFF 0xffffffffUL
 468 #endif
 469
 470 #ifdef NO_LONG_LONG
 471 #undef ULLong
 472 #ifdef Just_16
 473 #undef Pack_32
 474 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
 475  * This makes some inner loops simpler and sometimes saves work
 476  * during multiplications, but it often seems to make things slightly
 477  * slower.  Hence the default is now to store 32 bits per Long.
 478  */
 479 #endif
 480 #else   /* long long available */
 481 #ifndef Llong
 482 #define Llong long long
 483 #endif
 484 #ifndef ULLong
 485 #define ULLong unsigned Llong
 486 #endif
 487 #endif /* NO_LONG_LONG */
 488
 489 #ifndef MULTIPLE_THREADS
 490 #define ACQUIRE_DTOA_LOCK(n)    /*nothing*/
 491 #define FREE_DTOA_LOCK(n)       /*nothing*/
 492 #endif
 493
 494 #define Kmax 15
 495
 496 #ifdef __cplusplus
 497 extern "C" double sb_strtod(const char *s00, char **se);
 498 #endif
 499
 500  struct
 501 Bigint {
 502         struct Bigint *next;
 503         int k, maxwds, sign, wds;
 504         ULong x[1];
 505         };
 506
 507  typedef struct Bigint Bigint;
 508
 509  static Bigint *
 510 Balloc
 511 #ifdef KR_headers
 512         (k) int k;
 513 #else
 514         (int k)
 515 #endif
 516 {
 517         int x;
 518         size_t len;
 519         Bigint *rv;
 520
 521         x = 1 << k;
 522         len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
 523                 /sizeof(double);
 524         rv = ckd_malloc(len*sizeof(double));
 525         rv->k = k;
 526         rv->maxwds = x;
 527         rv->sign = rv->wds = 0;
 528         return rv;
 529 }
 530
 531  static void
 532 Bfree
 533 #ifdef KR_headers
 534         (v) Bigint *v;
 535 #else
 536         (Bigint *v)
 537 #endif
 538 {
 539         ckd_free(v);
 540 }
 541
 542 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
 543 y->wds*sizeof(Long) + 2*sizeof(int))
 544
 545  static Bigint *
 546 multadd
 547 #ifdef KR_headers
 548         (b, m, a) Bigint *b; int m, a;
 549 #else
 550         (Bigint *b, int m, int a)       /* multiply by m and add a */
 551 #endif
 552 {
 553         int i, wds;
 554 #ifdef ULLong
 555         ULong *x;
 556         ULLong carry, y;
 557 #else
 558         ULong carry, *x, y;
 559 #ifdef Pack_32
 560         ULong xi, z;
 561 #endif
 562 #endif
 563         Bigint *b1;
 564
 565         wds = b->wds;
 566         x = b->x;
 567         i = 0;
 568         carry = a;
 569         do {
 570 #ifdef ULLong
 571                 y = *x * (ULLong)m + carry;
 572                 carry = y >> 32;
 573                 *x++ = y & FFFFFFFF;
 574 #else
 575 #ifdef Pack_32
 576                 xi = *x;
 577                 y = (xi & 0xffff) * m + carry;
 578                 z = (xi >> 16) * m + (y >> 16);
 579                 carry = z >> 16;
 580                 *x++ = (z << 16) + (y & 0xffff);
 581 #else
 582                 y = *x * m + carry;
 583                 carry = y >> 16;
 584                 *x++ = y & 0xffff;
 585 #endif
 586 #endif
 587                 }
 588                 while(++i < wds);
 589         if (carry) {
 590                 if (wds >= b->maxwds) {
 591                         b1 = Balloc(b->k+1);
 592                         Bcopy(b1, b);
 593                         Bfree(b);
 594                         b = b1;
 595                         }
 596                 b->x[wds++] = carry;
 597                 b->wds = wds;
 598                 }
 599         return b;
 600         }
 601
 602  static Bigint *
 603 s2b
 604 #ifdef KR_headers
 605         (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
 606 #else
 607         (CONST char *s, int nd0, int nd, ULong y9)
 608 #endif
 609 {
 610         Bigint *b;
 611         int i, k;
 612         Long x, y;
 613
 614         x = (nd + 8) / 9;
 615         for(k = 0, y = 1; x > y; y <<= 1, k++) ;
 616 #ifdef Pack_32
 617         b = Balloc(k);
 618         b->x[0] = y9;
 619         b->wds = 1;
 620 #else
 621         b = Balloc(k+1);
 622         b->x[0] = y9 & 0xffff;
 623         b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
 624 #endif
 625
 626         i = 9;
 627         if (9 < nd0) {
 628                 s += 9;
 629                 do b = multadd(b, 10, *s++ - '0');
 630                         while(++i < nd0);
 631                 s++;
 632                 }
 633         else
 634                 s += 10;
 635         for(; i < nd; i++)
 636                 b = multadd(b, 10, *s++ - '0');
 637         return b;
 638         }
 639
 640  static int
 641 hi0bits
 642 #ifdef KR_headers
 643         (x) register ULong x;
 644 #else
 645         (register ULong x)
 646 #endif
 647 {
 648         register int k = 0;
 649
 650         if (!(x & 0xffff0000)) {
 651                 k = 16;
 652                 x <<= 16;
 653                 }
 654         if (!(x & 0xff000000)) {
 655                 k += 8;
 656                 x <<= 8;
 657                 }
 658         if (!(x & 0xf0000000)) {
 659                 k += 4;
 660                 x <<= 4;
 661                 }
 662         if (!(x & 0xc0000000)) {
 663                 k += 2;
 664                 x <<= 2;
 665                 }
 666         if (!(x & 0x80000000)) {
 667                 k++;
 668                 if (!(x & 0x40000000))
 669                         return 32;
 670                 }
 671         return k;
 672         }
 673
 674  static int
 675 lo0bits
 676 #ifdef KR_headers
 677         (y) ULong *y;
 678 #else
 679         (ULong *y)
 680 #endif
 681 {
 682         register int k;
 683         register ULong x = *y;
 684
 685         if (x & 7) {
 686                 if (x & 1)
 687                         return 0;
 688                 if (x & 2) {
 689                         *y = x >> 1;
 690                         return 1;
 691                         }
 692                 *y = x >> 2;
 693                 return 2;
 694                 }
 695         k = 0;
 696         if (!(x & 0xffff)) {
 697                 k = 16;
 698                 x >>= 16;
 699                 }
 700         if (!(x & 0xff)) {
 701                 k += 8;
 702                 x >>= 8;
 703                 }
 704         if (!(x & 0xf)) {
 705                 k += 4;
 706                 x >>= 4;
 707                 }
 708         if (!(x & 0x3)) {
 709                 k += 2;
 710                 x >>= 2;
 711                 }
 712         if (!(x & 1)) {
 713                 k++;
 714                 x >>= 1;
 715                 if (!x)
 716                         return 32;
 717                 }
 718         *y = x;
 719         return k;
 720         }
 721
 722  static Bigint *
 723 i2b
 724 #ifdef KR_headers
 725         (i) int i;
 726 #else
 727         (int i)
 728 #endif
 729 {
 730         Bigint *b;
 731
 732         b = Balloc(1);
 733         b->x[0] = i;
 734         b->wds = 1;
 735         return b;
 736         }
 737
 738  static Bigint *
 739 mult
 740 #ifdef KR_headers
 741         (a, b) Bigint *a, *b;
 742 #else
 743         (Bigint *a, Bigint *b)
 744 #endif
 745 {
 746         Bigint *c;
 747         int k, wa, wb, wc;
 748         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 749         ULong y;
 750 #ifdef ULLong
 751         ULLong carry, z;
 752 #else
 753         ULong carry, z;
 754 #ifdef Pack_32
 755         ULong z2;
 756 #endif
 757 #endif
 758
 759         if (a->wds < b->wds) {
 760                 c = a;
 761                 a = b;
 762                 b = c;
 763                 }
 764         k = a->k;
 765         wa = a->wds;
 766         wb = b->wds;
 767         wc = wa + wb;
 768         if (wc > a->maxwds)
 769                 k++;
 770         c = Balloc(k);
 771         for(x = c->x, xa = x + wc; x < xa; x++)
 772                 *x = 0;
 773         xa = a->x;
 774         xae = xa + wa;
 775         xb = b->x;
 776         xbe = xb + wb;
 777         xc0 = c->x;
 778 #ifdef ULLong
 779         for(; xb < xbe; xc0++) {
 780                 if ((y = *xb++)) {
 781                         x = xa;
 782                         xc = xc0;
 783                         carry = 0;
 784                         do {
 785                                 z = *x++ * (ULLong)y + *xc + carry;
 786                                 carry = z >> 32;
 787                                 *xc++ = z & FFFFFFFF;
 788                                 }
 789                                 while(x < xae);
 790                         *xc = carry;
 791                         }
 792                 }
 793 #else
 794 #ifdef Pack_32
 795         for(; xb < xbe; xb++, xc0++) {
 796                 if (y = *xb & 0xffff) {
 797                         x = xa;
 798                         xc = xc0;
 799                         carry = 0;
 800                         do {
 801                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
 802                                 carry = z >> 16;
 803                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
 804                                 carry = z2 >> 16;
 805                                 Storeinc(xc, z2, z);
 806                                 }
 807                                 while(x < xae);
 808                         *xc = carry;
 809                         }
 810                 if (y = *xb >> 16) {
 811                         x = xa;
 812                         xc = xc0;
 813                         carry = 0;
 814                         z2 = *xc;
 815                         do {
 816                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
 817                                 carry = z >> 16;
 818                                 Storeinc(xc, z, z2);
 819                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
 820                                 carry = z2 >> 16;
 821                                 }
 822                                 while(x < xae);
 823                         *xc = z2;
 824                         }
 825                 }
 826 #else
 827         for(; xb < xbe; xc0++) {
 828                 if (y = *xb++) {
 829                         x = xa;
 830                         xc = xc0;
 831                         carry = 0;
 832                         do {
 833                                 z = *x++ * y + *xc + carry;
 834                                 carry = z >> 16;
 835                                 *xc++ = z & 0xffff;
 836                                 }
 837                                 while(x < xae);
 838                         *xc = carry;
 839                         }
 840                 }
 841 #endif
 842 #endif
 843         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 844         c->wds = wc;
 845         return c;
 846         }
 847
 848  static Bigint *
 849 pow5mult
 850 #ifdef KR_headers
 851         (b, k) Bigint *b; int k;
 852 #else
 853         (Bigint *b, int k)
 854 #endif
 855 {
 856         Bigint *b1, *p5, *p51;
 857         int i;
 858         static int CONST p05[3] = { 5, 25, 125 };
 859
 860         if ((i = k & 3))
 861                 b = multadd(b, p05[i-1], 0);
 862
 863         if (!(k >>= 2))
 864                 return b;
 865
 866         p5 = i2b(625);
 867         for(;;) {
 868                 if (k & 1) {
 869                         b1 = mult(b, p5);
 870                         Bfree(b);
 871                         b = b1;
 872                 }
 873                 if (!(k >>= 1))
 874                         break;
 875                 p51 = mult(p5,p5);
 876                 Bfree(p5);
 877                 p5 = p51;
 878         }
 879         Bfree(p5);
 880         return b;
 881 }
 882
 883  static Bigint *
 884 lshift
 885 #ifdef KR_headers
 886         (b, k) Bigint *b; int k;
 887 #else
 888         (Bigint *b, int k)
 889 #endif
 890 {
 891         int i, k1, n, n1;
 892         Bigint *b1;
 893         ULong *x, *x1, *xe, z;
 894
 895 #ifdef Pack_32
 896         n = k >> 5;
 897 #else
 898         n = k >> 4;
 899 #endif
 900         k1 = b->k;
 901         n1 = n + b->wds + 1;
 902         for(i = b->maxwds; n1 > i; i <<= 1)
 903                 k1++;
 904         b1 = Balloc(k1);
 905         x1 = b1->x;
 906         for(i = 0; i < n; i++)
 907                 *x1++ = 0;
 908         x = b->x;
 909         xe = x + b->wds;
 910 #ifdef Pack_32
 911         if (k &= 0x1f) {
 912                 k1 = 32 - k;
 913                 z = 0;
 914                 do {
 915                         *x1++ = *x << k | z;
 916                         z = *x++ >> k1;
 917                         }
 918                         while(x < xe);
 919                 if ((*x1 = z))
 920                         ++n1;
 921                 }
 922 #else
 923         if (k &= 0xf) {
 924                 k1 = 16 - k;
 925                 z = 0;
 926                 do {
 927                         *x1++ = *x << k  & 0xffff | z;
 928                         z = *x++ >> k1;
 929                         }
 930                         while(x < xe);
 931                 if (*x1 = z)
 932                         ++n1;
 933                 }
 934 #endif
 935         else do
 936                 *x1++ = *x++;
 937                 while(x < xe);
 938         b1->wds = n1 - 1;
 939         Bfree(b);
 940         return b1;
 941         }
 942
 943  static int
 944 cmp
 945 #ifdef KR_headers
 946         (a, b) Bigint *a, *b;
 947 #else
 948         (Bigint *a, Bigint *b)
 949 #endif
 950 {
 951         ULong *xa, *xa0, *xb, *xb0;
 952         int i, j;
 953
 954         i = a->wds;
 955         j = b->wds;
 956 #ifdef DEBUG
 957         if (i > 1 && !a->x[i-1])
 958                 Bug("cmp called with a->x[a->wds-1] == 0");
 959         if (j > 1 && !b->x[j-1])
 960                 Bug("cmp called with b->x[b->wds-1] == 0");
 961 #endif
 962         if (i -= j)
 963                 return i;
 964         xa0 = a->x;
 965         xa = xa0 + j;
 966         xb0 = b->x;
 967         xb = xb0 + j;
 968         for(;;) {
 969                 if (*--xa != *--xb)
 970                         return *xa < *xb ? -1 : 1;
 971                 if (xa <= xa0)
 972                         break;
 973                 }
 974         return 0;
 975         }
 976
 977  static Bigint *
 978 diff
 979 #ifdef KR_headers
 980         (a, b) Bigint *a, *b;
 981 #else
 982         (Bigint *a, Bigint *b)
 983 #endif
 984 {
 985         Bigint *c;
 986         int i, wa, wb;
 987         ULong *xa, *xae, *xb, *xbe, *xc;
 988 #ifdef ULLong
 989         ULLong borrow, y;
 990 #else
 991         ULong borrow, y;
 992 #ifdef Pack_32
 993         ULong z;
 994 #endif
 995 #endif
 996
 997         i = cmp(a,b);
 998         if (!i) {
 999                 c = Balloc(0);
1000                 c->wds = 1;
1001                 c->x[0] = 0;
1002                 return c;
1003                 }
1004         if (i < 0) {
1005                 c = a;
1006                 a = b;
1007                 b = c;
1008                 i = 1;
1009                 }
1010         else
1011                 i = 0;
1012         c = Balloc(a->k);
1013         c->sign = i;
1014         wa = a->wds;
1015         xa = a->x;
1016         xae = xa + wa;
1017         wb = b->wds;
1018         xb = b->x;
1019         xbe = xb + wb;
1020         xc = c->x;
1021         borrow = 0;
1022 #ifdef ULLong
1023         do {
1024                 y = (ULLong)*xa++ - *xb++ - borrow;
1025                 borrow = y >> 32 & (ULong)1;
1026                 *xc++ = y & FFFFFFFF;
1027                 }
1028                 while(xb < xbe);
1029         while(xa < xae) {
1030                 y = *xa++ - borrow;
1031                 borrow = y >> 32 & (ULong)1;
1032                 *xc++ = y & FFFFFFFF;
1033                 }
1034 #else
1035 #ifdef Pack_32
1036         do {
1037                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1038                 borrow = (y & 0x10000) >> 16;
1039                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1040                 borrow = (z & 0x10000) >> 16;
1041                 Storeinc(xc, z, y);
1042                 }
1043                 while(xb < xbe);
1044         while(xa < xae) {
1045                 y = (*xa & 0xffff) - borrow;
1046                 borrow = (y & 0x10000) >> 16;
1047                 z = (*xa++ >> 16) - borrow;
1048                 borrow = (z & 0x10000) >> 16;
1049                 Storeinc(xc, z, y);
1050                 }
1051 #else
1052         do {
1053                 y = *xa++ - *xb++ - borrow;
1054                 borrow = (y & 0x10000) >> 16;
1055                 *xc++ = y & 0xffff;
1056                 }
1057                 while(xb < xbe);
1058         while(xa < xae) {
1059                 y = *xa++ - borrow;
1060                 borrow = (y & 0x10000) >> 16;
1061                 *xc++ = y & 0xffff;
1062                 }
1063 #endif
1064 #endif
1065         while(!*--xc)
1066                 wa--;
1067         c->wds = wa;
1068         return c;
1069         }
1070
1071  static double
1072 ulp
1073 #ifdef KR_headers
1074         (x) double x;
1075 #else
1076         (double x)
1077 #endif
1078 {
1079         register Long L;
1080         U a;
1081
1082         L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1083 #ifndef Avoid_Underflow
1084 #ifndef Sudden_Underflow
1085         if (L > 0) {
1086 #endif
1087 #endif
1088 #ifdef IBM
1089                 L |= Exp_msk1 >> 4;
1090 #endif
1091                 word0(a) = L;
1092                 word1(a) = 0;
1093 #ifndef Avoid_Underflow
1094 #ifndef Sudden_Underflow
1095                 }
1096         else {
1097                 L = -L >> Exp_shift;
1098                 if (L < Exp_shift) {
1099                         word0(a) = 0x80000 >> L;
1100                         word1(a) = 0;
1101                         }
1102                 else {
1103                         word0(a) = 0;
1104                         L -= Exp_shift;
1105                         word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1106                         }
1107                 }
1108 #endif
1109 #endif
1110         return dval(a);
1111         }
1112
1113  static double
1114 b2d
1115 #ifdef KR_headers
1116         (a, e) Bigint *a; int *e;
1117 #else
1118         (Bigint *a, int *e)
1119 #endif
1120 {
1121         ULong *xa, *xa0, w, y, z;
1122         int k;
1123         U d;
1124 #ifdef VAX
1125         ULong d0, d1;
1126 #else
1127 #define d0 word0(d)
1128 #define d1 word1(d)
1129 #endif
1130
1131         xa0 = a->x;
1132         xa = xa0 + a->wds;
1133         y = *--xa;
1134 #ifdef DEBUG
1135         if (!y) Bug("zero y in b2d");
1136 #endif
1137         k = hi0bits(y);
1138         *e = 32 - k;
1139 #ifdef Pack_32
1140         if (k < Ebits) {
1141                 d0 = Exp_1 | y >> (Ebits - k);
1142                 w = xa > xa0 ? *--xa : 0;
1143                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1144                 goto ret_d;
1145                 }
1146         z = xa > xa0 ? *--xa : 0;
1147         if (k -= Ebits) {
1148                 d0 = Exp_1 | y << k | z >> (32 - k);
1149                 y = xa > xa0 ? *--xa : 0;
1150                 d1 = z << k | y >> (32 - k);
1151                 }
1152         else {
1153                 d0 = Exp_1 | y;
1154                 d1 = z;
1155                 }
1156 #else
1157         if (k < Ebits + 16) {
1158                 z = xa > xa0 ? *--xa : 0;
1159                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1160                 w = xa > xa0 ? *--xa : 0;
1161                 y = xa > xa0 ? *--xa : 0;
1162                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1163                 goto ret_d;
1164                 }
1165         z = xa > xa0 ? *--xa : 0;
1166         w = xa > xa0 ? *--xa : 0;
1167         k -= Ebits + 16;
1168         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1169         y = xa > xa0 ? *--xa : 0;
1170         d1 = w << k + 16 | y << k;
1171 #endif
1172  ret_d:
1173 #ifdef VAX
1174         word0(d) = d0 >> 16 | d0 << 16;
1175         word1(d) = d1 >> 16 | d1 << 16;
1176 #else
1177 #undef d0
1178 #undef d1
1179 #endif
1180         return dval(d);
1181         }
1182
1183  static Bigint *
1184 d2b
1185 #ifdef KR_headers
1186         (d, e, bits) double d; int *e, *bits;
1187 #else
1188         (double _d, int *e, int *bits)
1189 #endif
1190 {
1191         Bigint *b;
1192         int de, k;
1193         ULong *x, y, z;
1194         U d;
1195 #ifndef Sudden_Underflow
1196         int i;
1197 #endif
1198 #ifdef VAX
1199         ULong d0, d1;
1200         d0 = word0(d) >> 16 | word0(d) << 16;
1201         d1 = word1(d) >> 16 | word1(d) << 16;
1202 #else
1203 #define d0 word0(d)
1204 #define d1 word1(d)
1205 #endif
1206         dval(d) = _d;
1207
1208 #ifdef Pack_32
1209         b = Balloc(1);
1210 #else
1211         b = Balloc(2);
1212 #endif
1213         x = b->x;
1214
1215         z = d0 & Frac_mask;
1216         d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
1217 #ifdef Sudden_Underflow
1218         de = (int)(d0 >> Exp_shift);
1219 #ifndef IBM
1220         z |= Exp_msk11;
1221 #endif
1222 #else
1223         if ((de = (int)(d0 >> Exp_shift)))
1224                 z |= Exp_msk1;
1225 #endif
1226 #ifdef Pack_32
1227         if ((y = d1)) {
1228                 if ((k = lo0bits(&y))) {
1229                         x[0] = y | z << (32 - k);
1230                         z >>= k;
1231                         }
1232                 else
1233                         x[0] = y;
1234 #ifndef Sudden_Underflow
1235                 i =
1236 #endif
1237                     b->wds = (x[1] = z) ? 2 : 1;
1238                 }
1239         else {
1240 #ifdef DEBUG
1241                 if (!z)
1242                         Bug("Zero passed to d2b");
1243 #endif
1244                 k = lo0bits(&z);
1245                 x[0] = z;
1246 #ifndef Sudden_Underflow
1247                 i =
1248 #endif
1249                     b->wds = 1;
1250                 k += 32;
1251                 }
1252 #else
1253         if (y = d1) {
1254                 if (k = lo0bits(&y))
1255                         if (k >= 16) {
1256                                 x[0] = y | z << 32 - k & 0xffff;
1257                                 x[1] = z >> k - 16 & 0xffff;
1258                                 x[2] = z >> k;
1259                                 i = 2;
1260                                 }
1261                         else {
1262                                 x[0] = y & 0xffff;
1263                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
1264                                 x[2] = z >> k & 0xffff;
1265                                 x[3] = z >> k+16;
1266                                 i = 3;
1267                                 }
1268                 else {
1269                         x[0] = y & 0xffff;
1270                         x[1] = y >> 16;
1271                         x[2] = z & 0xffff;
1272                         x[3] = z >> 16;
1273                         i = 3;
1274                         }
1275                 }
1276         else {
1277 #ifdef DEBUG
1278                 if (!z)
1279                         Bug("Zero passed to d2b");
1280 #endif
1281                 k = lo0bits(&z);
1282                 if (k >= 16) {
1283                         x[0] = z;
1284                         i = 0;
1285                         }
1286                 else {
1287                         x[0] = z & 0xffff;
1288                         x[1] = z >> 16;
1289                         i = 1;
1290                         }
1291                 k += 32;
1292                 }
1293         while(!x[i])
1294                 --i;
1295         b->wds = i + 1;
1296 #endif
1297 #ifndef Sudden_Underflow
1298         if (de) {
1299 #endif
1300 #ifdef IBM
1301                 *e = (de - Bias - (P-1) << 2) + k;
1302                 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1303 #else
1304                 *e = de - Bias - (P-1) + k;
1305                 *bits = P - k;
1306 #endif
1307 #ifndef Sudden_Underflow
1308                 }
1309         else {
1310                 *e = de - Bias - (P-1) + 1 + k;
1311 #ifdef Pack_32
1312                 *bits = 32*i - hi0bits(x[i-1]);
1313 #else
1314                 *bits = (i+2)*16 - hi0bits(x[i]);
1315 #endif
1316                 }
1317 #endif
1318         return b;
1319         }
1320 #undef d0
1321 #undef d1
1322
1323  static double
1324 ratio
1325 #ifdef KR_headers
1326         (a, b) Bigint *a, *b;
1327 #else
1328         (Bigint *a, Bigint *b)
1329 #endif
1330 {
1331         U da, db;
1332         int k, ka, kb;
1333
1334         dval(da) = b2d(a, &ka);
1335         dval(db) = b2d(b, &kb);
1336 #ifdef Pack_32
1337         k = ka - kb + 32*(a->wds - b->wds);
1338 #else
1339         k = ka - kb + 16*(a->wds - b->wds);
1340 #endif
1341 #ifdef IBM
1342         if (k > 0) {
1343                 word0(da) += (k >> 2)*Exp_msk1;
1344                 if (k &= 3)
1345                         dval(da) *= 1 << k;
1346                 }
1347         else {
1348                 k = -k;
1349                 word0(db) += (k >> 2)*Exp_msk1;
1350                 if (k &= 3)
1351                         dval(db) *= 1 << k;
1352                 }
1353 #else
1354         if (k > 0)
1355                 word0(da) += k*Exp_msk1;
1356         else {
1357                 k = -k;
1358                 word0(db) += k*Exp_msk1;
1359                 }
1360 #endif
1361         return dval(da) / dval(db);
1362         }
1363
1364  static CONST double
1365 tens[] = {
1366                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1367                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1368                 1e20, 1e21, 1e22
1369 #ifdef VAX
1370                 , 1e23, 1e24
1371 #endif
1372                 };
1373
1374  static CONST double
1375 #ifdef IEEE_Arith
1376 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1377 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1378 #ifdef Avoid_Underflow
1379                 9007199254740992.*9007199254740992.e-256
1380                 /* = 2^106 * 1e-53 */
1381 #else
1382                 1e-256
1383 #endif
1384                 };
1385 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1386 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
1387 #define Scale_Bit 0x10
1388 #define n_bigtens 5
1389 #else
1390 #ifdef IBM
1391 bigtens[] = { 1e16, 1e32, 1e64 };
1392 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1393 #define n_bigtens 3
1394 #else
1395 bigtens[] = { 1e16, 1e32 };
1396 static CONST double tinytens[] = { 1e-16, 1e-32 };
1397 #define n_bigtens 2
1398 #endif
1399 #endif
1400
1401 #ifdef INFNAN_CHECK
1402
1403 #ifndef NAN_WORD0
1404 #define NAN_WORD0 0x7ff80000
1405 #endif
1406
1407 #ifndef NAN_WORD1
1408 #define NAN_WORD1 0
1409 #endif
1410
1411  static int
1412 match
1413 #ifdef KR_headers
1414         (sp, t) char **sp, *t;
1415 #else
1416         (CONST char **sp, char *t)
1417 #endif
1418 {
1419         int c, d;
1420         CONST char *s = *sp;
1421
1422         while((d = *t++)) {
1423                 if ((c = *++s) >= 'A' && c <= 'Z')
1424                         c += 'a' - 'A';
1425                 if (c != d)
1426                         return 0;
1427                 }
1428         *sp = s + 1;
1429         return 1;
1430         }
1431
1432 #ifndef No_Hex_NaN
1433  static void
1434 hexnan
1435 #ifdef KR_headers
1436         (rvp, sp) double *rvp; CONST char **sp;
1437 #else
1438         (U *rvp, CONST char **sp)
1439 #endif
1440 {
1441         ULong c, x[2];
1442         CONST char *s;
1443         int havedig, udx0, xshift;
1444
1445         x[0] = x[1] = 0;
1446         havedig = xshift = 0;
1447         udx0 = 1;
1448         s = *sp;
1449         /* allow optional initial 0x or 0X */
1450         while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1451                 ++s;
1452         if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1453                 s += 2;
1454         while((c = *(CONST unsigned char*)++s)) {
1455                 if (c >= '0' && c <= '9')
1456                         c -= '0';
1457                 else if (c >= 'a' && c <= 'f')
1458                         c += 10 - 'a';
1459                 else if (c >= 'A' && c <= 'F')
1460                         c += 10 - 'A';
1461                 else if (c <= ' ') {
1462                         if (udx0 && havedig) {
1463                                 udx0 = 0;
1464                                 xshift = 1;
1465                                 }
1466                         continue;
1467                         }
1468 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1469                 else if (/*(*/ c == ')' && havedig) {
1470                         *sp = s + 1;
1471                         break;
1472                         }
1473                 else
1474                         return; /* invalid form: don't change *sp */
1475 #else
1476                 else {
1477                         do {
1478                                 if (/*(*/ c == ')') {
1479                                         *sp = s + 1;
1480                                         break;
1481                                         }
1482                         } while((c = *++s));
1483                         break;
1484                 }
1485 #endif
1486                 havedig = 1;
1487                 if (xshift) {
1488                         xshift = 0;
1489                         x[0] = x[1];
1490                         x[1] = 0;
1491                         }
1492                 if (udx0)
1493                         x[0] = (x[0] << 4) | (x[1] >> 28);
1494                 x[1] = (x[1] << 4) | c;
1495                 }
1496         if ((x[0] &= 0xfffff) || x[1]) {
1497                 word0(*rvp) = Exp_mask | x[0];
1498                 word1(*rvp) = x[1];
1499                 }
1500         }
1501 #endif /*No_Hex_NaN*/
1502 #endif /* INFNAN_CHECK */
1503
1504  double
1505 sb_strtod
1506 #ifdef KR_headers
1507         (s00, se) CONST char *s00; char **se;
1508 #else
1509         (CONST char *s00, char **se)
1510 #endif
1511 {
1512 #ifdef Avoid_Underflow
1513         int scale;
1514 #endif
1515         int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1516                  e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1517         CONST char *s, *s0, *s1;
1518         double aadj, adj;
1519         U rv, rv0, aadj1;
1520         Long L;
1521         ULong y, z;
1522         Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1523 #ifdef SET_INEXACT
1524         int inexact, oldinexact;
1525 #endif
1526 #ifdef Honor_FLT_ROUNDS
1527         int rounding;
1528 #endif
1529 #ifdef USE_LOCALE
1530         CONST char *s2;
1531 #endif
1532
1533         sign = nz0 = nz = 0;
1534         dval(rv) = 0.;
1535         for(s = s00;;s++) switch(*s) {
1536                 case '-':
1537                         sign = 1;
1538                         /* no break */
1539                 case '+':
1540                         if (*++s)
1541                                 goto break2;
1542                         /* no break */
1543                 case 0:
1544                         goto ret0;
1545                 case '\t':
1546                 case '\n':
1547                 case '\v':
1548                 case '\f':
1549                 case '\r':
1550                 case ' ':
1551                         continue;
1552                 default:
1553                         goto break2;
1554                 }
1555  break2:
1556         if (*s == '0') {
1557                 nz0 = 1;
1558                 while(*++s == '0') ;
1559                 if (!*s)
1560                         goto ret;
1561                 }
1562         s0 = s;
1563         y = z = 0;
1564         for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1565                 if (nd < 9)
1566                         y = 10*y + c - '0';
1567                 else if (nd < 16)
1568                         z = 10*z + c - '0';
1569         nd0 = nd;
1570 #ifdef USE_LOCALE
1571         s1 = localeconv()->decimal_point;
1572         if (c == *s1) {
1573                 c = '.';
1574                 if (*++s1) {
1575                         s2 = s;
1576                         for(;;) {
1577                                 if (*++s2 != *s1) {
1578                                         c = 0;
1579                                         break;
1580                                         }
1581                                 if (!*++s1) {
1582                                         s = s2;
1583                                         break;
1584                                         }
1585                                 }
1586                         }
1587                 }
1588 #endif
1589         if (c == '.') {
1590                 c = *++s;
1591                 if (!nd) {
1592                         for(; c == '0'; c = *++s)
1593                                 nz++;
1594                         if (c > '0' && c <= '9') {
1595                                 s0 = s;
1596                                 nf += nz;
1597                                 nz = 0;
1598                                 goto have_dig;
1599                                 }
1600                         goto dig_done;
1601                         }
1602                 for(; c >= '0' && c <= '9'; c = *++s) {
1603  have_dig:
1604                         nz++;
1605                         if (c -= '0') {
1606                                 nf += nz;
1607                                 for(i = 1; i < nz; i++)
1608                                         if (nd++ < 9)
1609                                                 y *= 10;
1610                                         else if (nd <= DBL_DIG + 1)
1611                                                 z *= 10;
1612                                 if (nd++ < 9)
1613                                         y = 10*y + c;
1614                                 else if (nd <= DBL_DIG + 1)
1615                                         z = 10*z + c;
1616                                 nz = 0;
1617                                 }
1618                         }
1619                 }
1620  dig_done:
1621         e = 0;
1622         if (c == 'e' || c == 'E') {
1623                 if (!nd && !nz && !nz0) {
1624                         goto ret0;
1625                         }
1626                 s00 = s;
1627                 esign = 0;
1628                 switch(c = *++s) {
1629                         case '-':
1630                                 esign = 1;
1631                         case '+':
1632                                 c = *++s;
1633                         }
1634                 if (c >= '0' && c <= '9') {
1635                         while(c == '0')
1636                                 c = *++s;
1637                         if (c > '0' && c <= '9') {
1638                                 L = c - '0';
1639                                 s1 = s;
1640                                 while((c = *++s) >= '0' && c <= '9')
1641                                         L = 10*L + c - '0';
1642                                 if (s - s1 > 8 || L > 19999)
1643                                         /* Avoid confusion from exponents
1644                                          * so large that e might overflow.
1645                                          */
1646                                         e = 19999; /* safe for 16 bit ints */
1647                                 else
1648                                         e = (int)L;
1649                                 if (esign)
1650                                         e = -e;
1651                                 }
1652                         else
1653                                 e = 0;
1654                         }
1655                 else
1656                         s = s00;
1657                 }
1658         if (!nd) {
1659                 if (!nz && !nz0) {
1660 #ifdef INFNAN_CHECK
1661                         /* Check for Nan and Infinity */
1662                         switch(c) {
1663                           case 'i':
1664                           case 'I':
1665                                 if (match(&s,"nf")) {
1666                                         --s;
1667                                         if (!match(&s,"inity"))
1668                                                 ++s;
1669                                         word0(rv) = 0x7ff00000;
1670                                         word1(rv) = 0;
1671                                         goto ret;
1672                                         }
1673                                 break;
1674                           case 'n':
1675                           case 'N':
1676                                 if (match(&s, "an")) {
1677                                         word0(rv) = NAN_WORD0;
1678                                         word1(rv) = NAN_WORD1;
1679 #ifndef No_Hex_NaN
1680                                         if (*s == '(') /*)*/
1681                                                 hexnan(&rv, &s);
1682 #endif
1683                                         goto ret;
1684                                         }
1685                           }
1686 #endif /* INFNAN_CHECK */
1687  ret0:
1688                         s = s00;
1689                         sign = 0;
1690                         }
1691                 goto ret;
1692                 }
1693         e1 = e -= nf;
1694
1695         /* Now we have nd0 digits, starting at s0, followed by a
1696          * decimal point, followed by nd-nd0 digits.  The number we're
1697          * after is the integer represented by those digits times
1698          * 10**e */
1699
1700         if (!nd0)
1701                 nd0 = nd;
1702         k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1703         dval(rv) = y;
1704         if (k > 9) {
1705 #ifdef SET_INEXACT
1706                 if (k > DBL_DIG)
1707                         oldinexact = get_inexact();
1708 #endif
1709                 dval(rv) = tens[k - 9] * dval(rv) + z;
1710                 }
1711         bd0 = 0;
1712         if (nd <= DBL_DIG
1713 #ifndef RND_PRODQUOT
1714 #ifndef Honor_FLT_ROUNDS
1715                 && Flt_Rounds == 1
1716 #endif
1717 #endif
1718                         ) {
1719                 if (!e)
1720                         goto ret;
1721                 if (e > 0) {
1722                         if (e <= Ten_pmax) {
1723 #ifdef VAX
1724                                 goto vax_ovfl_check;
1725 #else
1726 #ifdef Honor_FLT_ROUNDS
1727                                 /* round correctly FLT_ROUNDS = 2 or 3 */
1728                                 if (sign) {
1729                                         rv = -rv;
1730                                         sign = 0;
1731                                         }
1732 #endif
1733                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1734                                 goto ret;
1735 #endif
1736                                 }
1737                         i = DBL_DIG - nd;
1738                         if (e <= Ten_pmax + i) {
1739                                 /* A fancier test would sometimes let us do
1740                                  * this for larger i values.
1741                                  */
1742 #ifdef Honor_FLT_ROUNDS
1743                                 /* round correctly FLT_ROUNDS = 2 or 3 */
1744                                 if (sign) {
1745                                         rv = -rv;
1746                                         sign = 0;
1747                                         }
1748 #endif
1749                                 e -= i;
1750                                 dval(rv) *= tens[i];
1751 #ifdef VAX
1752                                 /* VAX exponent range is so narrow we must
1753                                  * worry about overflow here...
1754                                  */
1755  vax_ovfl_check:
1756                                 word0(rv) -= P*Exp_msk1;
1757                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1758                                 if ((word0(rv) & Exp_mask)
1759                                  > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1760                                         goto ovfl;
1761                                 word0(rv) += P*Exp_msk1;
1762 #else
1763                                 /* rv = */ rounded_product(dval(rv), tens[e]);
1764 #endif
1765                                 goto ret;
1766                                 }
1767                         }
1768 #ifndef Inaccurate_Divide
1769                 else if (e >= -Ten_pmax) {
1770 #ifdef Honor_FLT_ROUNDS
1771                         /* round correctly FLT_ROUNDS = 2 or 3 */
1772                         if (sign) {
1773                                 rv = -rv;
1774                                 sign = 0;
1775                                 }
1776 #endif
1777                         /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1778                         goto ret;
1779                         }
1780 #endif
1781                 }
1782         e1 += nd - k;
1783
1784 #ifdef IEEE_Arith
1785 #ifdef SET_INEXACT
1786         inexact = 1;
1787         if (k <= DBL_DIG)
1788                 oldinexact = get_inexact();
1789 #endif
1790 #ifdef Avoid_Underflow
1791         scale = 0;
1792 #endif
1793 #ifdef Honor_FLT_ROUNDS
1794         if ((rounding = Flt_Rounds) >= 2) {
1795                 if (sign)
1796                         rounding = rounding == 2 ? 0 : 2;
1797                 else
1798                         if (rounding != 2)
1799                                 rounding = 0;
1800                 }
1801 #endif
1802 #endif /*IEEE_Arith*/
1803
1804         /* Get starting approximation = rv * 10**e1 */
1805
1806         if (e1 > 0) {
1807                 if ((i = e1 & 15))
1808                         dval(rv) *= tens[i];
1809                 if (e1 &= ~15) {
1810                         if (e1 > DBL_MAX_10_EXP) {
1811  ovfl:
1812 #ifndef NO_ERRNO
1813                                 errno = ERANGE;
1814 #endif
1815                                 /* Can't trust HUGE_VAL */
1816 #ifdef IEEE_Arith
1817 #ifdef Honor_FLT_ROUNDS
1818                                 switch(rounding) {
1819                                   case 0: /* toward 0 */
1820                                   case 3: /* toward -infinity */
1821                                         word0(rv) = Big0;
1822                                         word1(rv) = Big1;
1823                                         break;
1824                                   default:
1825                                         word0(rv) = Exp_mask;
1826                                         word1(rv) = 0;
1827                                   }
1828 #else /*Honor_FLT_ROUNDS*/
1829                                 word0(rv) = Exp_mask;
1830                                 word1(rv) = 0;
1831 #endif /*Honor_FLT_ROUNDS*/
1832 #ifdef SET_INEXACT
1833                                 /* set overflow bit */
1834                                 dval(rv0) = 1e300;
1835                                 dval(rv0) *= dval(rv0);
1836 #endif
1837 #else /*IEEE_Arith*/
1838                                 word0(rv) = Big0;
1839                                 word1(rv) = Big1;
1840 #endif /*IEEE_Arith*/
1841                                 if (bd0)
1842                                         goto retfree;
1843                                 goto ret;
1844                                 }
1845                         e1 >>= 4;
1846                         for(j = 0; e1 > 1; j++, e1 >>= 1)
1847                                 if (e1 & 1)
1848                                         dval(rv) *= bigtens[j];
1849                 /* The last multiplication could overflow. */
1850                         word0(rv) -= P*Exp_msk1;
1851                         dval(rv) *= bigtens[j];
1852                         if ((z = word0(rv) & Exp_mask)
1853                          > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1854                                 goto ovfl;
1855                         if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1856                                 /* set to largest number */
1857                                 /* (Can't trust DBL_MAX) */
1858                                 word0(rv) = Big0;
1859                                 word1(rv) = Big1;
1860                                 }
1861                         else
1862                                 word0(rv) += P*Exp_msk1;
1863                         }
1864                 }
1865         else if (e1 < 0) {
1866                 e1 = -e1;
1867                 if ((i = e1 & 15))
1868                         dval(rv) /= tens[i];
1869                 if (e1 >>= 4) {
1870                         if (e1 >= 1 << n_bigtens)
1871                                 goto undfl;
1872 #ifdef Avoid_Underflow
1873                         if (e1 & Scale_Bit)
1874                                 scale = 2*P;
1875                         for(j = 0; e1 > 0; j++, e1 >>= 1)
1876                                 if (e1 & 1)
1877                                         dval(rv) *= tinytens[j];
1878                         if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1879                                                 >> Exp_shift)) > 0) {
1880                                 /* scaled rv is denormal; zap j low bits */
1881                                 if (j >= 32) {
1882                                         word1(rv) = 0;
1883                                         if (j >= 53)
1884                                          word0(rv) = (P+2)*Exp_msk1;
1885                                         else
1886                                                 word0(rv) &= 0xffffffff << (j-32);
1887                                         }
1888                                 else
1889                                         word1(rv) &= 0xffffffff << j;
1890                                 }
1891 #else
1892                         for(j = 0; e1 > 1; j++, e1 >>= 1)
1893                                 if (e1 & 1)
1894                                         dval(rv) *= tinytens[j];
1895                         /* The last multiplication could underflow. */
1896                         dval(rv0) = dval(rv);
1897                         dval(rv) *= tinytens[j];
1898                         if (!dval(rv)) {
1899                                 dval(rv) = 2.*dval(rv0);
1900                                 dval(rv) *= tinytens[j];
1901 #endif
1902                                 if (!dval(rv)) {
1903  undfl:
1904                                         dval(rv) = 0.;
1905 #ifndef NO_ERRNO
1906                                         errno = ERANGE;
1907 #endif
1908                                         if (bd0)
1909                                                 goto retfree;
1910                                         goto ret;
1911                                         }
1912 #ifndef Avoid_Underflow
1913                                 word0(rv) = Tiny0;
1914                                 word1(rv) = Tiny1;
1915                                 /* The refinement below will clean
1916                                  * this approximation up.
1917                                  */
1918                                 }
1919 #endif
1920                         }
1921                 }
1922
1923         /* Now the hard part -- adjusting rv to the correct value.*/
1924
1925         /* Put digits into bd: true value = bd * 10^e */
1926
1927         bd0 = s2b(s0, nd0, nd, y);
1928
1929         for(;;) {
1930                 bd = Balloc(bd0->k);
1931                 Bcopy(bd, bd0);
1932                 bb = d2b(dval(rv), &bbe, &bbbits);      /* rv = bb * 2^bbe */
1933                 bs = i2b(1);
1934
1935                 if (e >= 0) {
1936                         bb2 = bb5 = 0;
1937                         bd2 = bd5 = e;
1938                         }
1939                 else {
1940                         bb2 = bb5 = -e;
1941                         bd2 = bd5 = 0;
1942                         }
1943                 if (bbe >= 0)
1944                         bb2 += bbe;
1945                 else
1946                         bd2 -= bbe;
1947                 bs2 = bb2;
1948 #ifdef Honor_FLT_ROUNDS
1949                 if (rounding != 1)
1950                         bs2++;
1951 #endif
1952 #ifdef Avoid_Underflow
1953                 j = bbe - scale;
1954                 i = j + bbbits - 1;     /* logb(rv) */
1955                 if (i < Emin)   /* denormal */
1956                         j += P - Emin;
1957                 else
1958                         j = P + 1 - bbbits;
1959 #else /*Avoid_Underflow*/
1960 #ifdef Sudden_Underflow
1961 #ifdef IBM
1962                 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1963 #else
1964                 j = P + 1 - bbbits;
1965 #endif
1966 #else /*Sudden_Underflow*/
1967                 j = bbe;
1968                 i = j + bbbits - 1;     /* logb(rv) */
1969                 if (i < Emin)   /* denormal */
1970                         j += P - Emin;
1971                 else
1972                         j = P + 1 - bbbits;
1973 #endif /*Sudden_Underflow*/
1974 #endif /*Avoid_Underflow*/
1975                 bb2 += j;
1976                 bd2 += j;
1977 #ifdef Avoid_Underflow
1978                 bd2 += scale;
1979 #endif
1980                 i = bb2 < bd2 ? bb2 : bd2;
1981                 if (i > bs2)
1982                         i = bs2;
1983                 if (i > 0) {
1984                         bb2 -= i;
1985                         bd2 -= i;
1986                         bs2 -= i;
1987                         }
1988                 if (bb5 > 0) {
1989                         bs = pow5mult(bs, bb5);
1990                         bb1 = mult(bs, bb);
1991                         Bfree(bb);
1992                         bb = bb1;
1993                         }
1994                 if (bb2 > 0)
1995                         bb = lshift(bb, bb2);
1996                 if (bd5 > 0)
1997                         bd = pow5mult(bd, bd5);
1998                 if (bd2 > 0)
1999                         bd = lshift(bd, bd2);
2000                 if (bs2 > 0)
2001                         bs = lshift(bs, bs2);
2002                 delta = diff(bb, bd);
2003                 dsign = delta->sign;
2004                 delta->sign = 0;
2005                 i = cmp(delta, bs);
2006 #ifdef Honor_FLT_ROUNDS
2007                 if (rounding != 1) {
2008                         if (i < 0) {
2009                                 /* Error is less than an ulp */
2010                                 if (!delta->x[0] && delta->wds <= 1) {
2011                                         /* exact */
2012 #ifdef SET_INEXACT
2013                                         inexact = 0;
2014 #endif
2015                                         break;
2016                                         }
2017                                 if (rounding) {
2018                                         if (dsign) {
2019                                                 adj = 1.;
2020                                                 goto apply_adj;
2021                                                 }
2022                                         }
2023                                 else if (!dsign) {
2024                                         adj = -1.;
2025                                         if (!word1(rv)
2026                                          && !(word0(rv) & Frac_mask)) {
2027                                                 y = word0(rv) & Exp_mask;
2028 #ifdef Avoid_Underflow
2029                                                 if (!scale || y > 2*P*Exp_msk1)
2030 #else
2031                                                 if (y)
2032 #endif
2033                                                   {
2034                                                   delta = lshift(delta,Log2P);
2035                                                   if (cmp(delta, bs) <= 0)
2036                                                         adj = -0.5;
2037                                                   }
2038                                                 }
2039  apply_adj:
2040 #ifdef Avoid_Underflow
2041                                         if (scale && (y = word0(rv) & Exp_mask)
2042                                                 <= 2*P*Exp_msk1)
2043                                           word0(adj) += (2*P+1)*Exp_msk1 - y;
2044 #else
2045 #ifdef Sudden_Underflow
2046                                         if ((word0(rv) & Exp_mask) <=
2047                                                         P*Exp_msk1) {
2048                                                 word0(rv) += P*Exp_msk1;
2049                                                 dval(rv) += adj*ulp(dval(rv));
2050                                                 word0(rv) -= P*Exp_msk1;
2051                                                 }
2052                                         else
2053 #endif /*Sudden_Underflow*/
2054 #endif /*Avoid_Underflow*/
2055                                         dval(rv) += adj*ulp(dval(rv));
2056                                         }
2057                                 break;
2058                                 }
2059                         adj = ratio(delta, bs);
2060                         if (adj < 1.)
2061                                 adj = 1.;
2062                         if (adj <= 0x7ffffffe) {
2063                                 /* adj = rounding ? ceil(adj) : floor(adj); */
2064                                 y = adj;
2065                                 if (y != adj) {
2066                                         if (!((rounding>>1) ^ dsign))
2067                                                 y++;
2068                                         adj = y;
2069                                         }
2070                                 }
2071 #ifdef Avoid_Underflow
2072                         if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2073                                 word0(adj) += (2*P+1)*Exp_msk1 - y;
2074 #else
2075 #ifdef Sudden_Underflow
2076                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2077                                 word0(rv) += P*Exp_msk1;
2078                                 adj *= ulp(dval(rv));
2079                                 if (dsign)
2080                                         dval(rv) += adj;
2081                                 else
2082                                         dval(rv) -= adj;
2083                                 word0(rv) -= P*Exp_msk1;
2084                                 goto cont;
2085                                 }
2086 #endif /*Sudden_Underflow*/
2087 #endif /*Avoid_Underflow*/
2088                         adj *= ulp(dval(rv));
2089                         if (dsign)
2090                                 dval(rv) += adj;
2091                         else
2092                                 dval(rv) -= adj;
2093                         goto cont;
2094                         }
2095 #endif /*Honor_FLT_ROUNDS*/
2096
2097                 if (i < 0) {
2098                         /* Error is less than half an ulp -- check for
2099                          * special case of mantissa a power of two.
2100                          */
2101                         if (dsign || word1(rv) || word0(rv) & Bndry_mask
2102 #ifdef IEEE_Arith
2103 #ifdef Avoid_Underflow
2104                          || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
2105 #else
2106                          || (word0(rv) & Exp_mask) <= Exp_msk1
2107 #endif
2108 #endif
2109                                 ) {
2110 #ifdef SET_INEXACT
2111                                 if (!delta->x[0] && delta->wds <= 1)
2112                                         inexact = 0;
2113 #endif
2114                                 break;
2115                                 }
2116                         if (!delta->x[0] && delta->wds <= 1) {
2117                                 /* exact result */
2118 #ifdef SET_INEXACT
2119                                 inexact = 0;
2120 #endif
2121                                 break;
2122                                 }
2123                         delta = lshift(delta,Log2P);
2124                         if (cmp(delta, bs) > 0)
2125                                 goto drop_down;
2126                         break;
2127                         }
2128                 if (i == 0) {
2129                         /* exactly half-way between */
2130                         if (dsign) {
2131                                 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2132                                  &&  word1(rv) == (
2133 #ifdef Avoid_Underflow
2134                         (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2135                 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2136 #endif
2137                                                    0xffffffff)) {
2138                                         /*boundary case -- increment exponent*/
2139                                         word0(rv) = (word0(rv) & Exp_mask)
2140                                                 + Exp_msk1
2141 #ifdef IBM
2142                                                 | Exp_msk1 >> 4
2143 #endif
2144                                                 ;
2145                                         word1(rv) = 0;
2146 #ifdef Avoid_Underflow
2147                                         dsign = 0;
2148 #endif
2149                                         break;
2150                                         }
2151                                 }
2152                         else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2153  drop_down:
2154                                 /* boundary case -- decrement exponent */
2155 #ifdef Sudden_Underflow /*{{*/
2156                                 L = word0(rv) & Exp_mask;
2157 #ifdef IBM
2158                                 if (L <  Exp_msk1)
2159 #else
2160 #ifdef Avoid_Underflow
2161                                 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
2162 #else
2163                                 if (L <= Exp_msk1)
2164 #endif /*Avoid_Underflow*/
2165 #endif /*IBM*/
2166                                         goto undfl;
2167                                 L -= Exp_msk1;
2168 #else /*Sudden_Underflow}{*/
2169 #ifdef Avoid_Underflow
2170                                 if (scale) {
2171                                         L = word0(rv) & Exp_mask;
2172                                         if (L <= (2*P+1)*Exp_msk1) {
2173                                                 if (L > (P+2)*Exp_msk1)
2174                                                         /* round even ==> */
2175                                                         /* accept rv */
2176                                                         break;
2177                                                 /* rv = smallest denormal */
2178                                                 goto undfl;
2179                                                 }
2180                                         }
2181 #endif /*Avoid_Underflow*/
2182                                 L = (word0(rv) & Exp_mask) - Exp_msk1;
2183 #endif /*Sudden_Underflow}}*/
2184                                 word0(rv) = L | Bndry_mask1;
2185                                 word1(rv) = 0xffffffff;
2186 #ifdef IBM
2187                                 goto cont;
2188 #else
2189                                 break;
2190 #endif
2191                                 }
2192 #ifndef ROUND_BIASED
2193                         if (!(word1(rv) & LSB))
2194                                 break;
2195 #endif
2196                         if (dsign)
2197                                 dval(rv) += ulp(dval(rv));
2198 #ifndef ROUND_BIASED
2199                         else {
2200                                 dval(rv) -= ulp(dval(rv));
2201 #ifndef Sudden_Underflow
2202                                 if (!dval(rv))
2203                                         goto undfl;
2204 #endif
2205                                 }
2206 #ifdef Avoid_Underflow
2207                         dsign = 1 - dsign;
2208 #endif
2209 #endif
2210                         break;
2211                         }
2212                 if ((aadj = ratio(delta, bs)) <= 2.) {
2213                         if (dsign)
2214                                 aadj = dval(aadj1) = 1.;
2215                         else if (word1(rv) || word0(rv) & Bndry_mask) {
2216 #ifndef Sudden_Underflow
2217                                 if (word1(rv) == Tiny1 && !word0(rv))
2218                                         goto undfl;
2219 #endif
2220                                 aadj = 1.;
2221                                 dval(aadj1) = -1.;
2222                                 }
2223                         else {
2224                                 /* special case -- power of FLT_RADIX to be */
2225                                 /* rounded down... */
2226
2227                                 if (aadj < 2./FLT_RADIX)
2228                                         aadj = 1./FLT_RADIX;
2229                                 else
2230                                         aadj *= 0.5;
2231                                 dval(aadj1) = -aadj;
2232                                 }
2233                         }
2234                 else {
2235                         aadj *= 0.5;
2236                         dval(aadj1) = dsign ? aadj : -aadj;
2237 #ifdef Check_FLT_ROUNDS
2238                         switch(Rounding) {
2239                                 case 2: /* towards +infinity */
2240                                         dval(aadj1) -= 0.5;
2241                                         break;
2242                                 case 0: /* towards 0 */
2243                                 case 3: /* towards -infinity */
2244                                         dval(aadj1) += 0.5;
2245                                 }
2246 #else
2247                         if (Flt_Rounds == 0)
2248                                 dval(aadj1) += 0.5;
2249 #endif /*Check_FLT_ROUNDS*/
2250                         }
2251                 y = word0(rv) & Exp_mask;
2252
2253                 /* Check for overflow */
2254
2255                 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2256                         dval(rv0) = dval(rv);
2257                         word0(rv) -= P*Exp_msk1;
2258                         adj = dval(aadj1) * ulp(dval(rv));
2259                         dval(rv) += adj;
2260                         if ((word0(rv) & Exp_mask) >=
2261                                         Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2262                                 if (word0(rv0) == Big0 && word1(rv0) == Big1)
2263                                         goto ovfl;
2264                                 word0(rv) = Big0;
2265                                 word1(rv) = Big1;
2266                                 goto cont;
2267                                 }
2268                         else
2269                                 word0(rv) += P*Exp_msk1;
2270                         }
2271                 else {
2272 #ifdef Avoid_Underflow
2273                         if (scale && y <= 2*P*Exp_msk1) {
2274                                 if (aadj <= 0x7fffffff) {
2275                                         if ((z = (uint32)aadj) <= 0)
2276                                                 z = 1;
2277                                         aadj = z;
2278                                         dval(aadj1) = dsign ? aadj : -aadj;
2279                                         }
2280                                 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
2281                                 }
2282                         adj = dval(aadj1) * ulp(dval(rv));
2283                         dval(rv) += adj;
2284 #else
2285 #ifdef Sudden_Underflow
2286                         if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2287                                 dval(rv0) = dval(rv);
2288                                 word0(rv) += P*Exp_msk1;
2289                                 adj = aadj1 * ulp(dval(rv));
2290                                 dval(rv) += adj;
2291 #ifdef IBM
2292                                 if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
2293 #else
2294                                 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2295 #endif
2296                                         {
2297                                         if (word0(rv0) == Tiny0
2298                                          && word1(rv0) == Tiny1)
2299                                                 goto undfl;
2300                                         word0(rv) = Tiny0;
2301                                         word1(rv) = Tiny1;
2302                                         goto cont;
2303                                         }
2304                                 else
2305                                         word0(rv) -= P*Exp_msk1;
2306                                 }
2307                         else {
2308                                 adj = aadj1 * ulp(dval(rv));
2309                                 dval(rv) += adj;
2310                                 }
2311 #else /*Sudden_Underflow*/
2312                         /* Compute adj so that the IEEE rounding rules will
2313                          * correctly round rv + adj in some half-way cases.
2314                          * If rv * ulp(rv) is denormalized (i.e.,
2315                          * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2316                          * trouble from bits lost to denormalization;
2317                          * example: 1.2e-307 .
2318                          */
2319                         if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
2320                                 aadj1 = (double)(int)(aadj + 0.5);
2321                                 if (!dsign)
2322                                         aadj1 = -aadj1;
2323                                 }
2324                         adj = aadj1 * ulp(dval(rv));
2325                         dval(rv) += adj;
2326 #endif /*Sudden_Underflow*/
2327 #endif /*Avoid_Underflow*/
2328                         }
2329                 z = word0(rv) & Exp_mask;
2330 #ifndef SET_INEXACT
2331 #ifdef Avoid_Underflow
2332                 if (!scale)
2333 #endif
2334                 if (y == z) {
2335                         /* Can we stop now? */
2336                         L = (Long)aadj;
2337                         aadj -= L;
2338                         /* The tolerances below are conservative. */
2339                         if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2340                                 if (aadj < .4999999 || aadj > .5000001)
2341                                         break;
2342                                 }
2343                         else if (aadj < .4999999/FLT_RADIX)
2344                                 break;
2345                         }
2346 #endif
2347  cont:
2348                 Bfree(bb);
2349                 Bfree(bd);
2350                 Bfree(bs);
2351                 Bfree(delta);
2352                 }
2353 #ifdef SET_INEXACT
2354         if (inexact) {
2355                 if (!oldinexact) {
2356                         word0(rv0) = Exp_1 + (70 << Exp_shift);
2357                         word1(rv0) = 0;
2358                         dval(rv0) += 1.;
2359                         }
2360                 }
2361         else if (!oldinexact)
2362                 clear_inexact();
2363 #endif
2364 #ifdef Avoid_Underflow
2365         if (scale) {
2366                 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2367                 word1(rv0) = 0;
2368                 dval(rv) *= dval(rv0);
2369 #ifndef NO_ERRNO
2370                 /* try to avoid the bug of testing an 8087 register value */
2371                 if (word0(rv) == 0 && word1(rv) == 0)
2372                         errno = ERANGE;
2373 #endif
2374                 }
2375 #endif /* Avoid_Underflow */
2376 #ifdef SET_INEXACT
2377         if (inexact && !(word0(rv) & Exp_mask)) {
2378                 /* set underflow bit */
2379                 dval(rv0) = 1e-300;
2380                 dval(rv0) *= dval(rv0);
2381                 }
2382 #endif
2383  retfree:
2384         Bfree(bb);
2385         Bfree(bd);
2386         Bfree(bs);
2387         Bfree(bd0);
2388         Bfree(delta);
2389  ret:
2390         if (se)
2391                 *se = (char *)s;
2392         return sign ? -dval(rv) : dval(rv);
2393         }