PageRenderTime 327ms CodeModel.GetById 146ms app.highlight 156ms RepoModel.GetById 1ms app.codeStats 1ms

/js/lib/Socket.IO-node/support/expresso/deps/jscoverage/js/dtoa.c

http://github.com/onedayitwillmake/RealtimeMultiplayerNodeJs
C | 3165 lines | 2670 code | 139 blank | 356 comment | 505 complexity | 6dcfb87deaf161136497bc671dec37e8 MD5 | raw file

Large files files are truncated, but you can click here to view the full file

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

Large files files are truncated, but you can click here to view the full file