PageRenderTime 140ms CodeModel.GetById 29ms app.highlight 94ms RepoModel.GetById 1ms app.codeStats 1ms

/js/src/dtoa.c

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

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