1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/js/src/dtoa.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,3248 @@
1.4 +/* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
1.5 +/****************************************************************
1.6 + *
1.7 + * The author of this software is David M. Gay.
1.8 + *
1.9 + * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
1.10 + *
1.11 + * Permission to use, copy, modify, and distribute this software for any
1.12 + * purpose without fee is hereby granted, provided that this entire notice
1.13 + * is included in all copies of any software which is or includes a copy
1.14 + * or modification of this software and in all copies of the supporting
1.15 + * documentation for such software.
1.16 + *
1.17 + * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
1.18 + * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
1.19 + * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
1.20 + * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
1.21 + *
1.22 + ***************************************************************/
1.23 +
1.24 +/* Please send bug reports to David M. Gay (dmg at acm dot org,
1.25 + * with " at " changed at "@" and " dot " changed to "."). */
1.26 +
1.27 +/* On a machine with IEEE extended-precision registers, it is
1.28 + * necessary to specify double-precision (53-bit) rounding precision
1.29 + * before invoking strtod or dtoa. If the machine uses (the equivalent
1.30 + * of) Intel 80x87 arithmetic, the call
1.31 + * _control87(PC_53, MCW_PC);
1.32 + * does this with many compilers. Whether this or another call is
1.33 + * appropriate depends on the compiler; for this to work, it may be
1.34 + * necessary to #include "float.h" or another system-dependent header
1.35 + * file.
1.36 + */
1.37 +
1.38 +/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
1.39 + *
1.40 + * This strtod returns a nearest machine number to the input decimal
1.41 + * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
1.42 + * broken by the IEEE round-even rule. Otherwise ties are broken by
1.43 + * biased rounding (add half and chop).
1.44 + *
1.45 + * Inspired loosely by William D. Clinger's paper "How to Read Floating
1.46 + * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
1.47 + *
1.48 + * Modifications:
1.49 + *
1.50 + * 1. We only require IEEE, IBM, or VAX double-precision
1.51 + * arithmetic (not IEEE double-extended).
1.52 + * 2. We get by with floating-point arithmetic in a case that
1.53 + * Clinger missed -- when we're computing d * 10^n
1.54 + * for a small integer d and the integer n is not too
1.55 + * much larger than 22 (the maximum integer k for which
1.56 + * we can represent 10^k exactly), we may be able to
1.57 + * compute (d*10^k) * 10^(e-k) with just one roundoff.
1.58 + * 3. Rather than a bit-at-a-time adjustment of the binary
1.59 + * result in the hard case, we use floating-point
1.60 + * arithmetic to determine the adjustment to within
1.61 + * one bit; only in really hard cases do we need to
1.62 + * compute a second residual.
1.63 + * 4. Because of 3., we don't need a large table of powers of 10
1.64 + * for ten-to-e (just some small tables, e.g. of 10^k
1.65 + * for 0 <= k <= 22).
1.66 + */
1.67 +
1.68 +/*
1.69 + * #define IEEE_8087 for IEEE-arithmetic machines where the least
1.70 + * significant byte has the lowest address.
1.71 + * #define IEEE_MC68k for IEEE-arithmetic machines where the most
1.72 + * significant byte has the lowest address.
1.73 + * #define Long int on machines with 32-bit ints and 64-bit longs.
1.74 + * #define IBM for IBM mainframe-style floating-point arithmetic.
1.75 + * #define VAX for VAX-style floating-point arithmetic (D_floating).
1.76 + * #define No_leftright to omit left-right logic in fast floating-point
1.77 + * computation of dtoa.
1.78 + * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
1.79 + * and strtod and dtoa should round accordingly.
1.80 + * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
1.81 + * and Honor_FLT_ROUNDS is not #defined.
1.82 + * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
1.83 + * that use extended-precision instructions to compute rounded
1.84 + * products and quotients) with IBM.
1.85 + * #define ROUND_BIASED for IEEE-format with biased rounding.
1.86 + * #define Inaccurate_Divide for IEEE-format with correctly rounded
1.87 + * products but inaccurate quotients, e.g., for Intel i860.
1.88 + * #define NO_LONG_LONG on machines that do not have a "long long"
1.89 + * integer type (of >= 64 bits). On such machines, you can
1.90 + * #define Just_16 to store 16 bits per 32-bit Long when doing
1.91 + * high-precision integer arithmetic. Whether this speeds things
1.92 + * up or slows things down depends on the machine and the number
1.93 + * being converted. If long long is available and the name is
1.94 + * something other than "long long", #define Llong to be the name,
1.95 + * and if "unsigned Llong" does not work as an unsigned version of
1.96 + * Llong, #define #ULLong to be the corresponding unsigned type.
1.97 + * #define KR_headers for old-style C function headers.
1.98 + * #define Bad_float_h if your system lacks a float.h or if it does not
1.99 + * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
1.100 + * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
1.101 + * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
1.102 + * if memory is available and otherwise does something you deem
1.103 + * appropriate. If MALLOC is undefined, malloc will be invoked
1.104 + * directly -- and assumed always to succeed. Similarly, if you
1.105 + * want something other than the system's free() to be called to
1.106 + * recycle memory acquired from MALLOC, #define FREE to be the
1.107 + * name of the alternate routine. (Unless you #define
1.108 + * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only
1.109 + * called in pathological cases, e.g., in a dtoa call after a dtoa
1.110 + * return in mode 3 with thousands of digits requested.)
1.111 + * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
1.112 + * memory allocations from a private pool of memory when possible.
1.113 + * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
1.114 + * unless #defined to be a different length. This default length
1.115 + * suffices to get rid of MALLOC calls except for unusual cases,
1.116 + * such as decimal-to-binary conversion of a very long string of
1.117 + * digits. The longest string dtoa can return is about 751 bytes
1.118 + * long. For conversions by strtod of strings of 800 digits and
1.119 + * all dtoa conversions in single-threaded executions with 8-byte
1.120 + * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
1.121 + * pointers, PRIVATE_MEM >= 7112 appears adequate.
1.122 + * #define MULTIPLE_THREADS if the system offers preemptively scheduled
1.123 + * multiple threads. In this case, you must provide (or suitably
1.124 + * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
1.125 + * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
1.126 + * in pow5mult, ensures lazy evaluation of only one copy of high
1.127 + * powers of 5; omitting this lock would introduce a small
1.128 + * probability of wasting memory, but would otherwise be harmless.)
1.129 + * You must also invoke freedtoa(s) to free the value s returned by
1.130 + * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
1.131 + * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
1.132 + * avoids underflows on inputs whose result does not underflow.
1.133 + * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
1.134 + * floating-point numbers and flushes underflows to zero rather
1.135 + * than implementing gradual underflow, then you must also #define
1.136 + * Sudden_Underflow.
1.137 + * #define USE_LOCALE to use the current locale's decimal_point value.
1.138 + * #define SET_INEXACT if IEEE arithmetic is being used and extra
1.139 + * computation should be done to set the inexact flag when the
1.140 + * result is inexact and avoid setting inexact when the result
1.141 + * is exact. In this case, dtoa.c must be compiled in
1.142 + * an environment, perhaps provided by #include "dtoa.c" in a
1.143 + * suitable wrapper, that defines two functions,
1.144 + * int get_inexact(void);
1.145 + * void clear_inexact(void);
1.146 + * such that get_inexact() returns a nonzero value if the
1.147 + * inexact bit is already set, and clear_inexact() sets the
1.148 + * inexact bit to 0. When SET_INEXACT is #defined, strtod
1.149 + * also does extra computations to set the underflow and overflow
1.150 + * flags when appropriate (i.e., when the result is tiny and
1.151 + * inexact or when it is a numeric value rounded to +-infinity).
1.152 + * #define NO_ERRNO if strtod should not assign errno = ERANGE when
1.153 + * the result overflows to +-Infinity or underflows to 0.
1.154 + * #define NO_GLOBAL_STATE to avoid defining any non-const global or
1.155 + * static variables. Instead the necessary state is stored in an
1.156 + * opaque struct, DtoaState, a pointer to which must be passed to
1.157 + * every entry point. Two new functions are added to the API:
1.158 + * DtoaState *newdtoa(void);
1.159 + * void destroydtoa(DtoaState *);
1.160 + */
1.161 +
1.162 +#ifndef Long
1.163 +#define Long long
1.164 +#endif
1.165 +#ifndef ULong
1.166 +typedef unsigned Long ULong;
1.167 +#endif
1.168 +
1.169 +#ifdef DEBUG
1.170 +#include <stdio.h>
1.171 +#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
1.172 +#endif
1.173 +
1.174 +#include <stdlib.h>
1.175 +#include <string.h>
1.176 +
1.177 +#ifdef USE_LOCALE
1.178 +#include <locale.h>
1.179 +#endif
1.180 +
1.181 +#ifdef MALLOC
1.182 +#ifdef KR_headers
1.183 +extern char *MALLOC();
1.184 +#else
1.185 +extern void *MALLOC(size_t);
1.186 +#endif
1.187 +#else
1.188 +#define MALLOC malloc
1.189 +#endif
1.190 +
1.191 +#ifndef FREE
1.192 +#define FREE free
1.193 +#endif
1.194 +
1.195 +#ifndef Omit_Private_Memory
1.196 +#ifndef PRIVATE_MEM
1.197 +#define PRIVATE_MEM 2304
1.198 +#endif
1.199 +#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
1.200 +#endif
1.201 +
1.202 +#undef IEEE_Arith
1.203 +#undef Avoid_Underflow
1.204 +#ifdef IEEE_MC68k
1.205 +#define IEEE_Arith
1.206 +#endif
1.207 +#ifdef IEEE_8087
1.208 +#define IEEE_Arith
1.209 +#endif
1.210 +
1.211 +#include <errno.h>
1.212 +
1.213 +#ifdef Bad_float_h
1.214 +
1.215 +#ifdef IEEE_Arith
1.216 +#define DBL_DIG 15
1.217 +#define DBL_MAX_10_EXP 308
1.218 +#define DBL_MAX_EXP 1024
1.219 +#define FLT_RADIX 2
1.220 +#endif /*IEEE_Arith*/
1.221 +
1.222 +#ifdef IBM
1.223 +#define DBL_DIG 16
1.224 +#define DBL_MAX_10_EXP 75
1.225 +#define DBL_MAX_EXP 63
1.226 +#define FLT_RADIX 16
1.227 +#define DBL_MAX 7.2370055773322621e+75
1.228 +#endif
1.229 +
1.230 +#ifdef VAX
1.231 +#define DBL_DIG 16
1.232 +#define DBL_MAX_10_EXP 38
1.233 +#define DBL_MAX_EXP 127
1.234 +#define FLT_RADIX 2
1.235 +#define DBL_MAX 1.7014118346046923e+38
1.236 +#endif
1.237 +
1.238 +#ifndef LONG_MAX
1.239 +#define LONG_MAX 2147483647
1.240 +#endif
1.241 +
1.242 +#else /* ifndef Bad_float_h */
1.243 +#include <float.h>
1.244 +#endif /* Bad_float_h */
1.245 +
1.246 +#ifndef __MATH_H__
1.247 +#include <math.h>
1.248 +#endif
1.249 +
1.250 +#ifndef CONST
1.251 +#ifdef KR_headers
1.252 +#define CONST /* blank */
1.253 +#else
1.254 +#define CONST const
1.255 +#endif
1.256 +#endif
1.257 +
1.258 +#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
1.259 +Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
1.260 +#endif
1.261 +
1.262 +typedef union { double d; ULong L[2]; } U;
1.263 +
1.264 +#define dval(x) ((x).d)
1.265 +#ifdef IEEE_8087
1.266 +#define word0(x) ((x).L[1])
1.267 +#define word1(x) ((x).L[0])
1.268 +#else
1.269 +#define word0(x) ((x).L[0])
1.270 +#define word1(x) ((x).L[1])
1.271 +#endif
1.272 +
1.273 +/* The following definition of Storeinc is appropriate for MIPS processors.
1.274 + * An alternative that might be better on some machines is
1.275 + * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
1.276 + */
1.277 +#if defined(IEEE_8087) + defined(VAX)
1.278 +#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
1.279 +((unsigned short *)a)[0] = (unsigned short)c, a++)
1.280 +#else
1.281 +#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
1.282 +((unsigned short *)a)[1] = (unsigned short)c, a++)
1.283 +#endif
1.284 +
1.285 +/* #define P DBL_MANT_DIG */
1.286 +/* Ten_pmax = floor(P*log(2)/log(5)) */
1.287 +/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
1.288 +/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
1.289 +/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
1.290 +
1.291 +#ifdef IEEE_Arith
1.292 +#define Exp_shift 20
1.293 +#define Exp_shift1 20
1.294 +#define Exp_msk1 0x100000
1.295 +#define Exp_msk11 0x100000
1.296 +#define Exp_mask 0x7ff00000
1.297 +#define P 53
1.298 +#define Bias 1023
1.299 +#define Emin (-1022)
1.300 +#define Exp_1 0x3ff00000
1.301 +#define Exp_11 0x3ff00000
1.302 +#define Ebits 11
1.303 +#define Frac_mask 0xfffff
1.304 +#define Frac_mask1 0xfffff
1.305 +#define Ten_pmax 22
1.306 +#define Bletch 0x10
1.307 +#define Bndry_mask 0xfffff
1.308 +#define Bndry_mask1 0xfffff
1.309 +#define LSB 1
1.310 +#define Sign_bit 0x80000000
1.311 +#define Log2P 1
1.312 +#define Tiny0 0
1.313 +#define Tiny1 1
1.314 +#define Quick_max 14
1.315 +#define Int_max 14
1.316 +#ifndef NO_IEEE_Scale
1.317 +#define Avoid_Underflow
1.318 +#ifdef Flush_Denorm /* debugging option */
1.319 +#undef Sudden_Underflow
1.320 +#endif
1.321 +#endif
1.322 +
1.323 +#ifndef Flt_Rounds
1.324 +#ifdef FLT_ROUNDS
1.325 +#define Flt_Rounds FLT_ROUNDS
1.326 +#else
1.327 +#define Flt_Rounds 1
1.328 +#endif
1.329 +#endif /*Flt_Rounds*/
1.330 +
1.331 +#ifdef Honor_FLT_ROUNDS
1.332 +#define Rounding rounding
1.333 +#undef Check_FLT_ROUNDS
1.334 +#define Check_FLT_ROUNDS
1.335 +#else
1.336 +#define Rounding Flt_Rounds
1.337 +#endif
1.338 +
1.339 +#else /* ifndef IEEE_Arith */
1.340 +#undef Check_FLT_ROUNDS
1.341 +#undef Honor_FLT_ROUNDS
1.342 +#undef SET_INEXACT
1.343 +#undef Sudden_Underflow
1.344 +#define Sudden_Underflow
1.345 +#ifdef IBM
1.346 +#undef Flt_Rounds
1.347 +#define Flt_Rounds 0
1.348 +#define Exp_shift 24
1.349 +#define Exp_shift1 24
1.350 +#define Exp_msk1 0x1000000
1.351 +#define Exp_msk11 0x1000000
1.352 +#define Exp_mask 0x7f000000
1.353 +#define P 14
1.354 +#define Bias 65
1.355 +#define Exp_1 0x41000000
1.356 +#define Exp_11 0x41000000
1.357 +#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
1.358 +#define Frac_mask 0xffffff
1.359 +#define Frac_mask1 0xffffff
1.360 +#define Bletch 4
1.361 +#define Ten_pmax 22
1.362 +#define Bndry_mask 0xefffff
1.363 +#define Bndry_mask1 0xffffff
1.364 +#define LSB 1
1.365 +#define Sign_bit 0x80000000
1.366 +#define Log2P 4
1.367 +#define Tiny0 0x100000
1.368 +#define Tiny1 0
1.369 +#define Quick_max 14
1.370 +#define Int_max 15
1.371 +#else /* VAX */
1.372 +#undef Flt_Rounds
1.373 +#define Flt_Rounds 1
1.374 +#define Exp_shift 23
1.375 +#define Exp_shift1 7
1.376 +#define Exp_msk1 0x80
1.377 +#define Exp_msk11 0x800000
1.378 +#define Exp_mask 0x7f80
1.379 +#define P 56
1.380 +#define Bias 129
1.381 +#define Exp_1 0x40800000
1.382 +#define Exp_11 0x4080
1.383 +#define Ebits 8
1.384 +#define Frac_mask 0x7fffff
1.385 +#define Frac_mask1 0xffff007f
1.386 +#define Ten_pmax 24
1.387 +#define Bletch 2
1.388 +#define Bndry_mask 0xffff007f
1.389 +#define Bndry_mask1 0xffff007f
1.390 +#define LSB 0x10000
1.391 +#define Sign_bit 0x8000
1.392 +#define Log2P 1
1.393 +#define Tiny0 0x80
1.394 +#define Tiny1 0
1.395 +#define Quick_max 15
1.396 +#define Int_max 15
1.397 +#endif /* IBM, VAX */
1.398 +#endif /* IEEE_Arith */
1.399 +
1.400 +#ifndef IEEE_Arith
1.401 +#define ROUND_BIASED
1.402 +#endif
1.403 +
1.404 +#ifdef RND_PRODQUOT
1.405 +#define rounded_product(a,b) a = rnd_prod(a, b)
1.406 +#define rounded_quotient(a,b) a = rnd_quot(a, b)
1.407 +#ifdef KR_headers
1.408 +extern double rnd_prod(), rnd_quot();
1.409 +#else
1.410 +extern double rnd_prod(double, double), rnd_quot(double, double);
1.411 +#endif
1.412 +#else
1.413 +#define rounded_product(a,b) a *= b
1.414 +#define rounded_quotient(a,b) a /= b
1.415 +#endif
1.416 +
1.417 +#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
1.418 +#define Big1 0xffffffff
1.419 +
1.420 +#ifndef Pack_32
1.421 +#define Pack_32
1.422 +#endif
1.423 +
1.424 +#ifdef KR_headers
1.425 +#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
1.426 +#else
1.427 +#define FFFFFFFF 0xffffffffUL
1.428 +#endif
1.429 +
1.430 +#ifdef NO_LONG_LONG
1.431 +#undef ULLong
1.432 +#ifdef Just_16
1.433 +#undef Pack_32
1.434 +/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
1.435 + * This makes some inner loops simpler and sometimes saves work
1.436 + * during multiplications, but it often seems to make things slightly
1.437 + * slower. Hence the default is now to store 32 bits per Long.
1.438 + */
1.439 +#endif
1.440 +#else /* long long available */
1.441 +#ifndef Llong
1.442 +#define Llong long long
1.443 +#endif
1.444 +#ifndef ULLong
1.445 +#define ULLong unsigned Llong
1.446 +#endif
1.447 +#endif /* NO_LONG_LONG */
1.448 +
1.449 +#ifndef MULTIPLE_THREADS
1.450 +#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
1.451 +#define FREE_DTOA_LOCK(n) /*nothing*/
1.452 +#endif
1.453 +
1.454 +#define Kmax 7
1.455 +
1.456 + struct
1.457 +Bigint {
1.458 + struct Bigint *next;
1.459 + int k, maxwds, sign, wds;
1.460 + ULong x[1];
1.461 + };
1.462 +
1.463 + typedef struct Bigint Bigint;
1.464 +
1.465 +#ifdef NO_GLOBAL_STATE
1.466 +#ifdef MULTIPLE_THREADS
1.467 +#error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"
1.468 +#endif
1.469 + struct
1.470 +DtoaState {
1.471 +#define DECLARE_GLOBAL_STATE /* nothing */
1.472 +#else
1.473 +#define DECLARE_GLOBAL_STATE static
1.474 +#endif
1.475 +
1.476 + DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1];
1.477 + DECLARE_GLOBAL_STATE Bigint *p5s;
1.478 +#ifndef Omit_Private_Memory
1.479 + DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem];
1.480 + DECLARE_GLOBAL_STATE double *pmem_next
1.481 +#ifndef NO_GLOBAL_STATE
1.482 + = private_mem
1.483 +#endif
1.484 + ;
1.485 +#endif
1.486 +#ifdef NO_GLOBAL_STATE
1.487 + };
1.488 + typedef struct DtoaState DtoaState;
1.489 +#ifdef KR_headers
1.490 +#define STATE_PARAM state,
1.491 +#define STATE_PARAM_DECL DtoaState *state;
1.492 +#else
1.493 +#define STATE_PARAM DtoaState *state,
1.494 +#endif
1.495 +#define PASS_STATE state,
1.496 +#define GET_STATE(field) (state->field)
1.497 +
1.498 + static DtoaState *
1.499 +newdtoa(void)
1.500 +{
1.501 + DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState));
1.502 + if (state) {
1.503 + memset(state, 0, sizeof(DtoaState));
1.504 +#ifndef Omit_Private_Memory
1.505 + state->pmem_next = state->private_mem;
1.506 +#endif
1.507 + }
1.508 + return state;
1.509 +}
1.510 +
1.511 + static void
1.512 +destroydtoa
1.513 +#ifdef KR_headers
1.514 + (state) STATE_PARAM_DECL
1.515 +#else
1.516 + (DtoaState *state)
1.517 +#endif
1.518 +{
1.519 + int i;
1.520 + Bigint *v, *next;
1.521 +
1.522 + for (i = 0; i <= Kmax; i++) {
1.523 + for (v = GET_STATE(freelist)[i]; v; v = next) {
1.524 + next = v->next;
1.525 +#ifndef Omit_Private_Memory
1.526 + if ((double*)v < GET_STATE(private_mem) ||
1.527 + (double*)v >= GET_STATE(private_mem) + PRIVATE_mem)
1.528 +#endif
1.529 + FREE((void*)v);
1.530 + }
1.531 + }
1.532 + FREE((void *)state);
1.533 +}
1.534 +
1.535 +#else
1.536 +#define STATE_PARAM /* nothing */
1.537 +#define STATE_PARAM_DECL /* nothing */
1.538 +#define PASS_STATE /* nothing */
1.539 +#define GET_STATE(name) name
1.540 +#endif
1.541 +
1.542 + static Bigint *
1.543 +Balloc
1.544 +#ifdef KR_headers
1.545 + (STATE_PARAM k) STATE_PARAM_DECL int k;
1.546 +#else
1.547 + (STATE_PARAM int k)
1.548 +#endif
1.549 +{
1.550 + int x;
1.551 + Bigint *rv;
1.552 +#ifndef Omit_Private_Memory
1.553 + size_t len;
1.554 +#endif
1.555 +
1.556 + ACQUIRE_DTOA_LOCK(0);
1.557 + /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
1.558 + /* but this case seems very unlikely. */
1.559 + if (k <= Kmax && (rv = GET_STATE(freelist)[k]))
1.560 + GET_STATE(freelist)[k] = rv->next;
1.561 + else {
1.562 + x = 1 << k;
1.563 +#ifdef Omit_Private_Memory
1.564 + rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
1.565 +#else
1.566 + len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
1.567 + /sizeof(double);
1.568 + if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) {
1.569 + rv = (Bigint*)GET_STATE(pmem_next);
1.570 + GET_STATE(pmem_next) += len;
1.571 + }
1.572 + else
1.573 + rv = (Bigint*)MALLOC(len*sizeof(double));
1.574 +#endif
1.575 + rv->k = k;
1.576 + rv->maxwds = x;
1.577 + }
1.578 + FREE_DTOA_LOCK(0);
1.579 + rv->sign = rv->wds = 0;
1.580 + return rv;
1.581 + }
1.582 +
1.583 + static void
1.584 +Bfree
1.585 +#ifdef KR_headers
1.586 + (STATE_PARAM v) STATE_PARAM_DECL Bigint *v;
1.587 +#else
1.588 + (STATE_PARAM Bigint *v)
1.589 +#endif
1.590 +{
1.591 + if (v) {
1.592 + if (v->k > Kmax)
1.593 + FREE((void*)v);
1.594 + else {
1.595 + ACQUIRE_DTOA_LOCK(0);
1.596 + v->next = GET_STATE(freelist)[v->k];
1.597 + GET_STATE(freelist)[v->k] = v;
1.598 + FREE_DTOA_LOCK(0);
1.599 + }
1.600 + }
1.601 + }
1.602 +
1.603 +#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
1.604 +y->wds*sizeof(Long) + 2*sizeof(int))
1.605 +
1.606 + static Bigint *
1.607 +multadd
1.608 +#ifdef KR_headers
1.609 + (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a;
1.610 +#else
1.611 + (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */
1.612 +#endif
1.613 +{
1.614 + int i, wds;
1.615 +#ifdef ULLong
1.616 + ULong *x;
1.617 + ULLong carry, y;
1.618 +#else
1.619 + ULong carry, *x, y;
1.620 +#ifdef Pack_32
1.621 + ULong xi, z;
1.622 +#endif
1.623 +#endif
1.624 + Bigint *b1;
1.625 +
1.626 + wds = b->wds;
1.627 + x = b->x;
1.628 + i = 0;
1.629 + carry = a;
1.630 + do {
1.631 +#ifdef ULLong
1.632 + y = *x * (ULLong)m + carry;
1.633 + carry = y >> 32;
1.634 + *x++ = (ULong) y & FFFFFFFF;
1.635 +#else
1.636 +#ifdef Pack_32
1.637 + xi = *x;
1.638 + y = (xi & 0xffff) * m + carry;
1.639 + z = (xi >> 16) * m + (y >> 16);
1.640 + carry = z >> 16;
1.641 + *x++ = (z << 16) + (y & 0xffff);
1.642 +#else
1.643 + y = *x * m + carry;
1.644 + carry = y >> 16;
1.645 + *x++ = y & 0xffff;
1.646 +#endif
1.647 +#endif
1.648 + }
1.649 + while(++i < wds);
1.650 + if (carry) {
1.651 + if (wds >= b->maxwds) {
1.652 + b1 = Balloc(PASS_STATE b->k+1);
1.653 + Bcopy(b1, b);
1.654 + Bfree(PASS_STATE b);
1.655 + b = b1;
1.656 + }
1.657 + b->x[wds++] = (ULong) carry;
1.658 + b->wds = wds;
1.659 + }
1.660 + return b;
1.661 + }
1.662 +
1.663 + static Bigint *
1.664 +s2b
1.665 +#ifdef KR_headers
1.666 + (STATE_PARAM s, nd0, nd, y9) STATE_PARAM_DECL CONST char *s; int nd0, nd; ULong y9;
1.667 +#else
1.668 + (STATE_PARAM CONST char *s, int nd0, int nd, ULong y9)
1.669 +#endif
1.670 +{
1.671 + Bigint *b;
1.672 + int i, k;
1.673 + Long x, y;
1.674 +
1.675 + x = (nd + 8) / 9;
1.676 + for(k = 0, y = 1; x > y; y <<= 1, k++) ;
1.677 +#ifdef Pack_32
1.678 + b = Balloc(PASS_STATE k);
1.679 + b->x[0] = y9;
1.680 + b->wds = 1;
1.681 +#else
1.682 + b = Balloc(PASS_STATE k+1);
1.683 + b->x[0] = y9 & 0xffff;
1.684 + b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
1.685 +#endif
1.686 +
1.687 + i = 9;
1.688 + if (9 < nd0) {
1.689 + s += 9;
1.690 + do b = multadd(PASS_STATE b, 10, *s++ - '0');
1.691 + while(++i < nd0);
1.692 + s++;
1.693 + }
1.694 + else
1.695 + s += 10;
1.696 + for(; i < nd; i++)
1.697 + b = multadd(PASS_STATE b, 10, *s++ - '0');
1.698 + return b;
1.699 + }
1.700 +
1.701 + static int
1.702 +hi0bits
1.703 +#ifdef KR_headers
1.704 + (x) ULong x;
1.705 +#else
1.706 + (ULong x)
1.707 +#endif
1.708 +{
1.709 + int k = 0;
1.710 +
1.711 + if (!(x & 0xffff0000)) {
1.712 + k = 16;
1.713 + x <<= 16;
1.714 + }
1.715 + if (!(x & 0xff000000)) {
1.716 + k += 8;
1.717 + x <<= 8;
1.718 + }
1.719 + if (!(x & 0xf0000000)) {
1.720 + k += 4;
1.721 + x <<= 4;
1.722 + }
1.723 + if (!(x & 0xc0000000)) {
1.724 + k += 2;
1.725 + x <<= 2;
1.726 + }
1.727 + if (!(x & 0x80000000)) {
1.728 + k++;
1.729 + if (!(x & 0x40000000))
1.730 + return 32;
1.731 + }
1.732 + return k;
1.733 + }
1.734 +
1.735 + static int
1.736 +lo0bits
1.737 +#ifdef KR_headers
1.738 + (y) ULong *y;
1.739 +#else
1.740 + (ULong *y)
1.741 +#endif
1.742 +{
1.743 + int k;
1.744 + ULong x = *y;
1.745 +
1.746 + if (x & 7) {
1.747 + if (x & 1)
1.748 + return 0;
1.749 + if (x & 2) {
1.750 + *y = x >> 1;
1.751 + return 1;
1.752 + }
1.753 + *y = x >> 2;
1.754 + return 2;
1.755 + }
1.756 + k = 0;
1.757 + if (!(x & 0xffff)) {
1.758 + k = 16;
1.759 + x >>= 16;
1.760 + }
1.761 + if (!(x & 0xff)) {
1.762 + k += 8;
1.763 + x >>= 8;
1.764 + }
1.765 + if (!(x & 0xf)) {
1.766 + k += 4;
1.767 + x >>= 4;
1.768 + }
1.769 + if (!(x & 0x3)) {
1.770 + k += 2;
1.771 + x >>= 2;
1.772 + }
1.773 + if (!(x & 1)) {
1.774 + k++;
1.775 + x >>= 1;
1.776 + if (!x)
1.777 + return 32;
1.778 + }
1.779 + *y = x;
1.780 + return k;
1.781 + }
1.782 +
1.783 + static Bigint *
1.784 +i2b
1.785 +#ifdef KR_headers
1.786 + (STATE_PARAM i) STATE_PARAM_DECL int i;
1.787 +#else
1.788 + (STATE_PARAM int i)
1.789 +#endif
1.790 +{
1.791 + Bigint *b;
1.792 +
1.793 + b = Balloc(PASS_STATE 1);
1.794 + b->x[0] = i;
1.795 + b->wds = 1;
1.796 + return b;
1.797 + }
1.798 +
1.799 + static Bigint *
1.800 +mult
1.801 +#ifdef KR_headers
1.802 + (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
1.803 +#else
1.804 + (STATE_PARAM Bigint *a, Bigint *b)
1.805 +#endif
1.806 +{
1.807 + Bigint *c;
1.808 + int k, wa, wb, wc;
1.809 + ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
1.810 + ULong y;
1.811 +#ifdef ULLong
1.812 + ULLong carry, z;
1.813 +#else
1.814 + ULong carry, z;
1.815 +#ifdef Pack_32
1.816 + ULong z2;
1.817 +#endif
1.818 +#endif
1.819 +
1.820 + if (a->wds < b->wds) {
1.821 + c = a;
1.822 + a = b;
1.823 + b = c;
1.824 + }
1.825 + k = a->k;
1.826 + wa = a->wds;
1.827 + wb = b->wds;
1.828 + wc = wa + wb;
1.829 + if (wc > a->maxwds)
1.830 + k++;
1.831 + c = Balloc(PASS_STATE k);
1.832 + for(x = c->x, xa = x + wc; x < xa; x++)
1.833 + *x = 0;
1.834 + xa = a->x;
1.835 + xae = xa + wa;
1.836 + xb = b->x;
1.837 + xbe = xb + wb;
1.838 + xc0 = c->x;
1.839 +#ifdef ULLong
1.840 + for(; xb < xbe; xc0++) {
1.841 + if ((y = *xb++)) {
1.842 + x = xa;
1.843 + xc = xc0;
1.844 + carry = 0;
1.845 + do {
1.846 + z = *x++ * (ULLong)y + *xc + carry;
1.847 + carry = z >> 32;
1.848 + *xc++ = (ULong) z & FFFFFFFF;
1.849 + }
1.850 + while(x < xae);
1.851 + *xc = (ULong) carry;
1.852 + }
1.853 + }
1.854 +#else
1.855 +#ifdef Pack_32
1.856 + for(; xb < xbe; xb++, xc0++) {
1.857 + if (y = *xb & 0xffff) {
1.858 + x = xa;
1.859 + xc = xc0;
1.860 + carry = 0;
1.861 + do {
1.862 + z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
1.863 + carry = z >> 16;
1.864 + z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
1.865 + carry = z2 >> 16;
1.866 + Storeinc(xc, z2, z);
1.867 + }
1.868 + while(x < xae);
1.869 + *xc = carry;
1.870 + }
1.871 + if (y = *xb >> 16) {
1.872 + x = xa;
1.873 + xc = xc0;
1.874 + carry = 0;
1.875 + z2 = *xc;
1.876 + do {
1.877 + z = (*x & 0xffff) * y + (*xc >> 16) + carry;
1.878 + carry = z >> 16;
1.879 + Storeinc(xc, z, z2);
1.880 + z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
1.881 + carry = z2 >> 16;
1.882 + }
1.883 + while(x < xae);
1.884 + *xc = z2;
1.885 + }
1.886 + }
1.887 +#else
1.888 + for(; xb < xbe; xc0++) {
1.889 + if (y = *xb++) {
1.890 + x = xa;
1.891 + xc = xc0;
1.892 + carry = 0;
1.893 + do {
1.894 + z = *x++ * y + *xc + carry;
1.895 + carry = z >> 16;
1.896 + *xc++ = z & 0xffff;
1.897 + }
1.898 + while(x < xae);
1.899 + *xc = carry;
1.900 + }
1.901 + }
1.902 +#endif
1.903 +#endif
1.904 + for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
1.905 + c->wds = wc;
1.906 + return c;
1.907 + }
1.908 +
1.909 + static Bigint *
1.910 +pow5mult
1.911 +#ifdef KR_headers
1.912 + (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
1.913 +#else
1.914 + (STATE_PARAM Bigint *b, int k)
1.915 +#endif
1.916 +{
1.917 + Bigint *b1, *p5, *p51;
1.918 + int i;
1.919 + static CONST int p05[3] = { 5, 25, 125 };
1.920 +
1.921 + if ((i = k & 3))
1.922 + b = multadd(PASS_STATE b, p05[i-1], 0);
1.923 +
1.924 + if (!(k >>= 2))
1.925 + return b;
1.926 + if (!(p5 = GET_STATE(p5s))) {
1.927 + /* first time */
1.928 +#ifdef MULTIPLE_THREADS
1.929 + ACQUIRE_DTOA_LOCK(1);
1.930 + if (!(p5 = p5s)) {
1.931 + p5 = p5s = i2b(625);
1.932 + p5->next = 0;
1.933 + }
1.934 + FREE_DTOA_LOCK(1);
1.935 +#else
1.936 + p5 = GET_STATE(p5s) = i2b(PASS_STATE 625);
1.937 + p5->next = 0;
1.938 +#endif
1.939 + }
1.940 + for(;;) {
1.941 + if (k & 1) {
1.942 + b1 = mult(PASS_STATE b, p5);
1.943 + Bfree(PASS_STATE b);
1.944 + b = b1;
1.945 + }
1.946 + if (!(k >>= 1))
1.947 + break;
1.948 + if (!(p51 = p5->next)) {
1.949 +#ifdef MULTIPLE_THREADS
1.950 + ACQUIRE_DTOA_LOCK(1);
1.951 + if (!(p51 = p5->next)) {
1.952 + p51 = p5->next = mult(p5,p5);
1.953 + p51->next = 0;
1.954 + }
1.955 + FREE_DTOA_LOCK(1);
1.956 +#else
1.957 + p51 = p5->next = mult(PASS_STATE p5,p5);
1.958 + p51->next = 0;
1.959 +#endif
1.960 + }
1.961 + p5 = p51;
1.962 + }
1.963 + return b;
1.964 + }
1.965 +
1.966 + static Bigint *
1.967 +lshift
1.968 +#ifdef KR_headers
1.969 + (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
1.970 +#else
1.971 + (STATE_PARAM Bigint *b, int k)
1.972 +#endif
1.973 +{
1.974 + int i, k1, n, n1;
1.975 + Bigint *b1;
1.976 + ULong *x, *x1, *xe, z;
1.977 +
1.978 +#ifdef Pack_32
1.979 + n = k >> 5;
1.980 +#else
1.981 + n = k >> 4;
1.982 +#endif
1.983 + k1 = b->k;
1.984 + n1 = n + b->wds + 1;
1.985 + for(i = b->maxwds; n1 > i; i <<= 1)
1.986 + k1++;
1.987 + b1 = Balloc(PASS_STATE k1);
1.988 + x1 = b1->x;
1.989 + for(i = 0; i < n; i++)
1.990 + *x1++ = 0;
1.991 + x = b->x;
1.992 + xe = x + b->wds;
1.993 +#ifdef Pack_32
1.994 + if (k &= 0x1f) {
1.995 + k1 = 32 - k;
1.996 + z = 0;
1.997 + do {
1.998 + *x1++ = *x << k | z;
1.999 + z = *x++ >> k1;
1.1000 + }
1.1001 + while(x < xe);
1.1002 + if ((*x1 = z))
1.1003 + ++n1;
1.1004 + }
1.1005 +#else
1.1006 + if (k &= 0xf) {
1.1007 + k1 = 16 - k;
1.1008 + z = 0;
1.1009 + do {
1.1010 + *x1++ = *x << k & 0xffff | z;
1.1011 + z = *x++ >> k1;
1.1012 + }
1.1013 + while(x < xe);
1.1014 + if (*x1 = z)
1.1015 + ++n1;
1.1016 + }
1.1017 +#endif
1.1018 + else do
1.1019 + *x1++ = *x++;
1.1020 + while(x < xe);
1.1021 + b1->wds = n1 - 1;
1.1022 + Bfree(PASS_STATE b);
1.1023 + return b1;
1.1024 + }
1.1025 +
1.1026 + static int
1.1027 +cmp
1.1028 +#ifdef KR_headers
1.1029 + (a, b) Bigint *a, *b;
1.1030 +#else
1.1031 + (Bigint *a, Bigint *b)
1.1032 +#endif
1.1033 +{
1.1034 + ULong *xa, *xa0, *xb, *xb0;
1.1035 + int i, j;
1.1036 +
1.1037 + i = a->wds;
1.1038 + j = b->wds;
1.1039 +#ifdef DEBUG
1.1040 + if (i > 1 && !a->x[i-1])
1.1041 + Bug("cmp called with a->x[a->wds-1] == 0");
1.1042 + if (j > 1 && !b->x[j-1])
1.1043 + Bug("cmp called with b->x[b->wds-1] == 0");
1.1044 +#endif
1.1045 + if (i -= j)
1.1046 + return i;
1.1047 + xa0 = a->x;
1.1048 + xa = xa0 + j;
1.1049 + xb0 = b->x;
1.1050 + xb = xb0 + j;
1.1051 + for(;;) {
1.1052 + if (*--xa != *--xb)
1.1053 + return *xa < *xb ? -1 : 1;
1.1054 + if (xa <= xa0)
1.1055 + break;
1.1056 + }
1.1057 + return 0;
1.1058 + }
1.1059 +
1.1060 + static Bigint *
1.1061 +diff
1.1062 +#ifdef KR_headers
1.1063 + (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
1.1064 +#else
1.1065 + (STATE_PARAM Bigint *a, Bigint *b)
1.1066 +#endif
1.1067 +{
1.1068 + Bigint *c;
1.1069 + int i, wa, wb;
1.1070 + ULong *xa, *xae, *xb, *xbe, *xc;
1.1071 +#ifdef ULLong
1.1072 + ULLong borrow, y;
1.1073 +#else
1.1074 + ULong borrow, y;
1.1075 +#ifdef Pack_32
1.1076 + ULong z;
1.1077 +#endif
1.1078 +#endif
1.1079 +
1.1080 + i = cmp(a,b);
1.1081 + if (!i) {
1.1082 + c = Balloc(PASS_STATE 0);
1.1083 + c->wds = 1;
1.1084 + c->x[0] = 0;
1.1085 + return c;
1.1086 + }
1.1087 + if (i < 0) {
1.1088 + c = a;
1.1089 + a = b;
1.1090 + b = c;
1.1091 + i = 1;
1.1092 + }
1.1093 + else
1.1094 + i = 0;
1.1095 + c = Balloc(PASS_STATE a->k);
1.1096 + c->sign = i;
1.1097 + wa = a->wds;
1.1098 + xa = a->x;
1.1099 + xae = xa + wa;
1.1100 + wb = b->wds;
1.1101 + xb = b->x;
1.1102 + xbe = xb + wb;
1.1103 + xc = c->x;
1.1104 + borrow = 0;
1.1105 +#ifdef ULLong
1.1106 + do {
1.1107 + y = (ULLong)*xa++ - *xb++ - borrow;
1.1108 + borrow = y >> 32 & (ULong)1;
1.1109 + *xc++ = (ULong) y & FFFFFFFF;
1.1110 + }
1.1111 + while(xb < xbe);
1.1112 + while(xa < xae) {
1.1113 + y = *xa++ - borrow;
1.1114 + borrow = y >> 32 & (ULong)1;
1.1115 + *xc++ = (ULong) y & FFFFFFFF;
1.1116 + }
1.1117 +#else
1.1118 +#ifdef Pack_32
1.1119 + do {
1.1120 + y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1.1121 + borrow = (y & 0x10000) >> 16;
1.1122 + z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1.1123 + borrow = (z & 0x10000) >> 16;
1.1124 + Storeinc(xc, z, y);
1.1125 + }
1.1126 + while(xb < xbe);
1.1127 + while(xa < xae) {
1.1128 + y = (*xa & 0xffff) - borrow;
1.1129 + borrow = (y & 0x10000) >> 16;
1.1130 + z = (*xa++ >> 16) - borrow;
1.1131 + borrow = (z & 0x10000) >> 16;
1.1132 + Storeinc(xc, z, y);
1.1133 + }
1.1134 +#else
1.1135 + do {
1.1136 + y = *xa++ - *xb++ - borrow;
1.1137 + borrow = (y & 0x10000) >> 16;
1.1138 + *xc++ = y & 0xffff;
1.1139 + }
1.1140 + while(xb < xbe);
1.1141 + while(xa < xae) {
1.1142 + y = *xa++ - borrow;
1.1143 + borrow = (y & 0x10000) >> 16;
1.1144 + *xc++ = y & 0xffff;
1.1145 + }
1.1146 +#endif
1.1147 +#endif
1.1148 + while(!*--xc)
1.1149 + wa--;
1.1150 + c->wds = wa;
1.1151 + return c;
1.1152 + }
1.1153 +
1.1154 + static double
1.1155 +ulp
1.1156 +#ifdef KR_headers
1.1157 + (x) U x;
1.1158 +#else
1.1159 + (U x)
1.1160 +#endif
1.1161 +{
1.1162 + Long L;
1.1163 + U a;
1.1164 +
1.1165 + L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1.1166 +#ifndef Avoid_Underflow
1.1167 +#ifndef Sudden_Underflow
1.1168 + if (L > 0) {
1.1169 +#endif
1.1170 +#endif
1.1171 +#ifdef IBM
1.1172 + L |= Exp_msk1 >> 4;
1.1173 +#endif
1.1174 + word0(a) = L;
1.1175 + word1(a) = 0;
1.1176 +#ifndef Avoid_Underflow
1.1177 +#ifndef Sudden_Underflow
1.1178 + }
1.1179 + else {
1.1180 + L = -L >> Exp_shift;
1.1181 + if (L < Exp_shift) {
1.1182 + word0(a) = 0x80000 >> L;
1.1183 + word1(a) = 0;
1.1184 + }
1.1185 + else {
1.1186 + word0(a) = 0;
1.1187 + L -= Exp_shift;
1.1188 + word1(a) = L >= 31 ? 1 : 1 << 31 - L;
1.1189 + }
1.1190 + }
1.1191 +#endif
1.1192 +#endif
1.1193 + return dval(a);
1.1194 + }
1.1195 +
1.1196 + static double
1.1197 +b2d
1.1198 +#ifdef KR_headers
1.1199 + (a, e) Bigint *a; int *e;
1.1200 +#else
1.1201 + (Bigint *a, int *e)
1.1202 +#endif
1.1203 +{
1.1204 + ULong *xa, *xa0, w, y, z;
1.1205 + int k;
1.1206 + U d;
1.1207 +#ifdef VAX
1.1208 + ULong d0, d1;
1.1209 +#else
1.1210 +#define d0 word0(d)
1.1211 +#define d1 word1(d)
1.1212 +#endif
1.1213 +
1.1214 + xa0 = a->x;
1.1215 + xa = xa0 + a->wds;
1.1216 + y = *--xa;
1.1217 +#ifdef DEBUG
1.1218 + if (!y) Bug("zero y in b2d");
1.1219 +#endif
1.1220 + k = hi0bits(y);
1.1221 + *e = 32 - k;
1.1222 +#ifdef Pack_32
1.1223 + if (k < Ebits) {
1.1224 + d0 = Exp_1 | y >> (Ebits - k);
1.1225 + w = xa > xa0 ? *--xa : 0;
1.1226 + d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1.1227 + goto ret_d;
1.1228 + }
1.1229 + z = xa > xa0 ? *--xa : 0;
1.1230 + if (k -= Ebits) {
1.1231 + d0 = Exp_1 | y << k | z >> (32 - k);
1.1232 + y = xa > xa0 ? *--xa : 0;
1.1233 + d1 = z << k | y >> (32 - k);
1.1234 + }
1.1235 + else {
1.1236 + d0 = Exp_1 | y;
1.1237 + d1 = z;
1.1238 + }
1.1239 +#else
1.1240 + if (k < Ebits + 16) {
1.1241 + z = xa > xa0 ? *--xa : 0;
1.1242 + d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1.1243 + w = xa > xa0 ? *--xa : 0;
1.1244 + y = xa > xa0 ? *--xa : 0;
1.1245 + d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1.1246 + goto ret_d;
1.1247 + }
1.1248 + z = xa > xa0 ? *--xa : 0;
1.1249 + w = xa > xa0 ? *--xa : 0;
1.1250 + k -= Ebits + 16;
1.1251 + d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1.1252 + y = xa > xa0 ? *--xa : 0;
1.1253 + d1 = w << k + 16 | y << k;
1.1254 +#endif
1.1255 + ret_d:
1.1256 +#ifdef VAX
1.1257 + word0(d) = d0 >> 16 | d0 << 16;
1.1258 + word1(d) = d1 >> 16 | d1 << 16;
1.1259 +#else
1.1260 +#undef d0
1.1261 +#undef d1
1.1262 +#endif
1.1263 + return dval(d);
1.1264 + }
1.1265 +
1.1266 + static Bigint *
1.1267 +d2b
1.1268 +#ifdef KR_headers
1.1269 + (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits;
1.1270 +#else
1.1271 + (STATE_PARAM U d, int *e, int *bits)
1.1272 +#endif
1.1273 +{
1.1274 + Bigint *b;
1.1275 + int de, k;
1.1276 + ULong *x, y, z;
1.1277 +#ifndef Sudden_Underflow
1.1278 + int i;
1.1279 +#endif
1.1280 +#ifdef VAX
1.1281 + ULong d0, d1;
1.1282 + d0 = word0(d) >> 16 | word0(d) << 16;
1.1283 + d1 = word1(d) >> 16 | word1(d) << 16;
1.1284 +#else
1.1285 +#define d0 word0(d)
1.1286 +#define d1 word1(d)
1.1287 +#endif
1.1288 +
1.1289 +#ifdef Pack_32
1.1290 + b = Balloc(PASS_STATE 1);
1.1291 +#else
1.1292 + b = Balloc(PASS_STATE 2);
1.1293 +#endif
1.1294 + x = b->x;
1.1295 +
1.1296 + z = d0 & Frac_mask;
1.1297 + d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1.1298 +#ifdef Sudden_Underflow
1.1299 + de = (int)(d0 >> Exp_shift);
1.1300 +#ifndef IBM
1.1301 + z |= Exp_msk11;
1.1302 +#endif
1.1303 +#else
1.1304 + if ((de = (int)(d0 >> Exp_shift)))
1.1305 + z |= Exp_msk1;
1.1306 +#endif
1.1307 +#ifdef Pack_32
1.1308 + if ((y = d1)) {
1.1309 + if ((k = lo0bits(&y))) {
1.1310 + x[0] = y | z << (32 - k);
1.1311 + z >>= k;
1.1312 + }
1.1313 + else
1.1314 + x[0] = y;
1.1315 +#ifndef Sudden_Underflow
1.1316 + i =
1.1317 +#endif
1.1318 + b->wds = (x[1] = z) ? 2 : 1;
1.1319 + }
1.1320 + else {
1.1321 + k = lo0bits(&z);
1.1322 + x[0] = z;
1.1323 +#ifndef Sudden_Underflow
1.1324 + i =
1.1325 +#endif
1.1326 + b->wds = 1;
1.1327 + k += 32;
1.1328 + }
1.1329 +#else
1.1330 + if (y = d1) {
1.1331 + if (k = lo0bits(&y))
1.1332 + if (k >= 16) {
1.1333 + x[0] = y | z << 32 - k & 0xffff;
1.1334 + x[1] = z >> k - 16 & 0xffff;
1.1335 + x[2] = z >> k;
1.1336 + i = 2;
1.1337 + }
1.1338 + else {
1.1339 + x[0] = y & 0xffff;
1.1340 + x[1] = y >> 16 | z << 16 - k & 0xffff;
1.1341 + x[2] = z >> k & 0xffff;
1.1342 + x[3] = z >> k+16;
1.1343 + i = 3;
1.1344 + }
1.1345 + else {
1.1346 + x[0] = y & 0xffff;
1.1347 + x[1] = y >> 16;
1.1348 + x[2] = z & 0xffff;
1.1349 + x[3] = z >> 16;
1.1350 + i = 3;
1.1351 + }
1.1352 + }
1.1353 + else {
1.1354 +#ifdef DEBUG
1.1355 + if (!z)
1.1356 + Bug("Zero passed to d2b");
1.1357 +#endif
1.1358 + k = lo0bits(&z);
1.1359 + if (k >= 16) {
1.1360 + x[0] = z;
1.1361 + i = 0;
1.1362 + }
1.1363 + else {
1.1364 + x[0] = z & 0xffff;
1.1365 + x[1] = z >> 16;
1.1366 + i = 1;
1.1367 + }
1.1368 + k += 32;
1.1369 + }
1.1370 + while(!x[i])
1.1371 + --i;
1.1372 + b->wds = i + 1;
1.1373 +#endif
1.1374 +#ifndef Sudden_Underflow
1.1375 + if (de) {
1.1376 +#endif
1.1377 +#ifdef IBM
1.1378 + *e = (de - Bias - (P-1) << 2) + k;
1.1379 + *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1.1380 +#else
1.1381 + *e = de - Bias - (P-1) + k;
1.1382 + *bits = P - k;
1.1383 +#endif
1.1384 +#ifndef Sudden_Underflow
1.1385 + }
1.1386 + else {
1.1387 + *e = de - Bias - (P-1) + 1 + k;
1.1388 +#ifdef Pack_32
1.1389 + *bits = 32*i - hi0bits(x[i-1]);
1.1390 +#else
1.1391 + *bits = (i+2)*16 - hi0bits(x[i]);
1.1392 +#endif
1.1393 + }
1.1394 +#endif
1.1395 + return b;
1.1396 + }
1.1397 +#undef d0
1.1398 +#undef d1
1.1399 +
1.1400 + static double
1.1401 +ratio
1.1402 +#ifdef KR_headers
1.1403 + (a, b) Bigint *a, *b;
1.1404 +#else
1.1405 + (Bigint *a, Bigint *b)
1.1406 +#endif
1.1407 +{
1.1408 + U da, db;
1.1409 + int k, ka, kb;
1.1410 +
1.1411 + dval(da) = b2d(a, &ka);
1.1412 + dval(db) = b2d(b, &kb);
1.1413 +#ifdef Pack_32
1.1414 + k = ka - kb + 32*(a->wds - b->wds);
1.1415 +#else
1.1416 + k = ka - kb + 16*(a->wds - b->wds);
1.1417 +#endif
1.1418 +#ifdef IBM
1.1419 + if (k > 0) {
1.1420 + word0(da) += (k >> 2)*Exp_msk1;
1.1421 + if (k &= 3)
1.1422 + dval(da) *= 1 << k;
1.1423 + }
1.1424 + else {
1.1425 + k = -k;
1.1426 + word0(db) += (k >> 2)*Exp_msk1;
1.1427 + if (k &= 3)
1.1428 + dval(db) *= 1 << k;
1.1429 + }
1.1430 +#else
1.1431 + if (k > 0)
1.1432 + word0(da) += k*Exp_msk1;
1.1433 + else {
1.1434 + k = -k;
1.1435 + word0(db) += k*Exp_msk1;
1.1436 + }
1.1437 +#endif
1.1438 + return dval(da) / dval(db);
1.1439 + }
1.1440 +
1.1441 + static CONST double
1.1442 +tens[] = {
1.1443 + 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1.1444 + 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1.1445 + 1e20, 1e21, 1e22
1.1446 +#ifdef VAX
1.1447 + , 1e23, 1e24
1.1448 +#endif
1.1449 + };
1.1450 +
1.1451 + static CONST double
1.1452 +#ifdef IEEE_Arith
1.1453 +bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1.1454 +static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1.1455 +#ifdef Avoid_Underflow
1.1456 + 9007199254740992.*9007199254740992.e-256
1.1457 + /* = 2^106 * 1e-53 */
1.1458 +#else
1.1459 + 1e-256
1.1460 +#endif
1.1461 + };
1.1462 +/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1.1463 +/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1.1464 +#define Scale_Bit 0x10
1.1465 +#define n_bigtens 5
1.1466 +#else
1.1467 +#ifdef IBM
1.1468 +bigtens[] = { 1e16, 1e32, 1e64 };
1.1469 +static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1.1470 +#define n_bigtens 3
1.1471 +#else
1.1472 +bigtens[] = { 1e16, 1e32 };
1.1473 +static CONST double tinytens[] = { 1e-16, 1e-32 };
1.1474 +#define n_bigtens 2
1.1475 +#endif
1.1476 +#endif
1.1477 +
1.1478 + static double
1.1479 +_strtod
1.1480 +#ifdef KR_headers
1.1481 + (STATE_PARAM s00, se) STATE_PARAM_DECL CONST char *s00; char **se;
1.1482 +#else
1.1483 + (STATE_PARAM CONST char *s00, char **se)
1.1484 +#endif
1.1485 +{
1.1486 +#ifdef Avoid_Underflow
1.1487 + int scale;
1.1488 +#endif
1.1489 + int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1.1490 + e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1.1491 + CONST char *s, *s0, *s1;
1.1492 + double aadj, adj;
1.1493 + U aadj1, rv, rv0;
1.1494 + Long L;
1.1495 + ULong y, z;
1.1496 + Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1.1497 +#ifdef SET_INEXACT
1.1498 + int inexact, oldinexact;
1.1499 +#endif
1.1500 +#ifdef Honor_FLT_ROUNDS
1.1501 + int rounding;
1.1502 +#endif
1.1503 +#ifdef USE_LOCALE
1.1504 + CONST char *s2;
1.1505 +#endif
1.1506 +
1.1507 +#ifdef __GNUC__
1.1508 + delta = bb = bd = bs = 0;
1.1509 +#endif
1.1510 +
1.1511 + sign = nz0 = nz = 0;
1.1512 + dval(rv) = 0.;
1.1513 + for(s = s00;;s++) switch(*s) {
1.1514 + case '-':
1.1515 + sign = 1;
1.1516 + /* no break */
1.1517 + case '+':
1.1518 + if (*++s)
1.1519 + goto break2;
1.1520 + /* no break */
1.1521 + case 0:
1.1522 + goto ret0;
1.1523 + case '\t':
1.1524 + case '\n':
1.1525 + case '\v':
1.1526 + case '\f':
1.1527 + case '\r':
1.1528 + case ' ':
1.1529 + continue;
1.1530 + default:
1.1531 + goto break2;
1.1532 + }
1.1533 + break2:
1.1534 + if (*s == '0') {
1.1535 + nz0 = 1;
1.1536 + while(*++s == '0') ;
1.1537 + if (!*s)
1.1538 + goto ret;
1.1539 + }
1.1540 + s0 = s;
1.1541 + y = z = 0;
1.1542 + for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1.1543 + if (nd < 9)
1.1544 + y = 10*y + c - '0';
1.1545 + else if (nd < 16)
1.1546 + z = 10*z + c - '0';
1.1547 + nd0 = nd;
1.1548 +#ifdef USE_LOCALE
1.1549 + s1 = localeconv()->decimal_point;
1.1550 + if (c == *s1) {
1.1551 + c = '.';
1.1552 + if (*++s1) {
1.1553 + s2 = s;
1.1554 + for(;;) {
1.1555 + if (*++s2 != *s1) {
1.1556 + c = 0;
1.1557 + break;
1.1558 + }
1.1559 + if (!*++s1) {
1.1560 + s = s2;
1.1561 + break;
1.1562 + }
1.1563 + }
1.1564 + }
1.1565 + }
1.1566 +#endif
1.1567 + if (c == '.') {
1.1568 + c = *++s;
1.1569 + if (!nd) {
1.1570 + for(; c == '0'; c = *++s)
1.1571 + nz++;
1.1572 + if (c > '0' && c <= '9') {
1.1573 + s0 = s;
1.1574 + nf += nz;
1.1575 + nz = 0;
1.1576 + goto have_dig;
1.1577 + }
1.1578 + goto dig_done;
1.1579 + }
1.1580 + for(; c >= '0' && c <= '9'; c = *++s) {
1.1581 + have_dig:
1.1582 + nz++;
1.1583 + if (c -= '0') {
1.1584 + nf += nz;
1.1585 + for(i = 1; i < nz; i++)
1.1586 + if (nd++ < 9)
1.1587 + y *= 10;
1.1588 + else if (nd <= DBL_DIG + 1)
1.1589 + z *= 10;
1.1590 + if (nd++ < 9)
1.1591 + y = 10*y + c;
1.1592 + else if (nd <= DBL_DIG + 1)
1.1593 + z = 10*z + c;
1.1594 + nz = 0;
1.1595 + }
1.1596 + }
1.1597 + }
1.1598 + dig_done:
1.1599 + e = 0;
1.1600 + if (c == 'e' || c == 'E') {
1.1601 + if (!nd && !nz && !nz0) {
1.1602 + goto ret0;
1.1603 + }
1.1604 + s00 = s;
1.1605 + esign = 0;
1.1606 + switch(c = *++s) {
1.1607 + case '-':
1.1608 + esign = 1;
1.1609 + case '+':
1.1610 + c = *++s;
1.1611 + }
1.1612 + if (c >= '0' && c <= '9') {
1.1613 + while(c == '0')
1.1614 + c = *++s;
1.1615 + if (c > '0' && c <= '9') {
1.1616 + L = c - '0';
1.1617 + s1 = s;
1.1618 + while((c = *++s) >= '0' && c <= '9')
1.1619 + L = 10*L + c - '0';
1.1620 + if (s - s1 > 8 || L > 19999)
1.1621 + /* Avoid confusion from exponents
1.1622 + * so large that e might overflow.
1.1623 + */
1.1624 + e = 19999; /* safe for 16 bit ints */
1.1625 + else
1.1626 + e = (int)L;
1.1627 + if (esign)
1.1628 + e = -e;
1.1629 + }
1.1630 + else
1.1631 + e = 0;
1.1632 + }
1.1633 + else
1.1634 + s = s00;
1.1635 + }
1.1636 + if (!nd) {
1.1637 + if (!nz && !nz0) {
1.1638 + ret0:
1.1639 + s = s00;
1.1640 + sign = 0;
1.1641 + }
1.1642 + goto ret;
1.1643 + }
1.1644 + e1 = e -= nf;
1.1645 +
1.1646 + /* Now we have nd0 digits, starting at s0, followed by a
1.1647 + * decimal point, followed by nd-nd0 digits. The number we're
1.1648 + * after is the integer represented by those digits times
1.1649 + * 10**e */
1.1650 +
1.1651 + if (!nd0)
1.1652 + nd0 = nd;
1.1653 + k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1.1654 + dval(rv) = y;
1.1655 + if (k > 9) {
1.1656 +#ifdef SET_INEXACT
1.1657 + if (k > DBL_DIG)
1.1658 + oldinexact = get_inexact();
1.1659 +#endif
1.1660 + dval(rv) = tens[k - 9] * dval(rv) + z;
1.1661 + }
1.1662 + bd0 = 0;
1.1663 + if (nd <= DBL_DIG
1.1664 +#ifndef RND_PRODQUOT
1.1665 +#ifndef Honor_FLT_ROUNDS
1.1666 + && Flt_Rounds == 1
1.1667 +#endif
1.1668 +#endif
1.1669 + ) {
1.1670 + if (!e)
1.1671 + goto ret;
1.1672 + if (e > 0) {
1.1673 + if (e <= Ten_pmax) {
1.1674 +#ifdef VAX
1.1675 + goto vax_ovfl_check;
1.1676 +#else
1.1677 +#ifdef Honor_FLT_ROUNDS
1.1678 + /* round correctly FLT_ROUNDS = 2 or 3 */
1.1679 + if (sign) {
1.1680 + rv = -rv;
1.1681 + sign = 0;
1.1682 + }
1.1683 +#endif
1.1684 + /* rv = */ rounded_product(dval(rv), tens[e]);
1.1685 + goto ret;
1.1686 +#endif
1.1687 + }
1.1688 + i = DBL_DIG - nd;
1.1689 + if (e <= Ten_pmax + i) {
1.1690 + /* A fancier test would sometimes let us do
1.1691 + * this for larger i values.
1.1692 + */
1.1693 +#ifdef Honor_FLT_ROUNDS
1.1694 + /* round correctly FLT_ROUNDS = 2 or 3 */
1.1695 + if (sign) {
1.1696 + rv = -rv;
1.1697 + sign = 0;
1.1698 + }
1.1699 +#endif
1.1700 + e -= i;
1.1701 + dval(rv) *= tens[i];
1.1702 +#ifdef VAX
1.1703 + /* VAX exponent range is so narrow we must
1.1704 + * worry about overflow here...
1.1705 + */
1.1706 + vax_ovfl_check:
1.1707 + word0(rv) -= P*Exp_msk1;
1.1708 + /* rv = */ rounded_product(dval(rv), tens[e]);
1.1709 + if ((word0(rv) & Exp_mask)
1.1710 + > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1.1711 + goto ovfl;
1.1712 + word0(rv) += P*Exp_msk1;
1.1713 +#else
1.1714 + /* rv = */ rounded_product(dval(rv), tens[e]);
1.1715 +#endif
1.1716 + goto ret;
1.1717 + }
1.1718 + }
1.1719 +#ifndef Inaccurate_Divide
1.1720 + else if (e >= -Ten_pmax) {
1.1721 +#ifdef Honor_FLT_ROUNDS
1.1722 + /* round correctly FLT_ROUNDS = 2 or 3 */
1.1723 + if (sign) {
1.1724 + rv = -rv;
1.1725 + sign = 0;
1.1726 + }
1.1727 +#endif
1.1728 + /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1.1729 + goto ret;
1.1730 + }
1.1731 +#endif
1.1732 + }
1.1733 + e1 += nd - k;
1.1734 +
1.1735 +#ifdef IEEE_Arith
1.1736 +#ifdef SET_INEXACT
1.1737 + inexact = 1;
1.1738 + if (k <= DBL_DIG)
1.1739 + oldinexact = get_inexact();
1.1740 +#endif
1.1741 +#ifdef Avoid_Underflow
1.1742 + scale = 0;
1.1743 +#endif
1.1744 +#ifdef Honor_FLT_ROUNDS
1.1745 + if ((rounding = Flt_Rounds) >= 2) {
1.1746 + if (sign)
1.1747 + rounding = rounding == 2 ? 0 : 2;
1.1748 + else
1.1749 + if (rounding != 2)
1.1750 + rounding = 0;
1.1751 + }
1.1752 +#endif
1.1753 +#endif /*IEEE_Arith*/
1.1754 +
1.1755 + /* Get starting approximation = rv * 10**e1 */
1.1756 +
1.1757 + if (e1 > 0) {
1.1758 + if ((i = e1 & 15))
1.1759 + dval(rv) *= tens[i];
1.1760 + if (e1 &= ~15) {
1.1761 + if (e1 > DBL_MAX_10_EXP) {
1.1762 + ovfl:
1.1763 +#ifndef NO_ERRNO
1.1764 + errno = ERANGE;
1.1765 +#endif
1.1766 + /* Can't trust HUGE_VAL */
1.1767 +#ifdef IEEE_Arith
1.1768 +#ifdef Honor_FLT_ROUNDS
1.1769 + switch(rounding) {
1.1770 + case 0: /* toward 0 */
1.1771 + case 3: /* toward -infinity */
1.1772 + word0(rv) = Big0;
1.1773 + word1(rv) = Big1;
1.1774 + break;
1.1775 + default:
1.1776 + word0(rv) = Exp_mask;
1.1777 + word1(rv) = 0;
1.1778 + }
1.1779 +#else /*Honor_FLT_ROUNDS*/
1.1780 + word0(rv) = Exp_mask;
1.1781 + word1(rv) = 0;
1.1782 +#endif /*Honor_FLT_ROUNDS*/
1.1783 +#ifdef SET_INEXACT
1.1784 + /* set overflow bit */
1.1785 + dval(rv0) = 1e300;
1.1786 + dval(rv0) *= dval(rv0);
1.1787 +#endif
1.1788 +#else /*IEEE_Arith*/
1.1789 + word0(rv) = Big0;
1.1790 + word1(rv) = Big1;
1.1791 +#endif /*IEEE_Arith*/
1.1792 + if (bd0)
1.1793 + goto retfree;
1.1794 + goto ret;
1.1795 + }
1.1796 + e1 >>= 4;
1.1797 + for(j = 0; e1 > 1; j++, e1 >>= 1)
1.1798 + if (e1 & 1)
1.1799 + dval(rv) *= bigtens[j];
1.1800 + /* The last multiplication could overflow. */
1.1801 + word0(rv) -= P*Exp_msk1;
1.1802 + dval(rv) *= bigtens[j];
1.1803 + if ((z = word0(rv) & Exp_mask)
1.1804 + > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1.1805 + goto ovfl;
1.1806 + if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1.1807 + /* set to largest number */
1.1808 + /* (Can't trust DBL_MAX) */
1.1809 + word0(rv) = Big0;
1.1810 + word1(rv) = Big1;
1.1811 + }
1.1812 + else
1.1813 + word0(rv) += P*Exp_msk1;
1.1814 + }
1.1815 + }
1.1816 + else if (e1 < 0) {
1.1817 + e1 = -e1;
1.1818 + if ((i = e1 & 15))
1.1819 + dval(rv) /= tens[i];
1.1820 + if (e1 >>= 4) {
1.1821 + if (e1 >= 1 << n_bigtens)
1.1822 + goto undfl;
1.1823 +#ifdef Avoid_Underflow
1.1824 + if (e1 & Scale_Bit)
1.1825 + scale = 2*P;
1.1826 + for(j = 0; e1 > 0; j++, e1 >>= 1)
1.1827 + if (e1 & 1)
1.1828 + dval(rv) *= tinytens[j];
1.1829 + if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1.1830 + >> Exp_shift)) > 0) {
1.1831 + /* scaled rv is denormal; zap j low bits */
1.1832 + if (j >= 32) {
1.1833 + word1(rv) = 0;
1.1834 + if (j >= 53)
1.1835 + word0(rv) = (P+2)*Exp_msk1;
1.1836 + else
1.1837 + word0(rv) &= 0xffffffff << (j-32);
1.1838 + }
1.1839 + else
1.1840 + word1(rv) &= 0xffffffff << j;
1.1841 + }
1.1842 +#else
1.1843 + for(j = 0; e1 > 1; j++, e1 >>= 1)
1.1844 + if (e1 & 1)
1.1845 + dval(rv) *= tinytens[j];
1.1846 + /* The last multiplication could underflow. */
1.1847 + dval(rv0) = dval(rv);
1.1848 + dval(rv) *= tinytens[j];
1.1849 + if (!dval(rv)) {
1.1850 + dval(rv) = 2.*dval(rv0);
1.1851 + dval(rv) *= tinytens[j];
1.1852 +#endif
1.1853 + if (!dval(rv)) {
1.1854 + undfl:
1.1855 + dval(rv) = 0.;
1.1856 +#ifndef NO_ERRNO
1.1857 + errno = ERANGE;
1.1858 +#endif
1.1859 + if (bd0)
1.1860 + goto retfree;
1.1861 + goto ret;
1.1862 + }
1.1863 +#ifndef Avoid_Underflow
1.1864 + word0(rv) = Tiny0;
1.1865 + word1(rv) = Tiny1;
1.1866 + /* The refinement below will clean
1.1867 + * this approximation up.
1.1868 + */
1.1869 + }
1.1870 +#endif
1.1871 + }
1.1872 + }
1.1873 +
1.1874 + /* Now the hard part -- adjusting rv to the correct value.*/
1.1875 +
1.1876 + /* Put digits into bd: true value = bd * 10^e */
1.1877 +
1.1878 + bd0 = s2b(PASS_STATE s0, nd0, nd, y);
1.1879 +
1.1880 + for(;;) {
1.1881 + bd = Balloc(PASS_STATE bd0->k);
1.1882 + Bcopy(bd, bd0);
1.1883 + bb = d2b(PASS_STATE rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1.1884 + bs = i2b(PASS_STATE 1);
1.1885 +
1.1886 + if (e >= 0) {
1.1887 + bb2 = bb5 = 0;
1.1888 + bd2 = bd5 = e;
1.1889 + }
1.1890 + else {
1.1891 + bb2 = bb5 = -e;
1.1892 + bd2 = bd5 = 0;
1.1893 + }
1.1894 + if (bbe >= 0)
1.1895 + bb2 += bbe;
1.1896 + else
1.1897 + bd2 -= bbe;
1.1898 + bs2 = bb2;
1.1899 +#ifdef Honor_FLT_ROUNDS
1.1900 + if (rounding != 1)
1.1901 + bs2++;
1.1902 +#endif
1.1903 +#ifdef Avoid_Underflow
1.1904 + j = bbe - scale;
1.1905 + i = j + bbbits - 1; /* logb(rv) */
1.1906 + if (i < Emin) /* denormal */
1.1907 + j += P - Emin;
1.1908 + else
1.1909 + j = P + 1 - bbbits;
1.1910 +#else /*Avoid_Underflow*/
1.1911 +#ifdef Sudden_Underflow
1.1912 +#ifdef IBM
1.1913 + j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1.1914 +#else
1.1915 + j = P + 1 - bbbits;
1.1916 +#endif
1.1917 +#else /*Sudden_Underflow*/
1.1918 + j = bbe;
1.1919 + i = j + bbbits - 1; /* logb(rv) */
1.1920 + if (i < Emin) /* denormal */
1.1921 + j += P - Emin;
1.1922 + else
1.1923 + j = P + 1 - bbbits;
1.1924 +#endif /*Sudden_Underflow*/
1.1925 +#endif /*Avoid_Underflow*/
1.1926 + bb2 += j;
1.1927 + bd2 += j;
1.1928 +#ifdef Avoid_Underflow
1.1929 + bd2 += scale;
1.1930 +#endif
1.1931 + i = bb2 < bd2 ? bb2 : bd2;
1.1932 + if (i > bs2)
1.1933 + i = bs2;
1.1934 + if (i > 0) {
1.1935 + bb2 -= i;
1.1936 + bd2 -= i;
1.1937 + bs2 -= i;
1.1938 + }
1.1939 + if (bb5 > 0) {
1.1940 + bs = pow5mult(PASS_STATE bs, bb5);
1.1941 + bb1 = mult(PASS_STATE bs, bb);
1.1942 + Bfree(PASS_STATE bb);
1.1943 + bb = bb1;
1.1944 + }
1.1945 + if (bb2 > 0)
1.1946 + bb = lshift(PASS_STATE bb, bb2);
1.1947 + if (bd5 > 0)
1.1948 + bd = pow5mult(PASS_STATE bd, bd5);
1.1949 + if (bd2 > 0)
1.1950 + bd = lshift(PASS_STATE bd, bd2);
1.1951 + if (bs2 > 0)
1.1952 + bs = lshift(PASS_STATE bs, bs2);
1.1953 + delta = diff(PASS_STATE bb, bd);
1.1954 + dsign = delta->sign;
1.1955 + delta->sign = 0;
1.1956 + i = cmp(delta, bs);
1.1957 +#ifdef Honor_FLT_ROUNDS
1.1958 + if (rounding != 1) {
1.1959 + if (i < 0) {
1.1960 + /* Error is less than an ulp */
1.1961 + if (!delta->x[0] && delta->wds <= 1) {
1.1962 + /* exact */
1.1963 +#ifdef SET_INEXACT
1.1964 + inexact = 0;
1.1965 +#endif
1.1966 + break;
1.1967 + }
1.1968 + if (rounding) {
1.1969 + if (dsign) {
1.1970 + adj = 1.;
1.1971 + goto apply_adj;
1.1972 + }
1.1973 + }
1.1974 + else if (!dsign) {
1.1975 + adj = -1.;
1.1976 + if (!word1(rv)
1.1977 + && !(word0(rv) & Frac_mask)) {
1.1978 + y = word0(rv) & Exp_mask;
1.1979 +#ifdef Avoid_Underflow
1.1980 + if (!scale || y > 2*P*Exp_msk1)
1.1981 +#else
1.1982 + if (y)
1.1983 +#endif
1.1984 + {
1.1985 + delta = lshift(PASS_STATE delta,Log2P);
1.1986 + if (cmp(delta, bs) <= 0)
1.1987 + adj = -0.5;
1.1988 + }
1.1989 + }
1.1990 + apply_adj:
1.1991 +#ifdef Avoid_Underflow
1.1992 + if (scale && (y = word0(rv) & Exp_mask)
1.1993 + <= 2*P*Exp_msk1)
1.1994 + word0(adj) += (2*P+1)*Exp_msk1 - y;
1.1995 +#else
1.1996 +#ifdef Sudden_Underflow
1.1997 + if ((word0(rv) & Exp_mask) <=
1.1998 + P*Exp_msk1) {
1.1999 + word0(rv) += P*Exp_msk1;
1.2000 + dval(rv) += adj*ulp(rv);
1.2001 + word0(rv) -= P*Exp_msk1;
1.2002 + }
1.2003 + else
1.2004 +#endif /*Sudden_Underflow*/
1.2005 +#endif /*Avoid_Underflow*/
1.2006 + dval(rv) += adj*ulp(rv);
1.2007 + }
1.2008 + break;
1.2009 + }
1.2010 + adj = ratio(delta, bs);
1.2011 + if (adj < 1.)
1.2012 + adj = 1.;
1.2013 + if (adj <= 0x7ffffffe) {
1.2014 + /* adj = rounding ? ceil(adj) : floor(adj); */
1.2015 + y = adj;
1.2016 + if (y != adj) {
1.2017 + if (!((rounding>>1) ^ dsign))
1.2018 + y++;
1.2019 + adj = y;
1.2020 + }
1.2021 + }
1.2022 +#ifdef Avoid_Underflow
1.2023 + if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
1.2024 + word0(adj) += (2*P+1)*Exp_msk1 - y;
1.2025 +#else
1.2026 +#ifdef Sudden_Underflow
1.2027 + if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1.2028 + word0(rv) += P*Exp_msk1;
1.2029 + adj *= ulp(rv);
1.2030 + if (dsign)
1.2031 + dval(rv) += adj;
1.2032 + else
1.2033 + dval(rv) -= adj;
1.2034 + word0(rv) -= P*Exp_msk1;
1.2035 + goto cont;
1.2036 + }
1.2037 +#endif /*Sudden_Underflow*/
1.2038 +#endif /*Avoid_Underflow*/
1.2039 + adj *= ulp(rv);
1.2040 + if (dsign)
1.2041 + dval(rv) += adj;
1.2042 + else
1.2043 + dval(rv) -= adj;
1.2044 + goto cont;
1.2045 + }
1.2046 +#endif /*Honor_FLT_ROUNDS*/
1.2047 +
1.2048 + if (i < 0) {
1.2049 + /* Error is less than half an ulp -- check for
1.2050 + * special case of mantissa a power of two.
1.2051 + */
1.2052 + if (dsign || word1(rv) || word0(rv) & Bndry_mask
1.2053 +#ifdef IEEE_Arith
1.2054 +#ifdef Avoid_Underflow
1.2055 + || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
1.2056 +#else
1.2057 + || (word0(rv) & Exp_mask) <= Exp_msk1
1.2058 +#endif
1.2059 +#endif
1.2060 + ) {
1.2061 +#ifdef SET_INEXACT
1.2062 + if (!delta->x[0] && delta->wds <= 1)
1.2063 + inexact = 0;
1.2064 +#endif
1.2065 + break;
1.2066 + }
1.2067 + if (!delta->x[0] && delta->wds <= 1) {
1.2068 + /* exact result */
1.2069 +#ifdef SET_INEXACT
1.2070 + inexact = 0;
1.2071 +#endif
1.2072 + break;
1.2073 + }
1.2074 + delta = lshift(PASS_STATE delta,Log2P);
1.2075 + if (cmp(delta, bs) > 0)
1.2076 + goto drop_down;
1.2077 + break;
1.2078 + }
1.2079 + if (i == 0) {
1.2080 + /* exactly half-way between */
1.2081 + if (dsign) {
1.2082 + if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1.2083 + && word1(rv) == (
1.2084 +#ifdef Avoid_Underflow
1.2085 + (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
1.2086 + ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1.2087 +#endif
1.2088 + 0xffffffff)) {
1.2089 + /*boundary case -- increment exponent*/
1.2090 + word0(rv) = (word0(rv) & Exp_mask)
1.2091 + + Exp_msk1
1.2092 +#ifdef IBM
1.2093 + | Exp_msk1 >> 4
1.2094 +#endif
1.2095 + ;
1.2096 + word1(rv) = 0;
1.2097 +#ifdef Avoid_Underflow
1.2098 + dsign = 0;
1.2099 +#endif
1.2100 + break;
1.2101 + }
1.2102 + }
1.2103 + else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1.2104 + drop_down:
1.2105 + /* boundary case -- decrement exponent */
1.2106 +#ifdef Sudden_Underflow /*{{*/
1.2107 + L = word0(rv) & Exp_mask;
1.2108 +#ifdef IBM
1.2109 + if (L < Exp_msk1)
1.2110 +#else
1.2111 +#ifdef Avoid_Underflow
1.2112 + if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
1.2113 +#else
1.2114 + if (L <= Exp_msk1)
1.2115 +#endif /*Avoid_Underflow*/
1.2116 +#endif /*IBM*/
1.2117 + goto undfl;
1.2118 + L -= Exp_msk1;
1.2119 +#else /*Sudden_Underflow}{*/
1.2120 +#ifdef Avoid_Underflow
1.2121 + if (scale) {
1.2122 + L = word0(rv) & Exp_mask;
1.2123 + if (L <= (2*P+1)*Exp_msk1) {
1.2124 + if (L > (P+2)*Exp_msk1)
1.2125 + /* round even ==> */
1.2126 + /* accept rv */
1.2127 + break;
1.2128 + /* rv = smallest denormal */
1.2129 + goto undfl;
1.2130 + }
1.2131 + }
1.2132 +#endif /*Avoid_Underflow*/
1.2133 + L = (word0(rv) & Exp_mask) - Exp_msk1;
1.2134 +#endif /*Sudden_Underflow}}*/
1.2135 + word0(rv) = L | Bndry_mask1;
1.2136 + word1(rv) = 0xffffffff;
1.2137 +#ifdef IBM
1.2138 + goto cont;
1.2139 +#else
1.2140 + break;
1.2141 +#endif
1.2142 + }
1.2143 +#ifndef ROUND_BIASED
1.2144 + if (!(word1(rv) & LSB))
1.2145 + break;
1.2146 +#endif
1.2147 + if (dsign)
1.2148 + dval(rv) += ulp(rv);
1.2149 +#ifndef ROUND_BIASED
1.2150 + else {
1.2151 + dval(rv) -= ulp(rv);
1.2152 +#ifndef Sudden_Underflow
1.2153 + if (!dval(rv))
1.2154 + goto undfl;
1.2155 +#endif
1.2156 + }
1.2157 +#ifdef Avoid_Underflow
1.2158 + dsign = 1 - dsign;
1.2159 +#endif
1.2160 +#endif
1.2161 + break;
1.2162 + }
1.2163 + if ((aadj = ratio(delta, bs)) <= 2.) {
1.2164 + if (dsign)
1.2165 + aadj = dval(aadj1) = 1.;
1.2166 + else if (word1(rv) || word0(rv) & Bndry_mask) {
1.2167 +#ifndef Sudden_Underflow
1.2168 + if (word1(rv) == Tiny1 && !word0(rv))
1.2169 + goto undfl;
1.2170 +#endif
1.2171 + aadj = 1.;
1.2172 + dval(aadj1) = -1.;
1.2173 + }
1.2174 + else {
1.2175 + /* special case -- power of FLT_RADIX to be */
1.2176 + /* rounded down... */
1.2177 +
1.2178 + if (aadj < 2./FLT_RADIX)
1.2179 + aadj = 1./FLT_RADIX;
1.2180 + else
1.2181 + aadj *= 0.5;
1.2182 + dval(aadj1) = -aadj;
1.2183 + }
1.2184 + }
1.2185 + else {
1.2186 + aadj *= 0.5;
1.2187 + dval(aadj1) = dsign ? aadj : -aadj;
1.2188 +#ifdef Check_FLT_ROUNDS
1.2189 + switch(Rounding) {
1.2190 + case 2: /* towards +infinity */
1.2191 + dval(aadj1) -= 0.5;
1.2192 + break;
1.2193 + case 0: /* towards 0 */
1.2194 + case 3: /* towards -infinity */
1.2195 + dval(aadj1) += 0.5;
1.2196 + }
1.2197 +#else
1.2198 + if (Flt_Rounds == 0)
1.2199 + dval(aadj1) += 0.5;
1.2200 +#endif /*Check_FLT_ROUNDS*/
1.2201 + }
1.2202 + y = word0(rv) & Exp_mask;
1.2203 +
1.2204 + /* Check for overflow */
1.2205 +
1.2206 + if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1.2207 + dval(rv0) = dval(rv);
1.2208 + word0(rv) -= P*Exp_msk1;
1.2209 + adj = dval(aadj1) * ulp(rv);
1.2210 + dval(rv) += adj;
1.2211 + if ((word0(rv) & Exp_mask) >=
1.2212 + Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1.2213 + if (word0(rv0) == Big0 && word1(rv0) == Big1)
1.2214 + goto ovfl;
1.2215 + word0(rv) = Big0;
1.2216 + word1(rv) = Big1;
1.2217 + goto cont;
1.2218 + }
1.2219 + else
1.2220 + word0(rv) += P*Exp_msk1;
1.2221 + }
1.2222 + else {
1.2223 +#ifdef Avoid_Underflow
1.2224 + if (scale && y <= 2*P*Exp_msk1) {
1.2225 + if (aadj <= 0x7fffffff) {
1.2226 + if ((z = (ULong) aadj) <= 0)
1.2227 + z = 1;
1.2228 + aadj = z;
1.2229 + dval(aadj1) = dsign ? aadj : -aadj;
1.2230 + }
1.2231 + word0(aadj1) += (2*P+1)*Exp_msk1 - y;
1.2232 + }
1.2233 + adj = dval(aadj1) * ulp(rv);
1.2234 + dval(rv) += adj;
1.2235 +#else
1.2236 +#ifdef Sudden_Underflow
1.2237 + if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1.2238 + dval(rv0) = dval(rv);
1.2239 + word0(rv) += P*Exp_msk1;
1.2240 + adj = dval(aadj1) * ulp(rv);
1.2241 + dval(rv) += adj;
1.2242 +#ifdef IBM
1.2243 + if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1.2244 +#else
1.2245 + if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1.2246 +#endif
1.2247 + {
1.2248 + if (word0(rv0) == Tiny0
1.2249 + && word1(rv0) == Tiny1)
1.2250 + goto undfl;
1.2251 + word0(rv) = Tiny0;
1.2252 + word1(rv) = Tiny1;
1.2253 + goto cont;
1.2254 + }
1.2255 + else
1.2256 + word0(rv) -= P*Exp_msk1;
1.2257 + }
1.2258 + else {
1.2259 + adj = dval(aadj1) * ulp(rv);
1.2260 + dval(rv) += adj;
1.2261 + }
1.2262 +#else /*Sudden_Underflow*/
1.2263 + /* Compute adj so that the IEEE rounding rules will
1.2264 + * correctly round rv + adj in some half-way cases.
1.2265 + * If rv * ulp(rv) is denormalized (i.e.,
1.2266 + * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1.2267 + * trouble from bits lost to denormalization;
1.2268 + * example: 1.2e-307 .
1.2269 + */
1.2270 + if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1.2271 + dval(aadj1) = (double)(int)(aadj + 0.5);
1.2272 + if (!dsign)
1.2273 + dval(aadj1) = -dval(aadj1);
1.2274 + }
1.2275 + adj = dval(aadj1) * ulp(rv);
1.2276 + dval(rv) += adj;
1.2277 +#endif /*Sudden_Underflow*/
1.2278 +#endif /*Avoid_Underflow*/
1.2279 + }
1.2280 + z = word0(rv) & Exp_mask;
1.2281 +#ifndef SET_INEXACT
1.2282 +#ifdef Avoid_Underflow
1.2283 + if (!scale)
1.2284 +#endif
1.2285 + if (y == z) {
1.2286 + /* Can we stop now? */
1.2287 + L = (Long)aadj;
1.2288 + aadj -= L;
1.2289 + /* The tolerances below are conservative. */
1.2290 + if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1.2291 + if (aadj < .4999999 || aadj > .5000001)
1.2292 + break;
1.2293 + }
1.2294 + else if (aadj < .4999999/FLT_RADIX)
1.2295 + break;
1.2296 + }
1.2297 +#endif
1.2298 + cont:
1.2299 + Bfree(PASS_STATE bb);
1.2300 + Bfree(PASS_STATE bd);
1.2301 + Bfree(PASS_STATE bs);
1.2302 + Bfree(PASS_STATE delta);
1.2303 + }
1.2304 +#ifdef SET_INEXACT
1.2305 + if (inexact) {
1.2306 + if (!oldinexact) {
1.2307 + word0(rv0) = Exp_1 + (70 << Exp_shift);
1.2308 + word1(rv0) = 0;
1.2309 + dval(rv0) += 1.;
1.2310 + }
1.2311 + }
1.2312 + else if (!oldinexact)
1.2313 + clear_inexact();
1.2314 +#endif
1.2315 +#ifdef Avoid_Underflow
1.2316 + if (scale) {
1.2317 + word0(rv0) = Exp_1 - 2*P*Exp_msk1;
1.2318 + word1(rv0) = 0;
1.2319 + dval(rv) *= dval(rv0);
1.2320 +#ifndef NO_ERRNO
1.2321 + /* try to avoid the bug of testing an 8087 register value */
1.2322 + if (word0(rv) == 0 && word1(rv) == 0)
1.2323 + errno = ERANGE;
1.2324 +#endif
1.2325 + }
1.2326 +#endif /* Avoid_Underflow */
1.2327 +#ifdef SET_INEXACT
1.2328 + if (inexact && !(word0(rv) & Exp_mask)) {
1.2329 + /* set underflow bit */
1.2330 + dval(rv0) = 1e-300;
1.2331 + dval(rv0) *= dval(rv0);
1.2332 + }
1.2333 +#endif
1.2334 + retfree:
1.2335 + Bfree(PASS_STATE bb);
1.2336 + Bfree(PASS_STATE bd);
1.2337 + Bfree(PASS_STATE bs);
1.2338 + Bfree(PASS_STATE bd0);
1.2339 + Bfree(PASS_STATE delta);
1.2340 + ret:
1.2341 + if (se)
1.2342 + *se = (char *)s;
1.2343 + return sign ? -dval(rv) : dval(rv);
1.2344 + }
1.2345 +
1.2346 + static int
1.2347 +quorem
1.2348 +#ifdef KR_headers
1.2349 + (b, S) Bigint *b, *S;
1.2350 +#else
1.2351 + (Bigint *b, Bigint *S)
1.2352 +#endif
1.2353 +{
1.2354 + int n;
1.2355 + ULong *bx, *bxe, q, *sx, *sxe;
1.2356 +#ifdef ULLong
1.2357 + ULLong borrow, carry, y, ys;
1.2358 +#else
1.2359 + ULong borrow, carry, y, ys;
1.2360 +#ifdef Pack_32
1.2361 + ULong si, z, zs;
1.2362 +#endif
1.2363 +#endif
1.2364 +
1.2365 + n = S->wds;
1.2366 +#ifdef DEBUG
1.2367 + /*debug*/ if (b->wds > n)
1.2368 + /*debug*/ Bug("oversize b in quorem");
1.2369 +#endif
1.2370 + if (b->wds < n)
1.2371 + return 0;
1.2372 + sx = S->x;
1.2373 + sxe = sx + --n;
1.2374 + bx = b->x;
1.2375 + bxe = bx + n;
1.2376 + q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1.2377 +#ifdef DEBUG
1.2378 + /*debug*/ if (q > 9)
1.2379 + /*debug*/ Bug("oversized quotient in quorem");
1.2380 +#endif
1.2381 + if (q) {
1.2382 + borrow = 0;
1.2383 + carry = 0;
1.2384 + do {
1.2385 +#ifdef ULLong
1.2386 + ys = *sx++ * (ULLong)q + carry;
1.2387 + carry = ys >> 32;
1.2388 + y = *bx - (ys & FFFFFFFF) - borrow;
1.2389 + borrow = y >> 32 & (ULong)1;
1.2390 + *bx++ = (ULong) y & FFFFFFFF;
1.2391 +#else
1.2392 +#ifdef Pack_32
1.2393 + si = *sx++;
1.2394 + ys = (si & 0xffff) * q + carry;
1.2395 + zs = (si >> 16) * q + (ys >> 16);
1.2396 + carry = zs >> 16;
1.2397 + y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1.2398 + borrow = (y & 0x10000) >> 16;
1.2399 + z = (*bx >> 16) - (zs & 0xffff) - borrow;
1.2400 + borrow = (z & 0x10000) >> 16;
1.2401 + Storeinc(bx, z, y);
1.2402 +#else
1.2403 + ys = *sx++ * q + carry;
1.2404 + carry = ys >> 16;
1.2405 + y = *bx - (ys & 0xffff) - borrow;
1.2406 + borrow = (y & 0x10000) >> 16;
1.2407 + *bx++ = y & 0xffff;
1.2408 +#endif
1.2409 +#endif
1.2410 + }
1.2411 + while(sx <= sxe);
1.2412 + if (!*bxe) {
1.2413 + bx = b->x;
1.2414 + while(--bxe > bx && !*bxe)
1.2415 + --n;
1.2416 + b->wds = n;
1.2417 + }
1.2418 + }
1.2419 + if (cmp(b, S) >= 0) {
1.2420 + q++;
1.2421 + borrow = 0;
1.2422 + carry = 0;
1.2423 + bx = b->x;
1.2424 + sx = S->x;
1.2425 + do {
1.2426 +#ifdef ULLong
1.2427 + ys = *sx++ + carry;
1.2428 + carry = ys >> 32;
1.2429 + y = *bx - (ys & FFFFFFFF) - borrow;
1.2430 + borrow = y >> 32 & (ULong)1;
1.2431 + *bx++ = (ULong) y & FFFFFFFF;
1.2432 +#else
1.2433 +#ifdef Pack_32
1.2434 + si = *sx++;
1.2435 + ys = (si & 0xffff) + carry;
1.2436 + zs = (si >> 16) + (ys >> 16);
1.2437 + carry = zs >> 16;
1.2438 + y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1.2439 + borrow = (y & 0x10000) >> 16;
1.2440 + z = (*bx >> 16) - (zs & 0xffff) - borrow;
1.2441 + borrow = (z & 0x10000) >> 16;
1.2442 + Storeinc(bx, z, y);
1.2443 +#else
1.2444 + ys = *sx++ + carry;
1.2445 + carry = ys >> 16;
1.2446 + y = *bx - (ys & 0xffff) - borrow;
1.2447 + borrow = (y & 0x10000) >> 16;
1.2448 + *bx++ = y & 0xffff;
1.2449 +#endif
1.2450 +#endif
1.2451 + }
1.2452 + while(sx <= sxe);
1.2453 + bx = b->x;
1.2454 + bxe = bx + n;
1.2455 + if (!*bxe) {
1.2456 + while(--bxe > bx && !*bxe)
1.2457 + --n;
1.2458 + b->wds = n;
1.2459 + }
1.2460 + }
1.2461 + return q;
1.2462 + }
1.2463 +
1.2464 +#if !defined(MULTIPLE_THREADS) && !defined(NO_GLOBAL_STATE)
1.2465 +#define USE_DTOA_RESULT 1
1.2466 + static char *dtoa_result;
1.2467 +#endif
1.2468 +
1.2469 + static char *
1.2470 +#ifdef KR_headers
1.2471 +rv_alloc(STATE_PARAM i) STATE_PARAM_DECL int i;
1.2472 +#else
1.2473 +rv_alloc(STATE_PARAM int i)
1.2474 +#endif
1.2475 +{
1.2476 + int j, k, *r;
1.2477 +
1.2478 + j = sizeof(ULong);
1.2479 + for(k = 0;
1.2480 + sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned) i;
1.2481 + j <<= 1)
1.2482 + k++;
1.2483 + r = (int*)Balloc(PASS_STATE k);
1.2484 + *r = k;
1.2485 + return
1.2486 +#ifdef USE_DTOA_RESULT
1.2487 + dtoa_result =
1.2488 +#endif
1.2489 + (char *)(r+1);
1.2490 + }
1.2491 +
1.2492 + static char *
1.2493 +#ifdef KR_headers
1.2494 +nrv_alloc(STATE_PARAM s, rve, n) STATE_PARAM_DECL char *s, **rve; int n;
1.2495 +#else
1.2496 +nrv_alloc(STATE_PARAM CONST char *s, char **rve, int n)
1.2497 +#endif
1.2498 +{
1.2499 + char *rv, *t;
1.2500 +
1.2501 + t = rv = rv_alloc(PASS_STATE n);
1.2502 + while((*t = *s++)) t++;
1.2503 + if (rve)
1.2504 + *rve = t;
1.2505 + return rv;
1.2506 + }
1.2507 +
1.2508 +/* freedtoa(s) must be used to free values s returned by dtoa
1.2509 + * when MULTIPLE_THREADS is #defined. It should be used in all cases,
1.2510 + * but for consistency with earlier versions of dtoa, it is optional
1.2511 + * when MULTIPLE_THREADS is not defined.
1.2512 + */
1.2513 +
1.2514 + static void
1.2515 +#ifdef KR_headers
1.2516 +freedtoa(STATE_PARAM s) STATE_PARAM_DECL char *s;
1.2517 +#else
1.2518 +freedtoa(STATE_PARAM char *s)
1.2519 +#endif
1.2520 +{
1.2521 + Bigint *b = (Bigint *)((int *)s - 1);
1.2522 + b->maxwds = 1 << (b->k = *(int*)b);
1.2523 + Bfree(PASS_STATE b);
1.2524 +#ifdef USE_DTOA_RESULT
1.2525 + if (s == dtoa_result)
1.2526 + dtoa_result = 0;
1.2527 +#endif
1.2528 + }
1.2529 +
1.2530 +/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1.2531 + *
1.2532 + * Inspired by "How to Print Floating-Point Numbers Accurately" by
1.2533 + * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1.2534 + *
1.2535 + * Modifications:
1.2536 + * 1. Rather than iterating, we use a simple numeric overestimate
1.2537 + * to determine k = floor(log10(d)). We scale relevant
1.2538 + * quantities using O(log2(k)) rather than O(k) multiplications.
1.2539 + * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1.2540 + * try to generate digits strictly left to right. Instead, we
1.2541 + * compute with fewer bits and propagate the carry if necessary
1.2542 + * when rounding the final digit up. This is often faster.
1.2543 + * 3. Under the assumption that input will be rounded nearest,
1.2544 + * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1.2545 + * That is, we allow equality in stopping tests when the
1.2546 + * round-nearest rule will give the same floating-point value
1.2547 + * as would satisfaction of the stopping test with strict
1.2548 + * inequality.
1.2549 + * 4. We remove common factors of powers of 2 from relevant
1.2550 + * quantities.
1.2551 + * 5. When converting floating-point integers less than 1e16,
1.2552 + * we use floating-point arithmetic rather than resorting
1.2553 + * to multiple-precision integers.
1.2554 + * 6. When asked to produce fewer than 15 digits, we first try
1.2555 + * to get by with floating-point arithmetic; we resort to
1.2556 + * multiple-precision integer arithmetic only if we cannot
1.2557 + * guarantee that the floating-point calculation has given
1.2558 + * the correctly rounded result. For k requested digits and
1.2559 + * "uniformly" distributed input, the probability is
1.2560 + * something like 10^(k-15) that we must resort to the Long
1.2561 + * calculation.
1.2562 + */
1.2563 +
1.2564 + static char *
1.2565 +dtoa
1.2566 +#ifdef KR_headers
1.2567 + (STATE_PARAM d, mode, ndigits, decpt, sign, rve)
1.2568 + STATE_PARAM_DECL U d; int mode, ndigits, *decpt, *sign; char **rve;
1.2569 +#else
1.2570 + (STATE_PARAM U d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1.2571 +#endif
1.2572 +{
1.2573 + /* Arguments ndigits, decpt, sign are similar to those
1.2574 + of ecvt and fcvt; trailing zeros are suppressed from
1.2575 + the returned string. If not null, *rve is set to point
1.2576 + to the end of the return value. If d is +-Infinity or NaN,
1.2577 + then *decpt is set to 9999.
1.2578 +
1.2579 + mode:
1.2580 + 0 ==> shortest string that yields d when read in
1.2581 + and rounded to nearest.
1.2582 + 1 ==> like 0, but with Steele & White stopping rule;
1.2583 + e.g. with IEEE P754 arithmetic , mode 0 gives
1.2584 + 1e23 whereas mode 1 gives 9.999999999999999e22.
1.2585 + 2 ==> max(1,ndigits) significant digits. This gives a
1.2586 + return value similar to that of ecvt, except
1.2587 + that trailing zeros are suppressed.
1.2588 + 3 ==> through ndigits past the decimal point. This
1.2589 + gives a return value similar to that from fcvt,
1.2590 + except that trailing zeros are suppressed, and
1.2591 + ndigits can be negative.
1.2592 + 4,5 ==> similar to 2 and 3, respectively, but (in
1.2593 + round-nearest mode) with the tests of mode 0 to
1.2594 + possibly return a shorter string that rounds to d.
1.2595 + With IEEE arithmetic and compilation with
1.2596 + -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
1.2597 + as modes 2 and 3 when FLT_ROUNDS != 1.
1.2598 + 6-9 ==> Debugging modes similar to mode - 4: don't try
1.2599 + fast floating-point estimate (if applicable).
1.2600 +
1.2601 + Values of mode other than 0-9 are treated as mode 0.
1.2602 +
1.2603 + Sufficient space is allocated to the return value
1.2604 + to hold the suppressed trailing zeros.
1.2605 + */
1.2606 +
1.2607 + int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1.2608 + j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1.2609 + spec_case, try_quick;
1.2610 + Long L;
1.2611 +#ifndef Sudden_Underflow
1.2612 + int denorm;
1.2613 + ULong x;
1.2614 +#endif
1.2615 + Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1.2616 + U d2, eps;
1.2617 + double ds;
1.2618 + char *s, *s0;
1.2619 +#ifdef Honor_FLT_ROUNDS
1.2620 + int rounding;
1.2621 +#endif
1.2622 +#ifdef SET_INEXACT
1.2623 + int inexact, oldinexact;
1.2624 +#endif
1.2625 +
1.2626 +#ifdef __GNUC__
1.2627 + ilim = ilim1 = 0;
1.2628 + mlo = NULL;
1.2629 +#endif
1.2630 +
1.2631 +#ifdef USE_DTOA_RESULT
1.2632 + if (dtoa_result) {
1.2633 + freedtoa(PASS_STATE dtoa_result);
1.2634 + dtoa_result = 0;
1.2635 + }
1.2636 +#endif
1.2637 +
1.2638 + if (word0(d) & Sign_bit) {
1.2639 + /* set sign for everything, including 0's and NaNs */
1.2640 + *sign = 1;
1.2641 + word0(d) &= ~Sign_bit; /* clear sign bit */
1.2642 + }
1.2643 + else
1.2644 + *sign = 0;
1.2645 +
1.2646 +#if defined(IEEE_Arith) + defined(VAX)
1.2647 +#ifdef IEEE_Arith
1.2648 + if ((word0(d) & Exp_mask) == Exp_mask)
1.2649 +#else
1.2650 + if (word0(d) == 0x8000)
1.2651 +#endif
1.2652 + {
1.2653 + /* Infinity or NaN */
1.2654 + *decpt = 9999;
1.2655 +#ifdef IEEE_Arith
1.2656 + if (!word1(d) && !(word0(d) & 0xfffff))
1.2657 + return nrv_alloc(PASS_STATE "Infinity", rve, 8);
1.2658 +#endif
1.2659 + return nrv_alloc(PASS_STATE "NaN", rve, 3);
1.2660 + }
1.2661 +#endif
1.2662 +#ifdef IBM
1.2663 + dval(d) += 0; /* normalize */
1.2664 +#endif
1.2665 + if (!dval(d)) {
1.2666 + *decpt = 1;
1.2667 + return nrv_alloc(PASS_STATE "0", rve, 1);
1.2668 + }
1.2669 +
1.2670 +#ifdef SET_INEXACT
1.2671 + try_quick = oldinexact = get_inexact();
1.2672 + inexact = 1;
1.2673 +#endif
1.2674 +#ifdef Honor_FLT_ROUNDS
1.2675 + if ((rounding = Flt_Rounds) >= 2) {
1.2676 + if (*sign)
1.2677 + rounding = rounding == 2 ? 0 : 2;
1.2678 + else
1.2679 + if (rounding != 2)
1.2680 + rounding = 0;
1.2681 + }
1.2682 +#endif
1.2683 +
1.2684 + b = d2b(PASS_STATE d, &be, &bbits);
1.2685 +#ifdef Sudden_Underflow
1.2686 + i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1.2687 +#else
1.2688 + if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
1.2689 +#endif
1.2690 + dval(d2) = dval(d);
1.2691 + word0(d2) &= Frac_mask1;
1.2692 + word0(d2) |= Exp_11;
1.2693 +#ifdef IBM
1.2694 + if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1.2695 + dval(d2) /= 1 << j;
1.2696 +#endif
1.2697 +
1.2698 + /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1.2699 + * log10(x) = log(x) / log(10)
1.2700 + * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1.2701 + * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1.2702 + *
1.2703 + * This suggests computing an approximation k to log10(d) by
1.2704 + *
1.2705 + * k = (i - Bias)*0.301029995663981
1.2706 + * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1.2707 + *
1.2708 + * We want k to be too large rather than too small.
1.2709 + * The error in the first-order Taylor series approximation
1.2710 + * is in our favor, so we just round up the constant enough
1.2711 + * to compensate for any error in the multiplication of
1.2712 + * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1.2713 + * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1.2714 + * adding 1e-13 to the constant term more than suffices.
1.2715 + * Hence we adjust the constant term to 0.1760912590558.
1.2716 + * (We could get a more accurate k by invoking log10,
1.2717 + * but this is probably not worthwhile.)
1.2718 + */
1.2719 +
1.2720 + i -= Bias;
1.2721 +#ifdef IBM
1.2722 + i <<= 2;
1.2723 + i += j;
1.2724 +#endif
1.2725 +#ifndef Sudden_Underflow
1.2726 + denorm = 0;
1.2727 + }
1.2728 + else {
1.2729 + /* d is denormalized */
1.2730 +
1.2731 + i = bbits + be + (Bias + (P-1) - 1);
1.2732 + x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
1.2733 + : word1(d) << (32 - i);
1.2734 + dval(d2) = x;
1.2735 + word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1.2736 + i -= (Bias + (P-1) - 1) + 1;
1.2737 + denorm = 1;
1.2738 + }
1.2739 +#endif
1.2740 + ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1.2741 + k = (int)ds;
1.2742 + if (ds < 0. && ds != k)
1.2743 + k--; /* want k = floor(ds) */
1.2744 + k_check = 1;
1.2745 + if (k >= 0 && k <= Ten_pmax) {
1.2746 + if (dval(d) < tens[k])
1.2747 + k--;
1.2748 + k_check = 0;
1.2749 + }
1.2750 + j = bbits - i - 1;
1.2751 + if (j >= 0) {
1.2752 + b2 = 0;
1.2753 + s2 = j;
1.2754 + }
1.2755 + else {
1.2756 + b2 = -j;
1.2757 + s2 = 0;
1.2758 + }
1.2759 + if (k >= 0) {
1.2760 + b5 = 0;
1.2761 + s5 = k;
1.2762 + s2 += k;
1.2763 + }
1.2764 + else {
1.2765 + b2 -= k;
1.2766 + b5 = -k;
1.2767 + s5 = 0;
1.2768 + }
1.2769 + if (mode < 0 || mode > 9)
1.2770 + mode = 0;
1.2771 +
1.2772 +#ifndef SET_INEXACT
1.2773 +#ifdef Check_FLT_ROUNDS
1.2774 + try_quick = Rounding == 1;
1.2775 +#else
1.2776 + try_quick = 1;
1.2777 +#endif
1.2778 +#endif /*SET_INEXACT*/
1.2779 +
1.2780 + if (mode > 5) {
1.2781 + mode -= 4;
1.2782 + try_quick = 0;
1.2783 + }
1.2784 + leftright = 1;
1.2785 + switch(mode) {
1.2786 + case 0:
1.2787 + case 1:
1.2788 + ilim = ilim1 = -1;
1.2789 + i = 18;
1.2790 + ndigits = 0;
1.2791 + break;
1.2792 + case 2:
1.2793 + leftright = 0;
1.2794 + /* no break */
1.2795 + case 4:
1.2796 + if (ndigits <= 0)
1.2797 + ndigits = 1;
1.2798 + ilim = ilim1 = i = ndigits;
1.2799 + break;
1.2800 + case 3:
1.2801 + leftright = 0;
1.2802 + /* no break */
1.2803 + case 5:
1.2804 + i = ndigits + k + 1;
1.2805 + ilim = i;
1.2806 + ilim1 = i - 1;
1.2807 + if (i <= 0)
1.2808 + i = 1;
1.2809 + }
1.2810 + s = s0 = rv_alloc(PASS_STATE i);
1.2811 +
1.2812 +#ifdef Honor_FLT_ROUNDS
1.2813 + if (mode > 1 && rounding != 1)
1.2814 + leftright = 0;
1.2815 +#endif
1.2816 +
1.2817 + if (ilim >= 0 && ilim <= Quick_max && try_quick) {
1.2818 +
1.2819 + /* Try to get by with floating-point arithmetic. */
1.2820 +
1.2821 + i = 0;
1.2822 + dval(d2) = dval(d);
1.2823 + k0 = k;
1.2824 + ilim0 = ilim;
1.2825 + ieps = 2; /* conservative */
1.2826 + if (k > 0) {
1.2827 + ds = tens[k&0xf];
1.2828 + j = k >> 4;
1.2829 + if (j & Bletch) {
1.2830 + /* prevent overflows */
1.2831 + j &= Bletch - 1;
1.2832 + dval(d) /= bigtens[n_bigtens-1];
1.2833 + ieps++;
1.2834 + }
1.2835 + for(; j; j >>= 1, i++)
1.2836 + if (j & 1) {
1.2837 + ieps++;
1.2838 + ds *= bigtens[i];
1.2839 + }
1.2840 + dval(d) /= ds;
1.2841 + }
1.2842 + else if ((j1 = -k)) {
1.2843 + dval(d) *= tens[j1 & 0xf];
1.2844 + for(j = j1 >> 4; j; j >>= 1, i++)
1.2845 + if (j & 1) {
1.2846 + ieps++;
1.2847 + dval(d) *= bigtens[i];
1.2848 + }
1.2849 + }
1.2850 + if (k_check && dval(d) < 1. && ilim > 0) {
1.2851 + if (ilim1 <= 0)
1.2852 + goto fast_failed;
1.2853 + ilim = ilim1;
1.2854 + k--;
1.2855 + dval(d) *= 10.;
1.2856 + ieps++;
1.2857 + }
1.2858 + dval(eps) = ieps*dval(d) + 7.;
1.2859 + word0(eps) -= (P-1)*Exp_msk1;
1.2860 + if (ilim == 0) {
1.2861 + S = mhi = 0;
1.2862 + dval(d) -= 5.;
1.2863 + if (dval(d) > dval(eps))
1.2864 + goto one_digit;
1.2865 + if (dval(d) < -dval(eps))
1.2866 + goto no_digits;
1.2867 + goto fast_failed;
1.2868 + }
1.2869 +#ifndef No_leftright
1.2870 + if (leftright) {
1.2871 + /* Use Steele & White method of only
1.2872 + * generating digits needed.
1.2873 + */
1.2874 + dval(eps) = 0.5/tens[ilim-1] - dval(eps);
1.2875 + for(i = 0;;) {
1.2876 + L = (ULong) dval(d);
1.2877 + dval(d) -= L;
1.2878 + *s++ = '0' + (int)L;
1.2879 + if (dval(d) < dval(eps))
1.2880 + goto ret1;
1.2881 + if (1. - dval(d) < dval(eps))
1.2882 + goto bump_up;
1.2883 + if (++i >= ilim)
1.2884 + break;
1.2885 + dval(eps) *= 10.;
1.2886 + dval(d) *= 10.;
1.2887 + }
1.2888 + }
1.2889 + else {
1.2890 +#endif
1.2891 + /* Generate ilim digits, then fix them up. */
1.2892 + dval(eps) *= tens[ilim-1];
1.2893 + for(i = 1;; i++, dval(d) *= 10.) {
1.2894 + L = (Long)(dval(d));
1.2895 + if (!(dval(d) -= L))
1.2896 + ilim = i;
1.2897 + *s++ = '0' + (int)L;
1.2898 + if (i == ilim) {
1.2899 + if (dval(d) > 0.5 + dval(eps))
1.2900 + goto bump_up;
1.2901 + else if (dval(d) < 0.5 - dval(eps)) {
1.2902 + while(*--s == '0');
1.2903 + s++;
1.2904 + goto ret1;
1.2905 + }
1.2906 + break;
1.2907 + }
1.2908 + }
1.2909 +#ifndef No_leftright
1.2910 + }
1.2911 +#endif
1.2912 + fast_failed:
1.2913 + s = s0;
1.2914 + dval(d) = dval(d2);
1.2915 + k = k0;
1.2916 + ilim = ilim0;
1.2917 + }
1.2918 +
1.2919 + /* Do we have a "small" integer? */
1.2920 +
1.2921 + if (be >= 0 && k <= Int_max) {
1.2922 + /* Yes. */
1.2923 + ds = tens[k];
1.2924 + if (ndigits < 0 && ilim <= 0) {
1.2925 + S = mhi = 0;
1.2926 + if (ilim < 0 || dval(d) < 5*ds)
1.2927 + goto no_digits;
1.2928 + goto one_digit;
1.2929 + }
1.2930 + for(i = 1;; i++, dval(d) *= 10.) {
1.2931 + L = (Long)(dval(d) / ds);
1.2932 + dval(d) -= L*ds;
1.2933 +#ifdef Check_FLT_ROUNDS
1.2934 + /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1.2935 + if (dval(d) < 0) {
1.2936 + L--;
1.2937 + dval(d) += ds;
1.2938 + }
1.2939 +#endif
1.2940 + *s++ = '0' + (int)L;
1.2941 + if (!dval(d)) {
1.2942 +#ifdef SET_INEXACT
1.2943 + inexact = 0;
1.2944 +#endif
1.2945 + break;
1.2946 + }
1.2947 + if (i == ilim) {
1.2948 +#ifdef Honor_FLT_ROUNDS
1.2949 + if (mode > 1)
1.2950 + switch(rounding) {
1.2951 + case 0: goto ret1;
1.2952 + case 2: goto bump_up;
1.2953 + }
1.2954 +#endif
1.2955 + dval(d) += dval(d);
1.2956 + if (dval(d) > ds || (dval(d) == ds && L & 1)) {
1.2957 + bump_up:
1.2958 + while(*--s == '9')
1.2959 + if (s == s0) {
1.2960 + k++;
1.2961 + *s = '0';
1.2962 + break;
1.2963 + }
1.2964 + ++*s++;
1.2965 + }
1.2966 + break;
1.2967 + }
1.2968 + }
1.2969 + goto ret1;
1.2970 + }
1.2971 +
1.2972 + m2 = b2;
1.2973 + m5 = b5;
1.2974 + mhi = mlo = 0;
1.2975 + if (leftright) {
1.2976 + i =
1.2977 +#ifndef Sudden_Underflow
1.2978 + denorm ? be + (Bias + (P-1) - 1 + 1) :
1.2979 +#endif
1.2980 +#ifdef IBM
1.2981 + 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
1.2982 +#else
1.2983 + 1 + P - bbits;
1.2984 +#endif
1.2985 + b2 += i;
1.2986 + s2 += i;
1.2987 + mhi = i2b(PASS_STATE 1);
1.2988 + }
1.2989 + if (m2 > 0 && s2 > 0) {
1.2990 + i = m2 < s2 ? m2 : s2;
1.2991 + b2 -= i;
1.2992 + m2 -= i;
1.2993 + s2 -= i;
1.2994 + }
1.2995 + if (b5 > 0) {
1.2996 + if (leftright) {
1.2997 + if (m5 > 0) {
1.2998 + mhi = pow5mult(PASS_STATE mhi, m5);
1.2999 + b1 = mult(PASS_STATE mhi, b);
1.3000 + Bfree(PASS_STATE b);
1.3001 + b = b1;
1.3002 + }
1.3003 + if ((j = b5 - m5))
1.3004 + b = pow5mult(PASS_STATE b, j);
1.3005 + }
1.3006 + else
1.3007 + b = pow5mult(PASS_STATE b, b5);
1.3008 + }
1.3009 + S = i2b(PASS_STATE 1);
1.3010 + if (s5 > 0)
1.3011 + S = pow5mult(PASS_STATE S, s5);
1.3012 +
1.3013 + /* Check for special case that d is a normalized power of 2. */
1.3014 +
1.3015 + spec_case = 0;
1.3016 + if ((mode < 2 || leftright)
1.3017 +#ifdef Honor_FLT_ROUNDS
1.3018 + && rounding == 1
1.3019 +#endif
1.3020 + ) {
1.3021 + if (!word1(d) && !(word0(d) & Bndry_mask)
1.3022 +#ifndef Sudden_Underflow
1.3023 + && word0(d) & (Exp_mask & ~Exp_msk1)
1.3024 +#endif
1.3025 + ) {
1.3026 + /* The special case */
1.3027 + b2 += Log2P;
1.3028 + s2 += Log2P;
1.3029 + spec_case = 1;
1.3030 + }
1.3031 + }
1.3032 +
1.3033 + /* Arrange for convenient computation of quotients:
1.3034 + * shift left if necessary so divisor has 4 leading 0 bits.
1.3035 + *
1.3036 + * Perhaps we should just compute leading 28 bits of S once
1.3037 + * and for all and pass them and a shift to quorem, so it
1.3038 + * can do shifts and ors to compute the numerator for q.
1.3039 + */
1.3040 +#ifdef Pack_32
1.3041 + if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
1.3042 + i = 32 - i;
1.3043 +#else
1.3044 + if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
1.3045 + i = 16 - i;
1.3046 +#endif
1.3047 + if (i > 4) {
1.3048 + i -= 4;
1.3049 + b2 += i;
1.3050 + m2 += i;
1.3051 + s2 += i;
1.3052 + }
1.3053 + else if (i < 4) {
1.3054 + i += 28;
1.3055 + b2 += i;
1.3056 + m2 += i;
1.3057 + s2 += i;
1.3058 + }
1.3059 + if (b2 > 0)
1.3060 + b = lshift(PASS_STATE b, b2);
1.3061 + if (s2 > 0)
1.3062 + S = lshift(PASS_STATE S, s2);
1.3063 + if (k_check) {
1.3064 + if (cmp(b,S) < 0) {
1.3065 + k--;
1.3066 + b = multadd(PASS_STATE b, 10, 0); /* we botched the k estimate */
1.3067 + if (leftright)
1.3068 + mhi = multadd(PASS_STATE mhi, 10, 0);
1.3069 + ilim = ilim1;
1.3070 + }
1.3071 + }
1.3072 + if (ilim <= 0 && (mode == 3 || mode == 5)) {
1.3073 + if (ilim < 0 || cmp(b,S = multadd(PASS_STATE S,5,0)) < 0) {
1.3074 + /* no digits, fcvt style */
1.3075 + no_digits:
1.3076 + /* MOZILLA CHANGE: Always return a non-empty string. */
1.3077 + *s++ = '0';
1.3078 + k = 0;
1.3079 + goto ret;
1.3080 + }
1.3081 + one_digit:
1.3082 + *s++ = '1';
1.3083 + k++;
1.3084 + goto ret;
1.3085 + }
1.3086 + if (leftright) {
1.3087 + if (m2 > 0)
1.3088 + mhi = lshift(PASS_STATE mhi, m2);
1.3089 +
1.3090 + /* Compute mlo -- check for special case
1.3091 + * that d is a normalized power of 2.
1.3092 + */
1.3093 +
1.3094 + mlo = mhi;
1.3095 + if (spec_case) {
1.3096 + mhi = Balloc(PASS_STATE mhi->k);
1.3097 + Bcopy(mhi, mlo);
1.3098 + mhi = lshift(PASS_STATE mhi, Log2P);
1.3099 + }
1.3100 +
1.3101 + for(i = 1;;i++) {
1.3102 + dig = quorem(b,S) + '0';
1.3103 + /* Do we yet have the shortest decimal string
1.3104 + * that will round to d?
1.3105 + */
1.3106 + j = cmp(b, mlo);
1.3107 + delta = diff(PASS_STATE S, mhi);
1.3108 + j1 = delta->sign ? 1 : cmp(b, delta);
1.3109 + Bfree(PASS_STATE delta);
1.3110 +#ifndef ROUND_BIASED
1.3111 + if (j1 == 0 && mode != 1 && !(word1(d) & 1)
1.3112 +#ifdef Honor_FLT_ROUNDS
1.3113 + && rounding >= 1
1.3114 +#endif
1.3115 + ) {
1.3116 + if (dig == '9')
1.3117 + goto round_9_up;
1.3118 + if (j > 0)
1.3119 + dig++;
1.3120 +#ifdef SET_INEXACT
1.3121 + else if (!b->x[0] && b->wds <= 1)
1.3122 + inexact = 0;
1.3123 +#endif
1.3124 + *s++ = dig;
1.3125 + goto ret;
1.3126 + }
1.3127 +#endif
1.3128 + if (j < 0 || (j == 0 && mode != 1
1.3129 +#ifndef ROUND_BIASED
1.3130 + && !(word1(d) & 1)
1.3131 +#endif
1.3132 + )) {
1.3133 + if (!b->x[0] && b->wds <= 1) {
1.3134 +#ifdef SET_INEXACT
1.3135 + inexact = 0;
1.3136 +#endif
1.3137 + goto accept_dig;
1.3138 + }
1.3139 +#ifdef Honor_FLT_ROUNDS
1.3140 + if (mode > 1)
1.3141 + switch(rounding) {
1.3142 + case 0: goto accept_dig;
1.3143 + case 2: goto keep_dig;
1.3144 + }
1.3145 +#endif /*Honor_FLT_ROUNDS*/
1.3146 + if (j1 > 0) {
1.3147 + b = lshift(PASS_STATE b, 1);
1.3148 + j1 = cmp(b, S);
1.3149 + if ((j1 > 0 || (j1 == 0 && dig & 1))
1.3150 + && dig++ == '9')
1.3151 + goto round_9_up;
1.3152 + }
1.3153 + accept_dig:
1.3154 + *s++ = dig;
1.3155 + goto ret;
1.3156 + }
1.3157 + if (j1 > 0) {
1.3158 +#ifdef Honor_FLT_ROUNDS
1.3159 + if (!rounding)
1.3160 + goto accept_dig;
1.3161 +#endif
1.3162 + if (dig == '9') { /* possible if i == 1 */
1.3163 + round_9_up:
1.3164 + *s++ = '9';
1.3165 + goto roundoff;
1.3166 + }
1.3167 + *s++ = dig + 1;
1.3168 + goto ret;
1.3169 + }
1.3170 +#ifdef Honor_FLT_ROUNDS
1.3171 + keep_dig:
1.3172 +#endif
1.3173 + *s++ = dig;
1.3174 + if (i == ilim)
1.3175 + break;
1.3176 + b = multadd(PASS_STATE b, 10, 0);
1.3177 + if (mlo == mhi)
1.3178 + mlo = mhi = multadd(PASS_STATE mhi, 10, 0);
1.3179 + else {
1.3180 + mlo = multadd(PASS_STATE mlo, 10, 0);
1.3181 + mhi = multadd(PASS_STATE mhi, 10, 0);
1.3182 + }
1.3183 + }
1.3184 + }
1.3185 + else
1.3186 + for(i = 1;; i++) {
1.3187 + *s++ = dig = quorem(b,S) + '0';
1.3188 + if (!b->x[0] && b->wds <= 1) {
1.3189 +#ifdef SET_INEXACT
1.3190 + inexact = 0;
1.3191 +#endif
1.3192 + goto ret;
1.3193 + }
1.3194 + if (i >= ilim)
1.3195 + break;
1.3196 + b = multadd(PASS_STATE b, 10, 0);
1.3197 + }
1.3198 +
1.3199 + /* Round off last digit */
1.3200 +
1.3201 +#ifdef Honor_FLT_ROUNDS
1.3202 + switch(rounding) {
1.3203 + case 0: goto trimzeros;
1.3204 + case 2: goto roundoff;
1.3205 + }
1.3206 +#endif
1.3207 + b = lshift(PASS_STATE b, 1);
1.3208 + j = cmp(b, S);
1.3209 + if (j >= 0) { /* ECMA compatible rounding needed by Spidermonkey */
1.3210 + roundoff:
1.3211 + while(*--s == '9')
1.3212 + if (s == s0) {
1.3213 + k++;
1.3214 + *s++ = '1';
1.3215 + goto ret;
1.3216 + }
1.3217 + ++*s++;
1.3218 + }
1.3219 + else {
1.3220 +#ifdef Honor_FLT_ROUNDS
1.3221 + trimzeros:
1.3222 +#endif
1.3223 + while(*--s == '0');
1.3224 + s++;
1.3225 + }
1.3226 + ret:
1.3227 + Bfree(PASS_STATE S);
1.3228 + if (mhi) {
1.3229 + if (mlo && mlo != mhi)
1.3230 + Bfree(PASS_STATE mlo);
1.3231 + Bfree(PASS_STATE mhi);
1.3232 + }
1.3233 + ret1:
1.3234 +#ifdef SET_INEXACT
1.3235 + if (inexact) {
1.3236 + if (!oldinexact) {
1.3237 + word0(d) = Exp_1 + (70 << Exp_shift);
1.3238 + word1(d) = 0;
1.3239 + dval(d) += 1.;
1.3240 + }
1.3241 + }
1.3242 + else if (!oldinexact)
1.3243 + clear_inexact();
1.3244 +#endif
1.3245 + Bfree(PASS_STATE b);
1.3246 + *s = 0;
1.3247 + *decpt = k + 1;
1.3248 + if (rve)
1.3249 + *rve = s;
1.3250 + return s0;
1.3251 + }
