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