js/src/dtoa.c

Sat, 03 Jan 2015 20:18:00 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Sat, 03 Jan 2015 20:18:00 +0100
branch
TOR_BUG_3246
changeset 7
129ffea94266
permissions
-rw-r--r--

Conditionally enable double key logic according to:
private browsing mode or privacy.thirdparty.isolate preference and
implement in GetCookieStringCommon and FindCookie where it counts...
With some reservations of how to convince FindCookie users to test
condition and pass a nullptr when disabling double key logic.

michael@0 1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
michael@0 2 /****************************************************************
michael@0 3 *
michael@0 4 * The author of this software is David M. Gay.
michael@0 5 *
michael@0 6 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
michael@0 7 *
michael@0 8 * Permission to use, copy, modify, and distribute this software for any
michael@0 9 * purpose without fee is hereby granted, provided that this entire notice
michael@0 10 * is included in all copies of any software which is or includes a copy
michael@0 11 * or modification of this software and in all copies of the supporting
michael@0 12 * documentation for such software.
michael@0 13 *
michael@0 14 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
michael@0 15 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
michael@0 16 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
michael@0 17 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
michael@0 18 *
michael@0 19 ***************************************************************/
michael@0 20
michael@0 21 /* Please send bug reports to David M. Gay (dmg at acm dot org,
michael@0 22 * with " at " changed at "@" and " dot " changed to "."). */
michael@0 23
michael@0 24 /* On a machine with IEEE extended-precision registers, it is
michael@0 25 * necessary to specify double-precision (53-bit) rounding precision
michael@0 26 * before invoking strtod or dtoa. If the machine uses (the equivalent
michael@0 27 * of) Intel 80x87 arithmetic, the call
michael@0 28 * _control87(PC_53, MCW_PC);
michael@0 29 * does this with many compilers. Whether this or another call is
michael@0 30 * appropriate depends on the compiler; for this to work, it may be
michael@0 31 * necessary to #include "float.h" or another system-dependent header
michael@0 32 * file.
michael@0 33 */
michael@0 34
michael@0 35 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
michael@0 36 *
michael@0 37 * This strtod returns a nearest machine number to the input decimal
michael@0 38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
michael@0 39 * broken by the IEEE round-even rule. Otherwise ties are broken by
michael@0 40 * biased rounding (add half and chop).
michael@0 41 *
michael@0 42 * Inspired loosely by William D. Clinger's paper "How to Read Floating
michael@0 43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
michael@0 44 *
michael@0 45 * Modifications:
michael@0 46 *
michael@0 47 * 1. We only require IEEE, IBM, or VAX double-precision
michael@0 48 * arithmetic (not IEEE double-extended).
michael@0 49 * 2. We get by with floating-point arithmetic in a case that
michael@0 50 * Clinger missed -- when we're computing d * 10^n
michael@0 51 * for a small integer d and the integer n is not too
michael@0 52 * much larger than 22 (the maximum integer k for which
michael@0 53 * we can represent 10^k exactly), we may be able to
michael@0 54 * compute (d*10^k) * 10^(e-k) with just one roundoff.
michael@0 55 * 3. Rather than a bit-at-a-time adjustment of the binary
michael@0 56 * result in the hard case, we use floating-point
michael@0 57 * arithmetic to determine the adjustment to within
michael@0 58 * one bit; only in really hard cases do we need to
michael@0 59 * compute a second residual.
michael@0 60 * 4. Because of 3., we don't need a large table of powers of 10
michael@0 61 * for ten-to-e (just some small tables, e.g. of 10^k
michael@0 62 * for 0 <= k <= 22).
michael@0 63 */
michael@0 64
michael@0 65 /*
michael@0 66 * #define IEEE_8087 for IEEE-arithmetic machines where the least
michael@0 67 * significant byte has the lowest address.
michael@0 68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
michael@0 69 * significant byte has the lowest address.
michael@0 70 * #define Long int on machines with 32-bit ints and 64-bit longs.
michael@0 71 * #define IBM for IBM mainframe-style floating-point arithmetic.
michael@0 72 * #define VAX for VAX-style floating-point arithmetic (D_floating).
michael@0 73 * #define No_leftright to omit left-right logic in fast floating-point
michael@0 74 * computation of dtoa.
michael@0 75 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
michael@0 76 * and strtod and dtoa should round accordingly.
michael@0 77 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
michael@0 78 * and Honor_FLT_ROUNDS is not #defined.
michael@0 79 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
michael@0 80 * that use extended-precision instructions to compute rounded
michael@0 81 * products and quotients) with IBM.
michael@0 82 * #define ROUND_BIASED for IEEE-format with biased rounding.
michael@0 83 * #define Inaccurate_Divide for IEEE-format with correctly rounded
michael@0 84 * products but inaccurate quotients, e.g., for Intel i860.
michael@0 85 * #define NO_LONG_LONG on machines that do not have a "long long"
michael@0 86 * integer type (of >= 64 bits). On such machines, you can
michael@0 87 * #define Just_16 to store 16 bits per 32-bit Long when doing
michael@0 88 * high-precision integer arithmetic. Whether this speeds things
michael@0 89 * up or slows things down depends on the machine and the number
michael@0 90 * being converted. If long long is available and the name is
michael@0 91 * something other than "long long", #define Llong to be the name,
michael@0 92 * and if "unsigned Llong" does not work as an unsigned version of
michael@0 93 * Llong, #define #ULLong to be the corresponding unsigned type.
michael@0 94 * #define KR_headers for old-style C function headers.
michael@0 95 * #define Bad_float_h if your system lacks a float.h or if it does not
michael@0 96 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
michael@0 97 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
michael@0 98 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
michael@0 99 * if memory is available and otherwise does something you deem
michael@0 100 * appropriate. If MALLOC is undefined, malloc will be invoked
michael@0 101 * directly -- and assumed always to succeed. Similarly, if you
michael@0 102 * want something other than the system's free() to be called to
michael@0 103 * recycle memory acquired from MALLOC, #define FREE to be the
michael@0 104 * name of the alternate routine. (Unless you #define
michael@0 105 * NO_GLOBAL_STATE and call destroydtoa, FREE or free is only
michael@0 106 * called in pathological cases, e.g., in a dtoa call after a dtoa
michael@0 107 * return in mode 3 with thousands of digits requested.)
michael@0 108 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
michael@0 109 * memory allocations from a private pool of memory when possible.
michael@0 110 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
michael@0 111 * unless #defined to be a different length. This default length
michael@0 112 * suffices to get rid of MALLOC calls except for unusual cases,
michael@0 113 * such as decimal-to-binary conversion of a very long string of
michael@0 114 * digits. The longest string dtoa can return is about 751 bytes
michael@0 115 * long. For conversions by strtod of strings of 800 digits and
michael@0 116 * all dtoa conversions in single-threaded executions with 8-byte
michael@0 117 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
michael@0 118 * pointers, PRIVATE_MEM >= 7112 appears adequate.
michael@0 119 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
michael@0 120 * multiple threads. In this case, you must provide (or suitably
michael@0 121 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
michael@0 122 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
michael@0 123 * in pow5mult, ensures lazy evaluation of only one copy of high
michael@0 124 * powers of 5; omitting this lock would introduce a small
michael@0 125 * probability of wasting memory, but would otherwise be harmless.)
michael@0 126 * You must also invoke freedtoa(s) to free the value s returned by
michael@0 127 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
michael@0 128 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
michael@0 129 * avoids underflows on inputs whose result does not underflow.
michael@0 130 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
michael@0 131 * floating-point numbers and flushes underflows to zero rather
michael@0 132 * than implementing gradual underflow, then you must also #define
michael@0 133 * Sudden_Underflow.
michael@0 134 * #define USE_LOCALE to use the current locale's decimal_point value.
michael@0 135 * #define SET_INEXACT if IEEE arithmetic is being used and extra
michael@0 136 * computation should be done to set the inexact flag when the
michael@0 137 * result is inexact and avoid setting inexact when the result
michael@0 138 * is exact. In this case, dtoa.c must be compiled in
michael@0 139 * an environment, perhaps provided by #include "dtoa.c" in a
michael@0 140 * suitable wrapper, that defines two functions,
michael@0 141 * int get_inexact(void);
michael@0 142 * void clear_inexact(void);
michael@0 143 * such that get_inexact() returns a nonzero value if the
michael@0 144 * inexact bit is already set, and clear_inexact() sets the
michael@0 145 * inexact bit to 0. When SET_INEXACT is #defined, strtod
michael@0 146 * also does extra computations to set the underflow and overflow
michael@0 147 * flags when appropriate (i.e., when the result is tiny and
michael@0 148 * inexact or when it is a numeric value rounded to +-infinity).
michael@0 149 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
michael@0 150 * the result overflows to +-Infinity or underflows to 0.
michael@0 151 * #define NO_GLOBAL_STATE to avoid defining any non-const global or
michael@0 152 * static variables. Instead the necessary state is stored in an
michael@0 153 * opaque struct, DtoaState, a pointer to which must be passed to
michael@0 154 * every entry point. Two new functions are added to the API:
michael@0 155 * DtoaState *newdtoa(void);
michael@0 156 * void destroydtoa(DtoaState *);
michael@0 157 */
michael@0 158
michael@0 159 #ifndef Long
michael@0 160 #define Long long
michael@0 161 #endif
michael@0 162 #ifndef ULong
michael@0 163 typedef unsigned Long ULong;
michael@0 164 #endif
michael@0 165
michael@0 166 #ifdef DEBUG
michael@0 167 #include <stdio.h>
michael@0 168 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
michael@0 169 #endif
michael@0 170
michael@0 171 #include <stdlib.h>
michael@0 172 #include <string.h>
michael@0 173
michael@0 174 #ifdef USE_LOCALE
michael@0 175 #include <locale.h>
michael@0 176 #endif
michael@0 177
michael@0 178 #ifdef MALLOC
michael@0 179 #ifdef KR_headers
michael@0 180 extern char *MALLOC();
michael@0 181 #else
michael@0 182 extern void *MALLOC(size_t);
michael@0 183 #endif
michael@0 184 #else
michael@0 185 #define MALLOC malloc
michael@0 186 #endif
michael@0 187
michael@0 188 #ifndef FREE
michael@0 189 #define FREE free
michael@0 190 #endif
michael@0 191
michael@0 192 #ifndef Omit_Private_Memory
michael@0 193 #ifndef PRIVATE_MEM
michael@0 194 #define PRIVATE_MEM 2304
michael@0 195 #endif
michael@0 196 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
michael@0 197 #endif
michael@0 198
michael@0 199 #undef IEEE_Arith
michael@0 200 #undef Avoid_Underflow
michael@0 201 #ifdef IEEE_MC68k
michael@0 202 #define IEEE_Arith
michael@0 203 #endif
michael@0 204 #ifdef IEEE_8087
michael@0 205 #define IEEE_Arith
michael@0 206 #endif
michael@0 207
michael@0 208 #include <errno.h>
michael@0 209
michael@0 210 #ifdef Bad_float_h
michael@0 211
michael@0 212 #ifdef IEEE_Arith
michael@0 213 #define DBL_DIG 15
michael@0 214 #define DBL_MAX_10_EXP 308
michael@0 215 #define DBL_MAX_EXP 1024
michael@0 216 #define FLT_RADIX 2
michael@0 217 #endif /*IEEE_Arith*/
michael@0 218
michael@0 219 #ifdef IBM
michael@0 220 #define DBL_DIG 16
michael@0 221 #define DBL_MAX_10_EXP 75
michael@0 222 #define DBL_MAX_EXP 63
michael@0 223 #define FLT_RADIX 16
michael@0 224 #define DBL_MAX 7.2370055773322621e+75
michael@0 225 #endif
michael@0 226
michael@0 227 #ifdef VAX
michael@0 228 #define DBL_DIG 16
michael@0 229 #define DBL_MAX_10_EXP 38
michael@0 230 #define DBL_MAX_EXP 127
michael@0 231 #define FLT_RADIX 2
michael@0 232 #define DBL_MAX 1.7014118346046923e+38
michael@0 233 #endif
michael@0 234
michael@0 235 #ifndef LONG_MAX
michael@0 236 #define LONG_MAX 2147483647
michael@0 237 #endif
michael@0 238
michael@0 239 #else /* ifndef Bad_float_h */
michael@0 240 #include <float.h>
michael@0 241 #endif /* Bad_float_h */
michael@0 242
michael@0 243 #ifndef __MATH_H__
michael@0 244 #include <math.h>
michael@0 245 #endif
michael@0 246
michael@0 247 #ifndef CONST
michael@0 248 #ifdef KR_headers
michael@0 249 #define CONST /* blank */
michael@0 250 #else
michael@0 251 #define CONST const
michael@0 252 #endif
michael@0 253 #endif
michael@0 254
michael@0 255 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
michael@0 256 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
michael@0 257 #endif
michael@0 258
michael@0 259 typedef union { double d; ULong L[2]; } U;
michael@0 260
michael@0 261 #define dval(x) ((x).d)
michael@0 262 #ifdef IEEE_8087
michael@0 263 #define word0(x) ((x).L[1])
michael@0 264 #define word1(x) ((x).L[0])
michael@0 265 #else
michael@0 266 #define word0(x) ((x).L[0])
michael@0 267 #define word1(x) ((x).L[1])
michael@0 268 #endif
michael@0 269
michael@0 270 /* The following definition of Storeinc is appropriate for MIPS processors.
michael@0 271 * An alternative that might be better on some machines is
michael@0 272 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
michael@0 273 */
michael@0 274 #if defined(IEEE_8087) + defined(VAX)
michael@0 275 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
michael@0 276 ((unsigned short *)a)[0] = (unsigned short)c, a++)
michael@0 277 #else
michael@0 278 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
michael@0 279 ((unsigned short *)a)[1] = (unsigned short)c, a++)
michael@0 280 #endif
michael@0 281
michael@0 282 /* #define P DBL_MANT_DIG */
michael@0 283 /* Ten_pmax = floor(P*log(2)/log(5)) */
michael@0 284 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
michael@0 285 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
michael@0 286 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
michael@0 287
michael@0 288 #ifdef IEEE_Arith
michael@0 289 #define Exp_shift 20
michael@0 290 #define Exp_shift1 20
michael@0 291 #define Exp_msk1 0x100000
michael@0 292 #define Exp_msk11 0x100000
michael@0 293 #define Exp_mask 0x7ff00000
michael@0 294 #define P 53
michael@0 295 #define Bias 1023
michael@0 296 #define Emin (-1022)
michael@0 297 #define Exp_1 0x3ff00000
michael@0 298 #define Exp_11 0x3ff00000
michael@0 299 #define Ebits 11
michael@0 300 #define Frac_mask 0xfffff
michael@0 301 #define Frac_mask1 0xfffff
michael@0 302 #define Ten_pmax 22
michael@0 303 #define Bletch 0x10
michael@0 304 #define Bndry_mask 0xfffff
michael@0 305 #define Bndry_mask1 0xfffff
michael@0 306 #define LSB 1
michael@0 307 #define Sign_bit 0x80000000
michael@0 308 #define Log2P 1
michael@0 309 #define Tiny0 0
michael@0 310 #define Tiny1 1
michael@0 311 #define Quick_max 14
michael@0 312 #define Int_max 14
michael@0 313 #ifndef NO_IEEE_Scale
michael@0 314 #define Avoid_Underflow
michael@0 315 #ifdef Flush_Denorm /* debugging option */
michael@0 316 #undef Sudden_Underflow
michael@0 317 #endif
michael@0 318 #endif
michael@0 319
michael@0 320 #ifndef Flt_Rounds
michael@0 321 #ifdef FLT_ROUNDS
michael@0 322 #define Flt_Rounds FLT_ROUNDS
michael@0 323 #else
michael@0 324 #define Flt_Rounds 1
michael@0 325 #endif
michael@0 326 #endif /*Flt_Rounds*/
michael@0 327
michael@0 328 #ifdef Honor_FLT_ROUNDS
michael@0 329 #define Rounding rounding
michael@0 330 #undef Check_FLT_ROUNDS
michael@0 331 #define Check_FLT_ROUNDS
michael@0 332 #else
michael@0 333 #define Rounding Flt_Rounds
michael@0 334 #endif
michael@0 335
michael@0 336 #else /* ifndef IEEE_Arith */
michael@0 337 #undef Check_FLT_ROUNDS
michael@0 338 #undef Honor_FLT_ROUNDS
michael@0 339 #undef SET_INEXACT
michael@0 340 #undef Sudden_Underflow
michael@0 341 #define Sudden_Underflow
michael@0 342 #ifdef IBM
michael@0 343 #undef Flt_Rounds
michael@0 344 #define Flt_Rounds 0
michael@0 345 #define Exp_shift 24
michael@0 346 #define Exp_shift1 24
michael@0 347 #define Exp_msk1 0x1000000
michael@0 348 #define Exp_msk11 0x1000000
michael@0 349 #define Exp_mask 0x7f000000
michael@0 350 #define P 14
michael@0 351 #define Bias 65
michael@0 352 #define Exp_1 0x41000000
michael@0 353 #define Exp_11 0x41000000
michael@0 354 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
michael@0 355 #define Frac_mask 0xffffff
michael@0 356 #define Frac_mask1 0xffffff
michael@0 357 #define Bletch 4
michael@0 358 #define Ten_pmax 22
michael@0 359 #define Bndry_mask 0xefffff
michael@0 360 #define Bndry_mask1 0xffffff
michael@0 361 #define LSB 1
michael@0 362 #define Sign_bit 0x80000000
michael@0 363 #define Log2P 4
michael@0 364 #define Tiny0 0x100000
michael@0 365 #define Tiny1 0
michael@0 366 #define Quick_max 14
michael@0 367 #define Int_max 15
michael@0 368 #else /* VAX */
michael@0 369 #undef Flt_Rounds
michael@0 370 #define Flt_Rounds 1
michael@0 371 #define Exp_shift 23
michael@0 372 #define Exp_shift1 7
michael@0 373 #define Exp_msk1 0x80
michael@0 374 #define Exp_msk11 0x800000
michael@0 375 #define Exp_mask 0x7f80
michael@0 376 #define P 56
michael@0 377 #define Bias 129
michael@0 378 #define Exp_1 0x40800000
michael@0 379 #define Exp_11 0x4080
michael@0 380 #define Ebits 8
michael@0 381 #define Frac_mask 0x7fffff
michael@0 382 #define Frac_mask1 0xffff007f
michael@0 383 #define Ten_pmax 24
michael@0 384 #define Bletch 2
michael@0 385 #define Bndry_mask 0xffff007f
michael@0 386 #define Bndry_mask1 0xffff007f
michael@0 387 #define LSB 0x10000
michael@0 388 #define Sign_bit 0x8000
michael@0 389 #define Log2P 1
michael@0 390 #define Tiny0 0x80
michael@0 391 #define Tiny1 0
michael@0 392 #define Quick_max 15
michael@0 393 #define Int_max 15
michael@0 394 #endif /* IBM, VAX */
michael@0 395 #endif /* IEEE_Arith */
michael@0 396
michael@0 397 #ifndef IEEE_Arith
michael@0 398 #define ROUND_BIASED
michael@0 399 #endif
michael@0 400
michael@0 401 #ifdef RND_PRODQUOT
michael@0 402 #define rounded_product(a,b) a = rnd_prod(a, b)
michael@0 403 #define rounded_quotient(a,b) a = rnd_quot(a, b)
michael@0 404 #ifdef KR_headers
michael@0 405 extern double rnd_prod(), rnd_quot();
michael@0 406 #else
michael@0 407 extern double rnd_prod(double, double), rnd_quot(double, double);
michael@0 408 #endif
michael@0 409 #else
michael@0 410 #define rounded_product(a,b) a *= b
michael@0 411 #define rounded_quotient(a,b) a /= b
michael@0 412 #endif
michael@0 413
michael@0 414 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
michael@0 415 #define Big1 0xffffffff
michael@0 416
michael@0 417 #ifndef Pack_32
michael@0 418 #define Pack_32
michael@0 419 #endif
michael@0 420
michael@0 421 #ifdef KR_headers
michael@0 422 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
michael@0 423 #else
michael@0 424 #define FFFFFFFF 0xffffffffUL
michael@0 425 #endif
michael@0 426
michael@0 427 #ifdef NO_LONG_LONG
michael@0 428 #undef ULLong
michael@0 429 #ifdef Just_16
michael@0 430 #undef Pack_32
michael@0 431 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
michael@0 432 * This makes some inner loops simpler and sometimes saves work
michael@0 433 * during multiplications, but it often seems to make things slightly
michael@0 434 * slower. Hence the default is now to store 32 bits per Long.
michael@0 435 */
michael@0 436 #endif
michael@0 437 #else /* long long available */
michael@0 438 #ifndef Llong
michael@0 439 #define Llong long long
michael@0 440 #endif
michael@0 441 #ifndef ULLong
michael@0 442 #define ULLong unsigned Llong
michael@0 443 #endif
michael@0 444 #endif /* NO_LONG_LONG */
michael@0 445
michael@0 446 #ifndef MULTIPLE_THREADS
michael@0 447 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
michael@0 448 #define FREE_DTOA_LOCK(n) /*nothing*/
michael@0 449 #endif
michael@0 450
michael@0 451 #define Kmax 7
michael@0 452
michael@0 453 struct
michael@0 454 Bigint {
michael@0 455 struct Bigint *next;
michael@0 456 int k, maxwds, sign, wds;
michael@0 457 ULong x[1];
michael@0 458 };
michael@0 459
michael@0 460 typedef struct Bigint Bigint;
michael@0 461
michael@0 462 #ifdef NO_GLOBAL_STATE
michael@0 463 #ifdef MULTIPLE_THREADS
michael@0 464 #error "cannot have both NO_GLOBAL_STATE and MULTIPLE_THREADS"
michael@0 465 #endif
michael@0 466 struct
michael@0 467 DtoaState {
michael@0 468 #define DECLARE_GLOBAL_STATE /* nothing */
michael@0 469 #else
michael@0 470 #define DECLARE_GLOBAL_STATE static
michael@0 471 #endif
michael@0 472
michael@0 473 DECLARE_GLOBAL_STATE Bigint *freelist[Kmax+1];
michael@0 474 DECLARE_GLOBAL_STATE Bigint *p5s;
michael@0 475 #ifndef Omit_Private_Memory
michael@0 476 DECLARE_GLOBAL_STATE double private_mem[PRIVATE_mem];
michael@0 477 DECLARE_GLOBAL_STATE double *pmem_next
michael@0 478 #ifndef NO_GLOBAL_STATE
michael@0 479 = private_mem
michael@0 480 #endif
michael@0 481 ;
michael@0 482 #endif
michael@0 483 #ifdef NO_GLOBAL_STATE
michael@0 484 };
michael@0 485 typedef struct DtoaState DtoaState;
michael@0 486 #ifdef KR_headers
michael@0 487 #define STATE_PARAM state,
michael@0 488 #define STATE_PARAM_DECL DtoaState *state;
michael@0 489 #else
michael@0 490 #define STATE_PARAM DtoaState *state,
michael@0 491 #endif
michael@0 492 #define PASS_STATE state,
michael@0 493 #define GET_STATE(field) (state->field)
michael@0 494
michael@0 495 static DtoaState *
michael@0 496 newdtoa(void)
michael@0 497 {
michael@0 498 DtoaState *state = (DtoaState *) MALLOC(sizeof(DtoaState));
michael@0 499 if (state) {
michael@0 500 memset(state, 0, sizeof(DtoaState));
michael@0 501 #ifndef Omit_Private_Memory
michael@0 502 state->pmem_next = state->private_mem;
michael@0 503 #endif
michael@0 504 }
michael@0 505 return state;
michael@0 506 }
michael@0 507
michael@0 508 static void
michael@0 509 destroydtoa
michael@0 510 #ifdef KR_headers
michael@0 511 (state) STATE_PARAM_DECL
michael@0 512 #else
michael@0 513 (DtoaState *state)
michael@0 514 #endif
michael@0 515 {
michael@0 516 int i;
michael@0 517 Bigint *v, *next;
michael@0 518
michael@0 519 for (i = 0; i <= Kmax; i++) {
michael@0 520 for (v = GET_STATE(freelist)[i]; v; v = next) {
michael@0 521 next = v->next;
michael@0 522 #ifndef Omit_Private_Memory
michael@0 523 if ((double*)v < GET_STATE(private_mem) ||
michael@0 524 (double*)v >= GET_STATE(private_mem) + PRIVATE_mem)
michael@0 525 #endif
michael@0 526 FREE((void*)v);
michael@0 527 }
michael@0 528 }
michael@0 529 FREE((void *)state);
michael@0 530 }
michael@0 531
michael@0 532 #else
michael@0 533 #define STATE_PARAM /* nothing */
michael@0 534 #define STATE_PARAM_DECL /* nothing */
michael@0 535 #define PASS_STATE /* nothing */
michael@0 536 #define GET_STATE(name) name
michael@0 537 #endif
michael@0 538
michael@0 539 static Bigint *
michael@0 540 Balloc
michael@0 541 #ifdef KR_headers
michael@0 542 (STATE_PARAM k) STATE_PARAM_DECL int k;
michael@0 543 #else
michael@0 544 (STATE_PARAM int k)
michael@0 545 #endif
michael@0 546 {
michael@0 547 int x;
michael@0 548 Bigint *rv;
michael@0 549 #ifndef Omit_Private_Memory
michael@0 550 size_t len;
michael@0 551 #endif
michael@0 552
michael@0 553 ACQUIRE_DTOA_LOCK(0);
michael@0 554 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
michael@0 555 /* but this case seems very unlikely. */
michael@0 556 if (k <= Kmax && (rv = GET_STATE(freelist)[k]))
michael@0 557 GET_STATE(freelist)[k] = rv->next;
michael@0 558 else {
michael@0 559 x = 1 << k;
michael@0 560 #ifdef Omit_Private_Memory
michael@0 561 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
michael@0 562 #else
michael@0 563 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
michael@0 564 /sizeof(double);
michael@0 565 if (k <= Kmax && GET_STATE(pmem_next) - GET_STATE(private_mem) + len <= PRIVATE_mem) {
michael@0 566 rv = (Bigint*)GET_STATE(pmem_next);
michael@0 567 GET_STATE(pmem_next) += len;
michael@0 568 }
michael@0 569 else
michael@0 570 rv = (Bigint*)MALLOC(len*sizeof(double));
michael@0 571 #endif
michael@0 572 rv->k = k;
michael@0 573 rv->maxwds = x;
michael@0 574 }
michael@0 575 FREE_DTOA_LOCK(0);
michael@0 576 rv->sign = rv->wds = 0;
michael@0 577 return rv;
michael@0 578 }
michael@0 579
michael@0 580 static void
michael@0 581 Bfree
michael@0 582 #ifdef KR_headers
michael@0 583 (STATE_PARAM v) STATE_PARAM_DECL Bigint *v;
michael@0 584 #else
michael@0 585 (STATE_PARAM Bigint *v)
michael@0 586 #endif
michael@0 587 {
michael@0 588 if (v) {
michael@0 589 if (v->k > Kmax)
michael@0 590 FREE((void*)v);
michael@0 591 else {
michael@0 592 ACQUIRE_DTOA_LOCK(0);
michael@0 593 v->next = GET_STATE(freelist)[v->k];
michael@0 594 GET_STATE(freelist)[v->k] = v;
michael@0 595 FREE_DTOA_LOCK(0);
michael@0 596 }
michael@0 597 }
michael@0 598 }
michael@0 599
michael@0 600 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
michael@0 601 y->wds*sizeof(Long) + 2*sizeof(int))
michael@0 602
michael@0 603 static Bigint *
michael@0 604 multadd
michael@0 605 #ifdef KR_headers
michael@0 606 (STATE_PARAM b, m, a) STATE_PARAM_DECL Bigint *b; int m, a;
michael@0 607 #else
michael@0 608 (STATE_PARAM Bigint *b, int m, int a) /* multiply by m and add a */
michael@0 609 #endif
michael@0 610 {
michael@0 611 int i, wds;
michael@0 612 #ifdef ULLong
michael@0 613 ULong *x;
michael@0 614 ULLong carry, y;
michael@0 615 #else
michael@0 616 ULong carry, *x, y;
michael@0 617 #ifdef Pack_32
michael@0 618 ULong xi, z;
michael@0 619 #endif
michael@0 620 #endif
michael@0 621 Bigint *b1;
michael@0 622
michael@0 623 wds = b->wds;
michael@0 624 x = b->x;
michael@0 625 i = 0;
michael@0 626 carry = a;
michael@0 627 do {
michael@0 628 #ifdef ULLong
michael@0 629 y = *x * (ULLong)m + carry;
michael@0 630 carry = y >> 32;
michael@0 631 *x++ = (ULong) y & FFFFFFFF;
michael@0 632 #else
michael@0 633 #ifdef Pack_32
michael@0 634 xi = *x;
michael@0 635 y = (xi & 0xffff) * m + carry;
michael@0 636 z = (xi >> 16) * m + (y >> 16);
michael@0 637 carry = z >> 16;
michael@0 638 *x++ = (z << 16) + (y & 0xffff);
michael@0 639 #else
michael@0 640 y = *x * m + carry;
michael@0 641 carry = y >> 16;
michael@0 642 *x++ = y & 0xffff;
michael@0 643 #endif
michael@0 644 #endif
michael@0 645 }
michael@0 646 while(++i < wds);
michael@0 647 if (carry) {
michael@0 648 if (wds >= b->maxwds) {
michael@0 649 b1 = Balloc(PASS_STATE b->k+1);
michael@0 650 Bcopy(b1, b);
michael@0 651 Bfree(PASS_STATE b);
michael@0 652 b = b1;
michael@0 653 }
michael@0 654 b->x[wds++] = (ULong) carry;
michael@0 655 b->wds = wds;
michael@0 656 }
michael@0 657 return b;
michael@0 658 }
michael@0 659
michael@0 660 static Bigint *
michael@0 661 s2b
michael@0 662 #ifdef KR_headers
michael@0 663 (STATE_PARAM s, nd0, nd, y9) STATE_PARAM_DECL CONST char *s; int nd0, nd; ULong y9;
michael@0 664 #else
michael@0 665 (STATE_PARAM CONST char *s, int nd0, int nd, ULong y9)
michael@0 666 #endif
michael@0 667 {
michael@0 668 Bigint *b;
michael@0 669 int i, k;
michael@0 670 Long x, y;
michael@0 671
michael@0 672 x = (nd + 8) / 9;
michael@0 673 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
michael@0 674 #ifdef Pack_32
michael@0 675 b = Balloc(PASS_STATE k);
michael@0 676 b->x[0] = y9;
michael@0 677 b->wds = 1;
michael@0 678 #else
michael@0 679 b = Balloc(PASS_STATE k+1);
michael@0 680 b->x[0] = y9 & 0xffff;
michael@0 681 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
michael@0 682 #endif
michael@0 683
michael@0 684 i = 9;
michael@0 685 if (9 < nd0) {
michael@0 686 s += 9;
michael@0 687 do b = multadd(PASS_STATE b, 10, *s++ - '0');
michael@0 688 while(++i < nd0);
michael@0 689 s++;
michael@0 690 }
michael@0 691 else
michael@0 692 s += 10;
michael@0 693 for(; i < nd; i++)
michael@0 694 b = multadd(PASS_STATE b, 10, *s++ - '0');
michael@0 695 return b;
michael@0 696 }
michael@0 697
michael@0 698 static int
michael@0 699 hi0bits
michael@0 700 #ifdef KR_headers
michael@0 701 (x) ULong x;
michael@0 702 #else
michael@0 703 (ULong x)
michael@0 704 #endif
michael@0 705 {
michael@0 706 int k = 0;
michael@0 707
michael@0 708 if (!(x & 0xffff0000)) {
michael@0 709 k = 16;
michael@0 710 x <<= 16;
michael@0 711 }
michael@0 712 if (!(x & 0xff000000)) {
michael@0 713 k += 8;
michael@0 714 x <<= 8;
michael@0 715 }
michael@0 716 if (!(x & 0xf0000000)) {
michael@0 717 k += 4;
michael@0 718 x <<= 4;
michael@0 719 }
michael@0 720 if (!(x & 0xc0000000)) {
michael@0 721 k += 2;
michael@0 722 x <<= 2;
michael@0 723 }
michael@0 724 if (!(x & 0x80000000)) {
michael@0 725 k++;
michael@0 726 if (!(x & 0x40000000))
michael@0 727 return 32;
michael@0 728 }
michael@0 729 return k;
michael@0 730 }
michael@0 731
michael@0 732 static int
michael@0 733 lo0bits
michael@0 734 #ifdef KR_headers
michael@0 735 (y) ULong *y;
michael@0 736 #else
michael@0 737 (ULong *y)
michael@0 738 #endif
michael@0 739 {
michael@0 740 int k;
michael@0 741 ULong x = *y;
michael@0 742
michael@0 743 if (x & 7) {
michael@0 744 if (x & 1)
michael@0 745 return 0;
michael@0 746 if (x & 2) {
michael@0 747 *y = x >> 1;
michael@0 748 return 1;
michael@0 749 }
michael@0 750 *y = x >> 2;
michael@0 751 return 2;
michael@0 752 }
michael@0 753 k = 0;
michael@0 754 if (!(x & 0xffff)) {
michael@0 755 k = 16;
michael@0 756 x >>= 16;
michael@0 757 }
michael@0 758 if (!(x & 0xff)) {
michael@0 759 k += 8;
michael@0 760 x >>= 8;
michael@0 761 }
michael@0 762 if (!(x & 0xf)) {
michael@0 763 k += 4;
michael@0 764 x >>= 4;
michael@0 765 }
michael@0 766 if (!(x & 0x3)) {
michael@0 767 k += 2;
michael@0 768 x >>= 2;
michael@0 769 }
michael@0 770 if (!(x & 1)) {
michael@0 771 k++;
michael@0 772 x >>= 1;
michael@0 773 if (!x)
michael@0 774 return 32;
michael@0 775 }
michael@0 776 *y = x;
michael@0 777 return k;
michael@0 778 }
michael@0 779
michael@0 780 static Bigint *
michael@0 781 i2b
michael@0 782 #ifdef KR_headers
michael@0 783 (STATE_PARAM i) STATE_PARAM_DECL int i;
michael@0 784 #else
michael@0 785 (STATE_PARAM int i)
michael@0 786 #endif
michael@0 787 {
michael@0 788 Bigint *b;
michael@0 789
michael@0 790 b = Balloc(PASS_STATE 1);
michael@0 791 b->x[0] = i;
michael@0 792 b->wds = 1;
michael@0 793 return b;
michael@0 794 }
michael@0 795
michael@0 796 static Bigint *
michael@0 797 mult
michael@0 798 #ifdef KR_headers
michael@0 799 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
michael@0 800 #else
michael@0 801 (STATE_PARAM Bigint *a, Bigint *b)
michael@0 802 #endif
michael@0 803 {
michael@0 804 Bigint *c;
michael@0 805 int k, wa, wb, wc;
michael@0 806 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
michael@0 807 ULong y;
michael@0 808 #ifdef ULLong
michael@0 809 ULLong carry, z;
michael@0 810 #else
michael@0 811 ULong carry, z;
michael@0 812 #ifdef Pack_32
michael@0 813 ULong z2;
michael@0 814 #endif
michael@0 815 #endif
michael@0 816
michael@0 817 if (a->wds < b->wds) {
michael@0 818 c = a;
michael@0 819 a = b;
michael@0 820 b = c;
michael@0 821 }
michael@0 822 k = a->k;
michael@0 823 wa = a->wds;
michael@0 824 wb = b->wds;
michael@0 825 wc = wa + wb;
michael@0 826 if (wc > a->maxwds)
michael@0 827 k++;
michael@0 828 c = Balloc(PASS_STATE k);
michael@0 829 for(x = c->x, xa = x + wc; x < xa; x++)
michael@0 830 *x = 0;
michael@0 831 xa = a->x;
michael@0 832 xae = xa + wa;
michael@0 833 xb = b->x;
michael@0 834 xbe = xb + wb;
michael@0 835 xc0 = c->x;
michael@0 836 #ifdef ULLong
michael@0 837 for(; xb < xbe; xc0++) {
michael@0 838 if ((y = *xb++)) {
michael@0 839 x = xa;
michael@0 840 xc = xc0;
michael@0 841 carry = 0;
michael@0 842 do {
michael@0 843 z = *x++ * (ULLong)y + *xc + carry;
michael@0 844 carry = z >> 32;
michael@0 845 *xc++ = (ULong) z & FFFFFFFF;
michael@0 846 }
michael@0 847 while(x < xae);
michael@0 848 *xc = (ULong) carry;
michael@0 849 }
michael@0 850 }
michael@0 851 #else
michael@0 852 #ifdef Pack_32
michael@0 853 for(; xb < xbe; xb++, xc0++) {
michael@0 854 if (y = *xb & 0xffff) {
michael@0 855 x = xa;
michael@0 856 xc = xc0;
michael@0 857 carry = 0;
michael@0 858 do {
michael@0 859 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
michael@0 860 carry = z >> 16;
michael@0 861 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
michael@0 862 carry = z2 >> 16;
michael@0 863 Storeinc(xc, z2, z);
michael@0 864 }
michael@0 865 while(x < xae);
michael@0 866 *xc = carry;
michael@0 867 }
michael@0 868 if (y = *xb >> 16) {
michael@0 869 x = xa;
michael@0 870 xc = xc0;
michael@0 871 carry = 0;
michael@0 872 z2 = *xc;
michael@0 873 do {
michael@0 874 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
michael@0 875 carry = z >> 16;
michael@0 876 Storeinc(xc, z, z2);
michael@0 877 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
michael@0 878 carry = z2 >> 16;
michael@0 879 }
michael@0 880 while(x < xae);
michael@0 881 *xc = z2;
michael@0 882 }
michael@0 883 }
michael@0 884 #else
michael@0 885 for(; xb < xbe; xc0++) {
michael@0 886 if (y = *xb++) {
michael@0 887 x = xa;
michael@0 888 xc = xc0;
michael@0 889 carry = 0;
michael@0 890 do {
michael@0 891 z = *x++ * y + *xc + carry;
michael@0 892 carry = z >> 16;
michael@0 893 *xc++ = z & 0xffff;
michael@0 894 }
michael@0 895 while(x < xae);
michael@0 896 *xc = carry;
michael@0 897 }
michael@0 898 }
michael@0 899 #endif
michael@0 900 #endif
michael@0 901 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
michael@0 902 c->wds = wc;
michael@0 903 return c;
michael@0 904 }
michael@0 905
michael@0 906 static Bigint *
michael@0 907 pow5mult
michael@0 908 #ifdef KR_headers
michael@0 909 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
michael@0 910 #else
michael@0 911 (STATE_PARAM Bigint *b, int k)
michael@0 912 #endif
michael@0 913 {
michael@0 914 Bigint *b1, *p5, *p51;
michael@0 915 int i;
michael@0 916 static CONST int p05[3] = { 5, 25, 125 };
michael@0 917
michael@0 918 if ((i = k & 3))
michael@0 919 b = multadd(PASS_STATE b, p05[i-1], 0);
michael@0 920
michael@0 921 if (!(k >>= 2))
michael@0 922 return b;
michael@0 923 if (!(p5 = GET_STATE(p5s))) {
michael@0 924 /* first time */
michael@0 925 #ifdef MULTIPLE_THREADS
michael@0 926 ACQUIRE_DTOA_LOCK(1);
michael@0 927 if (!(p5 = p5s)) {
michael@0 928 p5 = p5s = i2b(625);
michael@0 929 p5->next = 0;
michael@0 930 }
michael@0 931 FREE_DTOA_LOCK(1);
michael@0 932 #else
michael@0 933 p5 = GET_STATE(p5s) = i2b(PASS_STATE 625);
michael@0 934 p5->next = 0;
michael@0 935 #endif
michael@0 936 }
michael@0 937 for(;;) {
michael@0 938 if (k & 1) {
michael@0 939 b1 = mult(PASS_STATE b, p5);
michael@0 940 Bfree(PASS_STATE b);
michael@0 941 b = b1;
michael@0 942 }
michael@0 943 if (!(k >>= 1))
michael@0 944 break;
michael@0 945 if (!(p51 = p5->next)) {
michael@0 946 #ifdef MULTIPLE_THREADS
michael@0 947 ACQUIRE_DTOA_LOCK(1);
michael@0 948 if (!(p51 = p5->next)) {
michael@0 949 p51 = p5->next = mult(p5,p5);
michael@0 950 p51->next = 0;
michael@0 951 }
michael@0 952 FREE_DTOA_LOCK(1);
michael@0 953 #else
michael@0 954 p51 = p5->next = mult(PASS_STATE p5,p5);
michael@0 955 p51->next = 0;
michael@0 956 #endif
michael@0 957 }
michael@0 958 p5 = p51;
michael@0 959 }
michael@0 960 return b;
michael@0 961 }
michael@0 962
michael@0 963 static Bigint *
michael@0 964 lshift
michael@0 965 #ifdef KR_headers
michael@0 966 (STATE_PARAM b, k) STATE_PARAM_DECL Bigint *b; int k;
michael@0 967 #else
michael@0 968 (STATE_PARAM Bigint *b, int k)
michael@0 969 #endif
michael@0 970 {
michael@0 971 int i, k1, n, n1;
michael@0 972 Bigint *b1;
michael@0 973 ULong *x, *x1, *xe, z;
michael@0 974
michael@0 975 #ifdef Pack_32
michael@0 976 n = k >> 5;
michael@0 977 #else
michael@0 978 n = k >> 4;
michael@0 979 #endif
michael@0 980 k1 = b->k;
michael@0 981 n1 = n + b->wds + 1;
michael@0 982 for(i = b->maxwds; n1 > i; i <<= 1)
michael@0 983 k1++;
michael@0 984 b1 = Balloc(PASS_STATE k1);
michael@0 985 x1 = b1->x;
michael@0 986 for(i = 0; i < n; i++)
michael@0 987 *x1++ = 0;
michael@0 988 x = b->x;
michael@0 989 xe = x + b->wds;
michael@0 990 #ifdef Pack_32
michael@0 991 if (k &= 0x1f) {
michael@0 992 k1 = 32 - k;
michael@0 993 z = 0;
michael@0 994 do {
michael@0 995 *x1++ = *x << k | z;
michael@0 996 z = *x++ >> k1;
michael@0 997 }
michael@0 998 while(x < xe);
michael@0 999 if ((*x1 = z))
michael@0 1000 ++n1;
michael@0 1001 }
michael@0 1002 #else
michael@0 1003 if (k &= 0xf) {
michael@0 1004 k1 = 16 - k;
michael@0 1005 z = 0;
michael@0 1006 do {
michael@0 1007 *x1++ = *x << k & 0xffff | z;
michael@0 1008 z = *x++ >> k1;
michael@0 1009 }
michael@0 1010 while(x < xe);
michael@0 1011 if (*x1 = z)
michael@0 1012 ++n1;
michael@0 1013 }
michael@0 1014 #endif
michael@0 1015 else do
michael@0 1016 *x1++ = *x++;
michael@0 1017 while(x < xe);
michael@0 1018 b1->wds = n1 - 1;
michael@0 1019 Bfree(PASS_STATE b);
michael@0 1020 return b1;
michael@0 1021 }
michael@0 1022
michael@0 1023 static int
michael@0 1024 cmp
michael@0 1025 #ifdef KR_headers
michael@0 1026 (a, b) Bigint *a, *b;
michael@0 1027 #else
michael@0 1028 (Bigint *a, Bigint *b)
michael@0 1029 #endif
michael@0 1030 {
michael@0 1031 ULong *xa, *xa0, *xb, *xb0;
michael@0 1032 int i, j;
michael@0 1033
michael@0 1034 i = a->wds;
michael@0 1035 j = b->wds;
michael@0 1036 #ifdef DEBUG
michael@0 1037 if (i > 1 && !a->x[i-1])
michael@0 1038 Bug("cmp called with a->x[a->wds-1] == 0");
michael@0 1039 if (j > 1 && !b->x[j-1])
michael@0 1040 Bug("cmp called with b->x[b->wds-1] == 0");
michael@0 1041 #endif
michael@0 1042 if (i -= j)
michael@0 1043 return i;
michael@0 1044 xa0 = a->x;
michael@0 1045 xa = xa0 + j;
michael@0 1046 xb0 = b->x;
michael@0 1047 xb = xb0 + j;
michael@0 1048 for(;;) {
michael@0 1049 if (*--xa != *--xb)
michael@0 1050 return *xa < *xb ? -1 : 1;
michael@0 1051 if (xa <= xa0)
michael@0 1052 break;
michael@0 1053 }
michael@0 1054 return 0;
michael@0 1055 }
michael@0 1056
michael@0 1057 static Bigint *
michael@0 1058 diff
michael@0 1059 #ifdef KR_headers
michael@0 1060 (STATE_PARAM a, b) STATE_PARAM_DECL Bigint *a, *b;
michael@0 1061 #else
michael@0 1062 (STATE_PARAM Bigint *a, Bigint *b)
michael@0 1063 #endif
michael@0 1064 {
michael@0 1065 Bigint *c;
michael@0 1066 int i, wa, wb;
michael@0 1067 ULong *xa, *xae, *xb, *xbe, *xc;
michael@0 1068 #ifdef ULLong
michael@0 1069 ULLong borrow, y;
michael@0 1070 #else
michael@0 1071 ULong borrow, y;
michael@0 1072 #ifdef Pack_32
michael@0 1073 ULong z;
michael@0 1074 #endif
michael@0 1075 #endif
michael@0 1076
michael@0 1077 i = cmp(a,b);
michael@0 1078 if (!i) {
michael@0 1079 c = Balloc(PASS_STATE 0);
michael@0 1080 c->wds = 1;
michael@0 1081 c->x[0] = 0;
michael@0 1082 return c;
michael@0 1083 }
michael@0 1084 if (i < 0) {
michael@0 1085 c = a;
michael@0 1086 a = b;
michael@0 1087 b = c;
michael@0 1088 i = 1;
michael@0 1089 }
michael@0 1090 else
michael@0 1091 i = 0;
michael@0 1092 c = Balloc(PASS_STATE a->k);
michael@0 1093 c->sign = i;
michael@0 1094 wa = a->wds;
michael@0 1095 xa = a->x;
michael@0 1096 xae = xa + wa;
michael@0 1097 wb = b->wds;
michael@0 1098 xb = b->x;
michael@0 1099 xbe = xb + wb;
michael@0 1100 xc = c->x;
michael@0 1101 borrow = 0;
michael@0 1102 #ifdef ULLong
michael@0 1103 do {
michael@0 1104 y = (ULLong)*xa++ - *xb++ - borrow;
michael@0 1105 borrow = y >> 32 & (ULong)1;
michael@0 1106 *xc++ = (ULong) y & FFFFFFFF;
michael@0 1107 }
michael@0 1108 while(xb < xbe);
michael@0 1109 while(xa < xae) {
michael@0 1110 y = *xa++ - borrow;
michael@0 1111 borrow = y >> 32 & (ULong)1;
michael@0 1112 *xc++ = (ULong) y & FFFFFFFF;
michael@0 1113 }
michael@0 1114 #else
michael@0 1115 #ifdef Pack_32
michael@0 1116 do {
michael@0 1117 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
michael@0 1118 borrow = (y & 0x10000) >> 16;
michael@0 1119 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
michael@0 1120 borrow = (z & 0x10000) >> 16;
michael@0 1121 Storeinc(xc, z, y);
michael@0 1122 }
michael@0 1123 while(xb < xbe);
michael@0 1124 while(xa < xae) {
michael@0 1125 y = (*xa & 0xffff) - borrow;
michael@0 1126 borrow = (y & 0x10000) >> 16;
michael@0 1127 z = (*xa++ >> 16) - borrow;
michael@0 1128 borrow = (z & 0x10000) >> 16;
michael@0 1129 Storeinc(xc, z, y);
michael@0 1130 }
michael@0 1131 #else
michael@0 1132 do {
michael@0 1133 y = *xa++ - *xb++ - borrow;
michael@0 1134 borrow = (y & 0x10000) >> 16;
michael@0 1135 *xc++ = y & 0xffff;
michael@0 1136 }
michael@0 1137 while(xb < xbe);
michael@0 1138 while(xa < xae) {
michael@0 1139 y = *xa++ - borrow;
michael@0 1140 borrow = (y & 0x10000) >> 16;
michael@0 1141 *xc++ = y & 0xffff;
michael@0 1142 }
michael@0 1143 #endif
michael@0 1144 #endif
michael@0 1145 while(!*--xc)
michael@0 1146 wa--;
michael@0 1147 c->wds = wa;
michael@0 1148 return c;
michael@0 1149 }
michael@0 1150
michael@0 1151 static double
michael@0 1152 ulp
michael@0 1153 #ifdef KR_headers
michael@0 1154 (x) U x;
michael@0 1155 #else
michael@0 1156 (U x)
michael@0 1157 #endif
michael@0 1158 {
michael@0 1159 Long L;
michael@0 1160 U a;
michael@0 1161
michael@0 1162 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
michael@0 1163 #ifndef Avoid_Underflow
michael@0 1164 #ifndef Sudden_Underflow
michael@0 1165 if (L > 0) {
michael@0 1166 #endif
michael@0 1167 #endif
michael@0 1168 #ifdef IBM
michael@0 1169 L |= Exp_msk1 >> 4;
michael@0 1170 #endif
michael@0 1171 word0(a) = L;
michael@0 1172 word1(a) = 0;
michael@0 1173 #ifndef Avoid_Underflow
michael@0 1174 #ifndef Sudden_Underflow
michael@0 1175 }
michael@0 1176 else {
michael@0 1177 L = -L >> Exp_shift;
michael@0 1178 if (L < Exp_shift) {
michael@0 1179 word0(a) = 0x80000 >> L;
michael@0 1180 word1(a) = 0;
michael@0 1181 }
michael@0 1182 else {
michael@0 1183 word0(a) = 0;
michael@0 1184 L -= Exp_shift;
michael@0 1185 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
michael@0 1186 }
michael@0 1187 }
michael@0 1188 #endif
michael@0 1189 #endif
michael@0 1190 return dval(a);
michael@0 1191 }
michael@0 1192
michael@0 1193 static double
michael@0 1194 b2d
michael@0 1195 #ifdef KR_headers
michael@0 1196 (a, e) Bigint *a; int *e;
michael@0 1197 #else
michael@0 1198 (Bigint *a, int *e)
michael@0 1199 #endif
michael@0 1200 {
michael@0 1201 ULong *xa, *xa0, w, y, z;
michael@0 1202 int k;
michael@0 1203 U d;
michael@0 1204 #ifdef VAX
michael@0 1205 ULong d0, d1;
michael@0 1206 #else
michael@0 1207 #define d0 word0(d)
michael@0 1208 #define d1 word1(d)
michael@0 1209 #endif
michael@0 1210
michael@0 1211 xa0 = a->x;
michael@0 1212 xa = xa0 + a->wds;
michael@0 1213 y = *--xa;
michael@0 1214 #ifdef DEBUG
michael@0 1215 if (!y) Bug("zero y in b2d");
michael@0 1216 #endif
michael@0 1217 k = hi0bits(y);
michael@0 1218 *e = 32 - k;
michael@0 1219 #ifdef Pack_32
michael@0 1220 if (k < Ebits) {
michael@0 1221 d0 = Exp_1 | y >> (Ebits - k);
michael@0 1222 w = xa > xa0 ? *--xa : 0;
michael@0 1223 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
michael@0 1224 goto ret_d;
michael@0 1225 }
michael@0 1226 z = xa > xa0 ? *--xa : 0;
michael@0 1227 if (k -= Ebits) {
michael@0 1228 d0 = Exp_1 | y << k | z >> (32 - k);
michael@0 1229 y = xa > xa0 ? *--xa : 0;
michael@0 1230 d1 = z << k | y >> (32 - k);
michael@0 1231 }
michael@0 1232 else {
michael@0 1233 d0 = Exp_1 | y;
michael@0 1234 d1 = z;
michael@0 1235 }
michael@0 1236 #else
michael@0 1237 if (k < Ebits + 16) {
michael@0 1238 z = xa > xa0 ? *--xa : 0;
michael@0 1239 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
michael@0 1240 w = xa > xa0 ? *--xa : 0;
michael@0 1241 y = xa > xa0 ? *--xa : 0;
michael@0 1242 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
michael@0 1243 goto ret_d;
michael@0 1244 }
michael@0 1245 z = xa > xa0 ? *--xa : 0;
michael@0 1246 w = xa > xa0 ? *--xa : 0;
michael@0 1247 k -= Ebits + 16;
michael@0 1248 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
michael@0 1249 y = xa > xa0 ? *--xa : 0;
michael@0 1250 d1 = w << k + 16 | y << k;
michael@0 1251 #endif
michael@0 1252 ret_d:
michael@0 1253 #ifdef VAX
michael@0 1254 word0(d) = d0 >> 16 | d0 << 16;
michael@0 1255 word1(d) = d1 >> 16 | d1 << 16;
michael@0 1256 #else
michael@0 1257 #undef d0
michael@0 1258 #undef d1
michael@0 1259 #endif
michael@0 1260 return dval(d);
michael@0 1261 }
michael@0 1262
michael@0 1263 static Bigint *
michael@0 1264 d2b
michael@0 1265 #ifdef KR_headers
michael@0 1266 (STATE_PARAM d, e, bits) STATE_PARAM_DECL U d; int *e, *bits;
michael@0 1267 #else
michael@0 1268 (STATE_PARAM U d, int *e, int *bits)
michael@0 1269 #endif
michael@0 1270 {
michael@0 1271 Bigint *b;
michael@0 1272 int de, k;
michael@0 1273 ULong *x, y, z;
michael@0 1274 #ifndef Sudden_Underflow
michael@0 1275 int i;
michael@0 1276 #endif
michael@0 1277 #ifdef VAX
michael@0 1278 ULong d0, d1;
michael@0 1279 d0 = word0(d) >> 16 | word0(d) << 16;
michael@0 1280 d1 = word1(d) >> 16 | word1(d) << 16;
michael@0 1281 #else
michael@0 1282 #define d0 word0(d)
michael@0 1283 #define d1 word1(d)
michael@0 1284 #endif
michael@0 1285
michael@0 1286 #ifdef Pack_32
michael@0 1287 b = Balloc(PASS_STATE 1);
michael@0 1288 #else
michael@0 1289 b = Balloc(PASS_STATE 2);
michael@0 1290 #endif
michael@0 1291 x = b->x;
michael@0 1292
michael@0 1293 z = d0 & Frac_mask;
michael@0 1294 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
michael@0 1295 #ifdef Sudden_Underflow
michael@0 1296 de = (int)(d0 >> Exp_shift);
michael@0 1297 #ifndef IBM
michael@0 1298 z |= Exp_msk11;
michael@0 1299 #endif
michael@0 1300 #else
michael@0 1301 if ((de = (int)(d0 >> Exp_shift)))
michael@0 1302 z |= Exp_msk1;
michael@0 1303 #endif
michael@0 1304 #ifdef Pack_32
michael@0 1305 if ((y = d1)) {
michael@0 1306 if ((k = lo0bits(&y))) {
michael@0 1307 x[0] = y | z << (32 - k);
michael@0 1308 z >>= k;
michael@0 1309 }
michael@0 1310 else
michael@0 1311 x[0] = y;
michael@0 1312 #ifndef Sudden_Underflow
michael@0 1313 i =
michael@0 1314 #endif
michael@0 1315 b->wds = (x[1] = z) ? 2 : 1;
michael@0 1316 }
michael@0 1317 else {
michael@0 1318 k = lo0bits(&z);
michael@0 1319 x[0] = z;
michael@0 1320 #ifndef Sudden_Underflow
michael@0 1321 i =
michael@0 1322 #endif
michael@0 1323 b->wds = 1;
michael@0 1324 k += 32;
michael@0 1325 }
michael@0 1326 #else
michael@0 1327 if (y = d1) {
michael@0 1328 if (k = lo0bits(&y))
michael@0 1329 if (k >= 16) {
michael@0 1330 x[0] = y | z << 32 - k & 0xffff;
michael@0 1331 x[1] = z >> k - 16 & 0xffff;
michael@0 1332 x[2] = z >> k;
michael@0 1333 i = 2;
michael@0 1334 }
michael@0 1335 else {
michael@0 1336 x[0] = y & 0xffff;
michael@0 1337 x[1] = y >> 16 | z << 16 - k & 0xffff;
michael@0 1338 x[2] = z >> k & 0xffff;
michael@0 1339 x[3] = z >> k+16;
michael@0 1340 i = 3;
michael@0 1341 }
michael@0 1342 else {
michael@0 1343 x[0] = y & 0xffff;
michael@0 1344 x[1] = y >> 16;
michael@0 1345 x[2] = z & 0xffff;
michael@0 1346 x[3] = z >> 16;
michael@0 1347 i = 3;
michael@0 1348 }
michael@0 1349 }
michael@0 1350 else {
michael@0 1351 #ifdef DEBUG
michael@0 1352 if (!z)
michael@0 1353 Bug("Zero passed to d2b");
michael@0 1354 #endif
michael@0 1355 k = lo0bits(&z);
michael@0 1356 if (k >= 16) {
michael@0 1357 x[0] = z;
michael@0 1358 i = 0;
michael@0 1359 }
michael@0 1360 else {
michael@0 1361 x[0] = z & 0xffff;
michael@0 1362 x[1] = z >> 16;
michael@0 1363 i = 1;
michael@0 1364 }
michael@0 1365 k += 32;
michael@0 1366 }
michael@0 1367 while(!x[i])
michael@0 1368 --i;
michael@0 1369 b->wds = i + 1;
michael@0 1370 #endif
michael@0 1371 #ifndef Sudden_Underflow
michael@0 1372 if (de) {
michael@0 1373 #endif
michael@0 1374 #ifdef IBM
michael@0 1375 *e = (de - Bias - (P-1) << 2) + k;
michael@0 1376 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
michael@0 1377 #else
michael@0 1378 *e = de - Bias - (P-1) + k;
michael@0 1379 *bits = P - k;
michael@0 1380 #endif
michael@0 1381 #ifndef Sudden_Underflow
michael@0 1382 }
michael@0 1383 else {
michael@0 1384 *e = de - Bias - (P-1) + 1 + k;
michael@0 1385 #ifdef Pack_32
michael@0 1386 *bits = 32*i - hi0bits(x[i-1]);
michael@0 1387 #else
michael@0 1388 *bits = (i+2)*16 - hi0bits(x[i]);
michael@0 1389 #endif
michael@0 1390 }
michael@0 1391 #endif
michael@0 1392 return b;
michael@0 1393 }
michael@0 1394 #undef d0
michael@0 1395 #undef d1
michael@0 1396
michael@0 1397 static double
michael@0 1398 ratio
michael@0 1399 #ifdef KR_headers
michael@0 1400 (a, b) Bigint *a, *b;
michael@0 1401 #else
michael@0 1402 (Bigint *a, Bigint *b)
michael@0 1403 #endif
michael@0 1404 {
michael@0 1405 U da, db;
michael@0 1406 int k, ka, kb;
michael@0 1407
michael@0 1408 dval(da) = b2d(a, &ka);
michael@0 1409 dval(db) = b2d(b, &kb);
michael@0 1410 #ifdef Pack_32
michael@0 1411 k = ka - kb + 32*(a->wds - b->wds);
michael@0 1412 #else
michael@0 1413 k = ka - kb + 16*(a->wds - b->wds);
michael@0 1414 #endif
michael@0 1415 #ifdef IBM
michael@0 1416 if (k > 0) {
michael@0 1417 word0(da) += (k >> 2)*Exp_msk1;
michael@0 1418 if (k &= 3)
michael@0 1419 dval(da) *= 1 << k;
michael@0 1420 }
michael@0 1421 else {
michael@0 1422 k = -k;
michael@0 1423 word0(db) += (k >> 2)*Exp_msk1;
michael@0 1424 if (k &= 3)
michael@0 1425 dval(db) *= 1 << k;
michael@0 1426 }
michael@0 1427 #else
michael@0 1428 if (k > 0)
michael@0 1429 word0(da) += k*Exp_msk1;
michael@0 1430 else {
michael@0 1431 k = -k;
michael@0 1432 word0(db) += k*Exp_msk1;
michael@0 1433 }
michael@0 1434 #endif
michael@0 1435 return dval(da) / dval(db);
michael@0 1436 }
michael@0 1437
michael@0 1438 static CONST double
michael@0 1439 tens[] = {
michael@0 1440 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
michael@0 1441 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
michael@0 1442 1e20, 1e21, 1e22
michael@0 1443 #ifdef VAX
michael@0 1444 , 1e23, 1e24
michael@0 1445 #endif
michael@0 1446 };
michael@0 1447
michael@0 1448 static CONST double
michael@0 1449 #ifdef IEEE_Arith
michael@0 1450 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
michael@0 1451 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
michael@0 1452 #ifdef Avoid_Underflow
michael@0 1453 9007199254740992.*9007199254740992.e-256
michael@0 1454 /* = 2^106 * 1e-53 */
michael@0 1455 #else
michael@0 1456 1e-256
michael@0 1457 #endif
michael@0 1458 };
michael@0 1459 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
michael@0 1460 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
michael@0 1461 #define Scale_Bit 0x10
michael@0 1462 #define n_bigtens 5
michael@0 1463 #else
michael@0 1464 #ifdef IBM
michael@0 1465 bigtens[] = { 1e16, 1e32, 1e64 };
michael@0 1466 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
michael@0 1467 #define n_bigtens 3
michael@0 1468 #else
michael@0 1469 bigtens[] = { 1e16, 1e32 };
michael@0 1470 static CONST double tinytens[] = { 1e-16, 1e-32 };
michael@0 1471 #define n_bigtens 2
michael@0 1472 #endif
michael@0 1473 #endif
michael@0 1474
michael@0 1475 static double
michael@0 1476 _strtod
michael@0 1477 #ifdef KR_headers
michael@0 1478 (STATE_PARAM s00, se) STATE_PARAM_DECL CONST char *s00; char **se;
michael@0 1479 #else
michael@0 1480 (STATE_PARAM CONST char *s00, char **se)
michael@0 1481 #endif
michael@0 1482 {
michael@0 1483 #ifdef Avoid_Underflow
michael@0 1484 int scale;
michael@0 1485 #endif
michael@0 1486 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
michael@0 1487 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
michael@0 1488 CONST char *s, *s0, *s1;
michael@0 1489 double aadj, adj;
michael@0 1490 U aadj1, rv, rv0;
michael@0 1491 Long L;
michael@0 1492 ULong y, z;
michael@0 1493 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
michael@0 1494 #ifdef SET_INEXACT
michael@0 1495 int inexact, oldinexact;
michael@0 1496 #endif
michael@0 1497 #ifdef Honor_FLT_ROUNDS
michael@0 1498 int rounding;
michael@0 1499 #endif
michael@0 1500 #ifdef USE_LOCALE
michael@0 1501 CONST char *s2;
michael@0 1502 #endif
michael@0 1503
michael@0 1504 #ifdef __GNUC__
michael@0 1505 delta = bb = bd = bs = 0;
michael@0 1506 #endif
michael@0 1507
michael@0 1508 sign = nz0 = nz = 0;
michael@0 1509 dval(rv) = 0.;
michael@0 1510 for(s = s00;;s++) switch(*s) {
michael@0 1511 case '-':
michael@0 1512 sign = 1;
michael@0 1513 /* no break */
michael@0 1514 case '+':
michael@0 1515 if (*++s)
michael@0 1516 goto break2;
michael@0 1517 /* no break */
michael@0 1518 case 0:
michael@0 1519 goto ret0;
michael@0 1520 case '\t':
michael@0 1521 case '\n':
michael@0 1522 case '\v':
michael@0 1523 case '\f':
michael@0 1524 case '\r':
michael@0 1525 case ' ':
michael@0 1526 continue;
michael@0 1527 default:
michael@0 1528 goto break2;
michael@0 1529 }
michael@0 1530 break2:
michael@0 1531 if (*s == '0') {
michael@0 1532 nz0 = 1;
michael@0 1533 while(*++s == '0') ;
michael@0 1534 if (!*s)
michael@0 1535 goto ret;
michael@0 1536 }
michael@0 1537 s0 = s;
michael@0 1538 y = z = 0;
michael@0 1539 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
michael@0 1540 if (nd < 9)
michael@0 1541 y = 10*y + c - '0';
michael@0 1542 else if (nd < 16)
michael@0 1543 z = 10*z + c - '0';
michael@0 1544 nd0 = nd;
michael@0 1545 #ifdef USE_LOCALE
michael@0 1546 s1 = localeconv()->decimal_point;
michael@0 1547 if (c == *s1) {
michael@0 1548 c = '.';
michael@0 1549 if (*++s1) {
michael@0 1550 s2 = s;
michael@0 1551 for(;;) {
michael@0 1552 if (*++s2 != *s1) {
michael@0 1553 c = 0;
michael@0 1554 break;
michael@0 1555 }
michael@0 1556 if (!*++s1) {
michael@0 1557 s = s2;
michael@0 1558 break;
michael@0 1559 }
michael@0 1560 }
michael@0 1561 }
michael@0 1562 }
michael@0 1563 #endif
michael@0 1564 if (c == '.') {
michael@0 1565 c = *++s;
michael@0 1566 if (!nd) {
michael@0 1567 for(; c == '0'; c = *++s)
michael@0 1568 nz++;
michael@0 1569 if (c > '0' && c <= '9') {
michael@0 1570 s0 = s;
michael@0 1571 nf += nz;
michael@0 1572 nz = 0;
michael@0 1573 goto have_dig;
michael@0 1574 }
michael@0 1575 goto dig_done;
michael@0 1576 }
michael@0 1577 for(; c >= '0' && c <= '9'; c = *++s) {
michael@0 1578 have_dig:
michael@0 1579 nz++;
michael@0 1580 if (c -= '0') {
michael@0 1581 nf += nz;
michael@0 1582 for(i = 1; i < nz; i++)
michael@0 1583 if (nd++ < 9)
michael@0 1584 y *= 10;
michael@0 1585 else if (nd <= DBL_DIG + 1)
michael@0 1586 z *= 10;
michael@0 1587 if (nd++ < 9)
michael@0 1588 y = 10*y + c;
michael@0 1589 else if (nd <= DBL_DIG + 1)
michael@0 1590 z = 10*z + c;
michael@0 1591 nz = 0;
michael@0 1592 }
michael@0 1593 }
michael@0 1594 }
michael@0 1595 dig_done:
michael@0 1596 e = 0;
michael@0 1597 if (c == 'e' || c == 'E') {
michael@0 1598 if (!nd && !nz && !nz0) {
michael@0 1599 goto ret0;
michael@0 1600 }
michael@0 1601 s00 = s;
michael@0 1602 esign = 0;
michael@0 1603 switch(c = *++s) {
michael@0 1604 case '-':
michael@0 1605 esign = 1;
michael@0 1606 case '+':
michael@0 1607 c = *++s;
michael@0 1608 }
michael@0 1609 if (c >= '0' && c <= '9') {
michael@0 1610 while(c == '0')
michael@0 1611 c = *++s;
michael@0 1612 if (c > '0' && c <= '9') {
michael@0 1613 L = c - '0';
michael@0 1614 s1 = s;
michael@0 1615 while((c = *++s) >= '0' && c <= '9')
michael@0 1616 L = 10*L + c - '0';
michael@0 1617 if (s - s1 > 8 || L > 19999)
michael@0 1618 /* Avoid confusion from exponents
michael@0 1619 * so large that e might overflow.
michael@0 1620 */
michael@0 1621 e = 19999; /* safe for 16 bit ints */
michael@0 1622 else
michael@0 1623 e = (int)L;
michael@0 1624 if (esign)
michael@0 1625 e = -e;
michael@0 1626 }
michael@0 1627 else
michael@0 1628 e = 0;
michael@0 1629 }
michael@0 1630 else
michael@0 1631 s = s00;
michael@0 1632 }
michael@0 1633 if (!nd) {
michael@0 1634 if (!nz && !nz0) {
michael@0 1635 ret0:
michael@0 1636 s = s00;
michael@0 1637 sign = 0;
michael@0 1638 }
michael@0 1639 goto ret;
michael@0 1640 }
michael@0 1641 e1 = e -= nf;
michael@0 1642
michael@0 1643 /* Now we have nd0 digits, starting at s0, followed by a
michael@0 1644 * decimal point, followed by nd-nd0 digits. The number we're
michael@0 1645 * after is the integer represented by those digits times
michael@0 1646 * 10**e */
michael@0 1647
michael@0 1648 if (!nd0)
michael@0 1649 nd0 = nd;
michael@0 1650 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
michael@0 1651 dval(rv) = y;
michael@0 1652 if (k > 9) {
michael@0 1653 #ifdef SET_INEXACT
michael@0 1654 if (k > DBL_DIG)
michael@0 1655 oldinexact = get_inexact();
michael@0 1656 #endif
michael@0 1657 dval(rv) = tens[k - 9] * dval(rv) + z;
michael@0 1658 }
michael@0 1659 bd0 = 0;
michael@0 1660 if (nd <= DBL_DIG
michael@0 1661 #ifndef RND_PRODQUOT
michael@0 1662 #ifndef Honor_FLT_ROUNDS
michael@0 1663 && Flt_Rounds == 1
michael@0 1664 #endif
michael@0 1665 #endif
michael@0 1666 ) {
michael@0 1667 if (!e)
michael@0 1668 goto ret;
michael@0 1669 if (e > 0) {
michael@0 1670 if (e <= Ten_pmax) {
michael@0 1671 #ifdef VAX
michael@0 1672 goto vax_ovfl_check;
michael@0 1673 #else
michael@0 1674 #ifdef Honor_FLT_ROUNDS
michael@0 1675 /* round correctly FLT_ROUNDS = 2 or 3 */
michael@0 1676 if (sign) {
michael@0 1677 rv = -rv;
michael@0 1678 sign = 0;
michael@0 1679 }
michael@0 1680 #endif
michael@0 1681 /* rv = */ rounded_product(dval(rv), tens[e]);
michael@0 1682 goto ret;
michael@0 1683 #endif
michael@0 1684 }
michael@0 1685 i = DBL_DIG - nd;
michael@0 1686 if (e <= Ten_pmax + i) {
michael@0 1687 /* A fancier test would sometimes let us do
michael@0 1688 * this for larger i values.
michael@0 1689 */
michael@0 1690 #ifdef Honor_FLT_ROUNDS
michael@0 1691 /* round correctly FLT_ROUNDS = 2 or 3 */
michael@0 1692 if (sign) {
michael@0 1693 rv = -rv;
michael@0 1694 sign = 0;
michael@0 1695 }
michael@0 1696 #endif
michael@0 1697 e -= i;
michael@0 1698 dval(rv) *= tens[i];
michael@0 1699 #ifdef VAX
michael@0 1700 /* VAX exponent range is so narrow we must
michael@0 1701 * worry about overflow here...
michael@0 1702 */
michael@0 1703 vax_ovfl_check:
michael@0 1704 word0(rv) -= P*Exp_msk1;
michael@0 1705 /* rv = */ rounded_product(dval(rv), tens[e]);
michael@0 1706 if ((word0(rv) & Exp_mask)
michael@0 1707 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
michael@0 1708 goto ovfl;
michael@0 1709 word0(rv) += P*Exp_msk1;
michael@0 1710 #else
michael@0 1711 /* rv = */ rounded_product(dval(rv), tens[e]);
michael@0 1712 #endif
michael@0 1713 goto ret;
michael@0 1714 }
michael@0 1715 }
michael@0 1716 #ifndef Inaccurate_Divide
michael@0 1717 else if (e >= -Ten_pmax) {
michael@0 1718 #ifdef Honor_FLT_ROUNDS
michael@0 1719 /* round correctly FLT_ROUNDS = 2 or 3 */
michael@0 1720 if (sign) {
michael@0 1721 rv = -rv;
michael@0 1722 sign = 0;
michael@0 1723 }
michael@0 1724 #endif
michael@0 1725 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
michael@0 1726 goto ret;
michael@0 1727 }
michael@0 1728 #endif
michael@0 1729 }
michael@0 1730 e1 += nd - k;
michael@0 1731
michael@0 1732 #ifdef IEEE_Arith
michael@0 1733 #ifdef SET_INEXACT
michael@0 1734 inexact = 1;
michael@0 1735 if (k <= DBL_DIG)
michael@0 1736 oldinexact = get_inexact();
michael@0 1737 #endif
michael@0 1738 #ifdef Avoid_Underflow
michael@0 1739 scale = 0;
michael@0 1740 #endif
michael@0 1741 #ifdef Honor_FLT_ROUNDS
michael@0 1742 if ((rounding = Flt_Rounds) >= 2) {
michael@0 1743 if (sign)
michael@0 1744 rounding = rounding == 2 ? 0 : 2;
michael@0 1745 else
michael@0 1746 if (rounding != 2)
michael@0 1747 rounding = 0;
michael@0 1748 }
michael@0 1749 #endif
michael@0 1750 #endif /*IEEE_Arith*/
michael@0 1751
michael@0 1752 /* Get starting approximation = rv * 10**e1 */
michael@0 1753
michael@0 1754 if (e1 > 0) {
michael@0 1755 if ((i = e1 & 15))
michael@0 1756 dval(rv) *= tens[i];
michael@0 1757 if (e1 &= ~15) {
michael@0 1758 if (e1 > DBL_MAX_10_EXP) {
michael@0 1759 ovfl:
michael@0 1760 #ifndef NO_ERRNO
michael@0 1761 errno = ERANGE;
michael@0 1762 #endif
michael@0 1763 /* Can't trust HUGE_VAL */
michael@0 1764 #ifdef IEEE_Arith
michael@0 1765 #ifdef Honor_FLT_ROUNDS
michael@0 1766 switch(rounding) {
michael@0 1767 case 0: /* toward 0 */
michael@0 1768 case 3: /* toward -infinity */
michael@0 1769 word0(rv) = Big0;
michael@0 1770 word1(rv) = Big1;
michael@0 1771 break;
michael@0 1772 default:
michael@0 1773 word0(rv) = Exp_mask;
michael@0 1774 word1(rv) = 0;
michael@0 1775 }
michael@0 1776 #else /*Honor_FLT_ROUNDS*/
michael@0 1777 word0(rv) = Exp_mask;
michael@0 1778 word1(rv) = 0;
michael@0 1779 #endif /*Honor_FLT_ROUNDS*/
michael@0 1780 #ifdef SET_INEXACT
michael@0 1781 /* set overflow bit */
michael@0 1782 dval(rv0) = 1e300;
michael@0 1783 dval(rv0) *= dval(rv0);
michael@0 1784 #endif
michael@0 1785 #else /*IEEE_Arith*/
michael@0 1786 word0(rv) = Big0;
michael@0 1787 word1(rv) = Big1;
michael@0 1788 #endif /*IEEE_Arith*/
michael@0 1789 if (bd0)
michael@0 1790 goto retfree;
michael@0 1791 goto ret;
michael@0 1792 }
michael@0 1793 e1 >>= 4;
michael@0 1794 for(j = 0; e1 > 1; j++, e1 >>= 1)
michael@0 1795 if (e1 & 1)
michael@0 1796 dval(rv) *= bigtens[j];
michael@0 1797 /* The last multiplication could overflow. */
michael@0 1798 word0(rv) -= P*Exp_msk1;
michael@0 1799 dval(rv) *= bigtens[j];
michael@0 1800 if ((z = word0(rv) & Exp_mask)
michael@0 1801 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
michael@0 1802 goto ovfl;
michael@0 1803 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
michael@0 1804 /* set to largest number */
michael@0 1805 /* (Can't trust DBL_MAX) */
michael@0 1806 word0(rv) = Big0;
michael@0 1807 word1(rv) = Big1;
michael@0 1808 }
michael@0 1809 else
michael@0 1810 word0(rv) += P*Exp_msk1;
michael@0 1811 }
michael@0 1812 }
michael@0 1813 else if (e1 < 0) {
michael@0 1814 e1 = -e1;
michael@0 1815 if ((i = e1 & 15))
michael@0 1816 dval(rv) /= tens[i];
michael@0 1817 if (e1 >>= 4) {
michael@0 1818 if (e1 >= 1 << n_bigtens)
michael@0 1819 goto undfl;
michael@0 1820 #ifdef Avoid_Underflow
michael@0 1821 if (e1 & Scale_Bit)
michael@0 1822 scale = 2*P;
michael@0 1823 for(j = 0; e1 > 0; j++, e1 >>= 1)
michael@0 1824 if (e1 & 1)
michael@0 1825 dval(rv) *= tinytens[j];
michael@0 1826 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
michael@0 1827 >> Exp_shift)) > 0) {
michael@0 1828 /* scaled rv is denormal; zap j low bits */
michael@0 1829 if (j >= 32) {
michael@0 1830 word1(rv) = 0;
michael@0 1831 if (j >= 53)
michael@0 1832 word0(rv) = (P+2)*Exp_msk1;
michael@0 1833 else
michael@0 1834 word0(rv) &= 0xffffffff << (j-32);
michael@0 1835 }
michael@0 1836 else
michael@0 1837 word1(rv) &= 0xffffffff << j;
michael@0 1838 }
michael@0 1839 #else
michael@0 1840 for(j = 0; e1 > 1; j++, e1 >>= 1)
michael@0 1841 if (e1 & 1)
michael@0 1842 dval(rv) *= tinytens[j];
michael@0 1843 /* The last multiplication could underflow. */
michael@0 1844 dval(rv0) = dval(rv);
michael@0 1845 dval(rv) *= tinytens[j];
michael@0 1846 if (!dval(rv)) {
michael@0 1847 dval(rv) = 2.*dval(rv0);
michael@0 1848 dval(rv) *= tinytens[j];
michael@0 1849 #endif
michael@0 1850 if (!dval(rv)) {
michael@0 1851 undfl:
michael@0 1852 dval(rv) = 0.;
michael@0 1853 #ifndef NO_ERRNO
michael@0 1854 errno = ERANGE;
michael@0 1855 #endif
michael@0 1856 if (bd0)
michael@0 1857 goto retfree;
michael@0 1858 goto ret;
michael@0 1859 }
michael@0 1860 #ifndef Avoid_Underflow
michael@0 1861 word0(rv) = Tiny0;
michael@0 1862 word1(rv) = Tiny1;
michael@0 1863 /* The refinement below will clean
michael@0 1864 * this approximation up.
michael@0 1865 */
michael@0 1866 }
michael@0 1867 #endif
michael@0 1868 }
michael@0 1869 }
michael@0 1870
michael@0 1871 /* Now the hard part -- adjusting rv to the correct value.*/
michael@0 1872
michael@0 1873 /* Put digits into bd: true value = bd * 10^e */
michael@0 1874
michael@0 1875 bd0 = s2b(PASS_STATE s0, nd0, nd, y);
michael@0 1876
michael@0 1877 for(;;) {
michael@0 1878 bd = Balloc(PASS_STATE bd0->k);
michael@0 1879 Bcopy(bd, bd0);
michael@0 1880 bb = d2b(PASS_STATE rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
michael@0 1881 bs = i2b(PASS_STATE 1);
michael@0 1882
michael@0 1883 if (e >= 0) {
michael@0 1884 bb2 = bb5 = 0;
michael@0 1885 bd2 = bd5 = e;
michael@0 1886 }
michael@0 1887 else {
michael@0 1888 bb2 = bb5 = -e;
michael@0 1889 bd2 = bd5 = 0;
michael@0 1890 }
michael@0 1891 if (bbe >= 0)
michael@0 1892 bb2 += bbe;
michael@0 1893 else
michael@0 1894 bd2 -= bbe;
michael@0 1895 bs2 = bb2;
michael@0 1896 #ifdef Honor_FLT_ROUNDS
michael@0 1897 if (rounding != 1)
michael@0 1898 bs2++;
michael@0 1899 #endif
michael@0 1900 #ifdef Avoid_Underflow
michael@0 1901 j = bbe - scale;
michael@0 1902 i = j + bbbits - 1; /* logb(rv) */
michael@0 1903 if (i < Emin) /* denormal */
michael@0 1904 j += P - Emin;
michael@0 1905 else
michael@0 1906 j = P + 1 - bbbits;
michael@0 1907 #else /*Avoid_Underflow*/
michael@0 1908 #ifdef Sudden_Underflow
michael@0 1909 #ifdef IBM
michael@0 1910 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
michael@0 1911 #else
michael@0 1912 j = P + 1 - bbbits;
michael@0 1913 #endif
michael@0 1914 #else /*Sudden_Underflow*/
michael@0 1915 j = bbe;
michael@0 1916 i = j + bbbits - 1; /* logb(rv) */
michael@0 1917 if (i < Emin) /* denormal */
michael@0 1918 j += P - Emin;
michael@0 1919 else
michael@0 1920 j = P + 1 - bbbits;
michael@0 1921 #endif /*Sudden_Underflow*/
michael@0 1922 #endif /*Avoid_Underflow*/
michael@0 1923 bb2 += j;
michael@0 1924 bd2 += j;
michael@0 1925 #ifdef Avoid_Underflow
michael@0 1926 bd2 += scale;
michael@0 1927 #endif
michael@0 1928 i = bb2 < bd2 ? bb2 : bd2;
michael@0 1929 if (i > bs2)
michael@0 1930 i = bs2;
michael@0 1931 if (i > 0) {
michael@0 1932 bb2 -= i;
michael@0 1933 bd2 -= i;
michael@0 1934 bs2 -= i;
michael@0 1935 }
michael@0 1936 if (bb5 > 0) {
michael@0 1937 bs = pow5mult(PASS_STATE bs, bb5);
michael@0 1938 bb1 = mult(PASS_STATE bs, bb);
michael@0 1939 Bfree(PASS_STATE bb);
michael@0 1940 bb = bb1;
michael@0 1941 }
michael@0 1942 if (bb2 > 0)
michael@0 1943 bb = lshift(PASS_STATE bb, bb2);
michael@0 1944 if (bd5 > 0)
michael@0 1945 bd = pow5mult(PASS_STATE bd, bd5);
michael@0 1946 if (bd2 > 0)
michael@0 1947 bd = lshift(PASS_STATE bd, bd2);
michael@0 1948 if (bs2 > 0)
michael@0 1949 bs = lshift(PASS_STATE bs, bs2);
michael@0 1950 delta = diff(PASS_STATE bb, bd);
michael@0 1951 dsign = delta->sign;
michael@0 1952 delta->sign = 0;
michael@0 1953 i = cmp(delta, bs);
michael@0 1954 #ifdef Honor_FLT_ROUNDS
michael@0 1955 if (rounding != 1) {
michael@0 1956 if (i < 0) {
michael@0 1957 /* Error is less than an ulp */
michael@0 1958 if (!delta->x[0] && delta->wds <= 1) {
michael@0 1959 /* exact */
michael@0 1960 #ifdef SET_INEXACT
michael@0 1961 inexact = 0;
michael@0 1962 #endif
michael@0 1963 break;
michael@0 1964 }
michael@0 1965 if (rounding) {
michael@0 1966 if (dsign) {
michael@0 1967 adj = 1.;
michael@0 1968 goto apply_adj;
michael@0 1969 }
michael@0 1970 }
michael@0 1971 else if (!dsign) {
michael@0 1972 adj = -1.;
michael@0 1973 if (!word1(rv)
michael@0 1974 && !(word0(rv) & Frac_mask)) {
michael@0 1975 y = word0(rv) & Exp_mask;
michael@0 1976 #ifdef Avoid_Underflow
michael@0 1977 if (!scale || y > 2*P*Exp_msk1)
michael@0 1978 #else
michael@0 1979 if (y)
michael@0 1980 #endif
michael@0 1981 {
michael@0 1982 delta = lshift(PASS_STATE delta,Log2P);
michael@0 1983 if (cmp(delta, bs) <= 0)
michael@0 1984 adj = -0.5;
michael@0 1985 }
michael@0 1986 }
michael@0 1987 apply_adj:
michael@0 1988 #ifdef Avoid_Underflow
michael@0 1989 if (scale && (y = word0(rv) & Exp_mask)
michael@0 1990 <= 2*P*Exp_msk1)
michael@0 1991 word0(adj) += (2*P+1)*Exp_msk1 - y;
michael@0 1992 #else
michael@0 1993 #ifdef Sudden_Underflow
michael@0 1994 if ((word0(rv) & Exp_mask) <=
michael@0 1995 P*Exp_msk1) {
michael@0 1996 word0(rv) += P*Exp_msk1;
michael@0 1997 dval(rv) += adj*ulp(rv);
michael@0 1998 word0(rv) -= P*Exp_msk1;
michael@0 1999 }
michael@0 2000 else
michael@0 2001 #endif /*Sudden_Underflow*/
michael@0 2002 #endif /*Avoid_Underflow*/
michael@0 2003 dval(rv) += adj*ulp(rv);
michael@0 2004 }
michael@0 2005 break;
michael@0 2006 }
michael@0 2007 adj = ratio(delta, bs);
michael@0 2008 if (adj < 1.)
michael@0 2009 adj = 1.;
michael@0 2010 if (adj <= 0x7ffffffe) {
michael@0 2011 /* adj = rounding ? ceil(adj) : floor(adj); */
michael@0 2012 y = adj;
michael@0 2013 if (y != adj) {
michael@0 2014 if (!((rounding>>1) ^ dsign))
michael@0 2015 y++;
michael@0 2016 adj = y;
michael@0 2017 }
michael@0 2018 }
michael@0 2019 #ifdef Avoid_Underflow
michael@0 2020 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
michael@0 2021 word0(adj) += (2*P+1)*Exp_msk1 - y;
michael@0 2022 #else
michael@0 2023 #ifdef Sudden_Underflow
michael@0 2024 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
michael@0 2025 word0(rv) += P*Exp_msk1;
michael@0 2026 adj *= ulp(rv);
michael@0 2027 if (dsign)
michael@0 2028 dval(rv) += adj;
michael@0 2029 else
michael@0 2030 dval(rv) -= adj;
michael@0 2031 word0(rv) -= P*Exp_msk1;
michael@0 2032 goto cont;
michael@0 2033 }
michael@0 2034 #endif /*Sudden_Underflow*/
michael@0 2035 #endif /*Avoid_Underflow*/
michael@0 2036 adj *= ulp(rv);
michael@0 2037 if (dsign)
michael@0 2038 dval(rv) += adj;
michael@0 2039 else
michael@0 2040 dval(rv) -= adj;
michael@0 2041 goto cont;
michael@0 2042 }
michael@0 2043 #endif /*Honor_FLT_ROUNDS*/
michael@0 2044
michael@0 2045 if (i < 0) {
michael@0 2046 /* Error is less than half an ulp -- check for
michael@0 2047 * special case of mantissa a power of two.
michael@0 2048 */
michael@0 2049 if (dsign || word1(rv) || word0(rv) & Bndry_mask
michael@0 2050 #ifdef IEEE_Arith
michael@0 2051 #ifdef Avoid_Underflow
michael@0 2052 || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
michael@0 2053 #else
michael@0 2054 || (word0(rv) & Exp_mask) <= Exp_msk1
michael@0 2055 #endif
michael@0 2056 #endif
michael@0 2057 ) {
michael@0 2058 #ifdef SET_INEXACT
michael@0 2059 if (!delta->x[0] && delta->wds <= 1)
michael@0 2060 inexact = 0;
michael@0 2061 #endif
michael@0 2062 break;
michael@0 2063 }
michael@0 2064 if (!delta->x[0] && delta->wds <= 1) {
michael@0 2065 /* exact result */
michael@0 2066 #ifdef SET_INEXACT
michael@0 2067 inexact = 0;
michael@0 2068 #endif
michael@0 2069 break;
michael@0 2070 }
michael@0 2071 delta = lshift(PASS_STATE delta,Log2P);
michael@0 2072 if (cmp(delta, bs) > 0)
michael@0 2073 goto drop_down;
michael@0 2074 break;
michael@0 2075 }
michael@0 2076 if (i == 0) {
michael@0 2077 /* exactly half-way between */
michael@0 2078 if (dsign) {
michael@0 2079 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
michael@0 2080 && word1(rv) == (
michael@0 2081 #ifdef Avoid_Underflow
michael@0 2082 (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
michael@0 2083 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
michael@0 2084 #endif
michael@0 2085 0xffffffff)) {
michael@0 2086 /*boundary case -- increment exponent*/
michael@0 2087 word0(rv) = (word0(rv) & Exp_mask)
michael@0 2088 + Exp_msk1
michael@0 2089 #ifdef IBM
michael@0 2090 | Exp_msk1 >> 4
michael@0 2091 #endif
michael@0 2092 ;
michael@0 2093 word1(rv) = 0;
michael@0 2094 #ifdef Avoid_Underflow
michael@0 2095 dsign = 0;
michael@0 2096 #endif
michael@0 2097 break;
michael@0 2098 }
michael@0 2099 }
michael@0 2100 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
michael@0 2101 drop_down:
michael@0 2102 /* boundary case -- decrement exponent */
michael@0 2103 #ifdef Sudden_Underflow /*{{*/
michael@0 2104 L = word0(rv) & Exp_mask;
michael@0 2105 #ifdef IBM
michael@0 2106 if (L < Exp_msk1)
michael@0 2107 #else
michael@0 2108 #ifdef Avoid_Underflow
michael@0 2109 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
michael@0 2110 #else
michael@0 2111 if (L <= Exp_msk1)
michael@0 2112 #endif /*Avoid_Underflow*/
michael@0 2113 #endif /*IBM*/
michael@0 2114 goto undfl;
michael@0 2115 L -= Exp_msk1;
michael@0 2116 #else /*Sudden_Underflow}{*/
michael@0 2117 #ifdef Avoid_Underflow
michael@0 2118 if (scale) {
michael@0 2119 L = word0(rv) & Exp_mask;
michael@0 2120 if (L <= (2*P+1)*Exp_msk1) {
michael@0 2121 if (L > (P+2)*Exp_msk1)
michael@0 2122 /* round even ==> */
michael@0 2123 /* accept rv */
michael@0 2124 break;
michael@0 2125 /* rv = smallest denormal */
michael@0 2126 goto undfl;
michael@0 2127 }
michael@0 2128 }
michael@0 2129 #endif /*Avoid_Underflow*/
michael@0 2130 L = (word0(rv) & Exp_mask) - Exp_msk1;
michael@0 2131 #endif /*Sudden_Underflow}}*/
michael@0 2132 word0(rv) = L | Bndry_mask1;
michael@0 2133 word1(rv) = 0xffffffff;
michael@0 2134 #ifdef IBM
michael@0 2135 goto cont;
michael@0 2136 #else
michael@0 2137 break;
michael@0 2138 #endif
michael@0 2139 }
michael@0 2140 #ifndef ROUND_BIASED
michael@0 2141 if (!(word1(rv) & LSB))
michael@0 2142 break;
michael@0 2143 #endif
michael@0 2144 if (dsign)
michael@0 2145 dval(rv) += ulp(rv);
michael@0 2146 #ifndef ROUND_BIASED
michael@0 2147 else {
michael@0 2148 dval(rv) -= ulp(rv);
michael@0 2149 #ifndef Sudden_Underflow
michael@0 2150 if (!dval(rv))
michael@0 2151 goto undfl;
michael@0 2152 #endif
michael@0 2153 }
michael@0 2154 #ifdef Avoid_Underflow
michael@0 2155 dsign = 1 - dsign;
michael@0 2156 #endif
michael@0 2157 #endif
michael@0 2158 break;
michael@0 2159 }
michael@0 2160 if ((aadj = ratio(delta, bs)) <= 2.) {
michael@0 2161 if (dsign)
michael@0 2162 aadj = dval(aadj1) = 1.;
michael@0 2163 else if (word1(rv) || word0(rv) & Bndry_mask) {
michael@0 2164 #ifndef Sudden_Underflow
michael@0 2165 if (word1(rv) == Tiny1 && !word0(rv))
michael@0 2166 goto undfl;
michael@0 2167 #endif
michael@0 2168 aadj = 1.;
michael@0 2169 dval(aadj1) = -1.;
michael@0 2170 }
michael@0 2171 else {
michael@0 2172 /* special case -- power of FLT_RADIX to be */
michael@0 2173 /* rounded down... */
michael@0 2174
michael@0 2175 if (aadj < 2./FLT_RADIX)
michael@0 2176 aadj = 1./FLT_RADIX;
michael@0 2177 else
michael@0 2178 aadj *= 0.5;
michael@0 2179 dval(aadj1) = -aadj;
michael@0 2180 }
michael@0 2181 }
michael@0 2182 else {
michael@0 2183 aadj *= 0.5;
michael@0 2184 dval(aadj1) = dsign ? aadj : -aadj;
michael@0 2185 #ifdef Check_FLT_ROUNDS
michael@0 2186 switch(Rounding) {
michael@0 2187 case 2: /* towards +infinity */
michael@0 2188 dval(aadj1) -= 0.5;
michael@0 2189 break;
michael@0 2190 case 0: /* towards 0 */
michael@0 2191 case 3: /* towards -infinity */
michael@0 2192 dval(aadj1) += 0.5;
michael@0 2193 }
michael@0 2194 #else
michael@0 2195 if (Flt_Rounds == 0)
michael@0 2196 dval(aadj1) += 0.5;
michael@0 2197 #endif /*Check_FLT_ROUNDS*/
michael@0 2198 }
michael@0 2199 y = word0(rv) & Exp_mask;
michael@0 2200
michael@0 2201 /* Check for overflow */
michael@0 2202
michael@0 2203 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
michael@0 2204 dval(rv0) = dval(rv);
michael@0 2205 word0(rv) -= P*Exp_msk1;
michael@0 2206 adj = dval(aadj1) * ulp(rv);
michael@0 2207 dval(rv) += adj;
michael@0 2208 if ((word0(rv) & Exp_mask) >=
michael@0 2209 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
michael@0 2210 if (word0(rv0) == Big0 && word1(rv0) == Big1)
michael@0 2211 goto ovfl;
michael@0 2212 word0(rv) = Big0;
michael@0 2213 word1(rv) = Big1;
michael@0 2214 goto cont;
michael@0 2215 }
michael@0 2216 else
michael@0 2217 word0(rv) += P*Exp_msk1;
michael@0 2218 }
michael@0 2219 else {
michael@0 2220 #ifdef Avoid_Underflow
michael@0 2221 if (scale && y <= 2*P*Exp_msk1) {
michael@0 2222 if (aadj <= 0x7fffffff) {
michael@0 2223 if ((z = (ULong) aadj) <= 0)
michael@0 2224 z = 1;
michael@0 2225 aadj = z;
michael@0 2226 dval(aadj1) = dsign ? aadj : -aadj;
michael@0 2227 }
michael@0 2228 word0(aadj1) += (2*P+1)*Exp_msk1 - y;
michael@0 2229 }
michael@0 2230 adj = dval(aadj1) * ulp(rv);
michael@0 2231 dval(rv) += adj;
michael@0 2232 #else
michael@0 2233 #ifdef Sudden_Underflow
michael@0 2234 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
michael@0 2235 dval(rv0) = dval(rv);
michael@0 2236 word0(rv) += P*Exp_msk1;
michael@0 2237 adj = dval(aadj1) * ulp(rv);
michael@0 2238 dval(rv) += adj;
michael@0 2239 #ifdef IBM
michael@0 2240 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
michael@0 2241 #else
michael@0 2242 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
michael@0 2243 #endif
michael@0 2244 {
michael@0 2245 if (word0(rv0) == Tiny0
michael@0 2246 && word1(rv0) == Tiny1)
michael@0 2247 goto undfl;
michael@0 2248 word0(rv) = Tiny0;
michael@0 2249 word1(rv) = Tiny1;
michael@0 2250 goto cont;
michael@0 2251 }
michael@0 2252 else
michael@0 2253 word0(rv) -= P*Exp_msk1;
michael@0 2254 }
michael@0 2255 else {
michael@0 2256 adj = dval(aadj1) * ulp(rv);
michael@0 2257 dval(rv) += adj;
michael@0 2258 }
michael@0 2259 #else /*Sudden_Underflow*/
michael@0 2260 /* Compute adj so that the IEEE rounding rules will
michael@0 2261 * correctly round rv + adj in some half-way cases.
michael@0 2262 * If rv * ulp(rv) is denormalized (i.e.,
michael@0 2263 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
michael@0 2264 * trouble from bits lost to denormalization;
michael@0 2265 * example: 1.2e-307 .
michael@0 2266 */
michael@0 2267 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
michael@0 2268 dval(aadj1) = (double)(int)(aadj + 0.5);
michael@0 2269 if (!dsign)
michael@0 2270 dval(aadj1) = -dval(aadj1);
michael@0 2271 }
michael@0 2272 adj = dval(aadj1) * ulp(rv);
michael@0 2273 dval(rv) += adj;
michael@0 2274 #endif /*Sudden_Underflow*/
michael@0 2275 #endif /*Avoid_Underflow*/
michael@0 2276 }
michael@0 2277 z = word0(rv) & Exp_mask;
michael@0 2278 #ifndef SET_INEXACT
michael@0 2279 #ifdef Avoid_Underflow
michael@0 2280 if (!scale)
michael@0 2281 #endif
michael@0 2282 if (y == z) {
michael@0 2283 /* Can we stop now? */
michael@0 2284 L = (Long)aadj;
michael@0 2285 aadj -= L;
michael@0 2286 /* The tolerances below are conservative. */
michael@0 2287 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
michael@0 2288 if (aadj < .4999999 || aadj > .5000001)
michael@0 2289 break;
michael@0 2290 }
michael@0 2291 else if (aadj < .4999999/FLT_RADIX)
michael@0 2292 break;
michael@0 2293 }
michael@0 2294 #endif
michael@0 2295 cont:
michael@0 2296 Bfree(PASS_STATE bb);
michael@0 2297 Bfree(PASS_STATE bd);
michael@0 2298 Bfree(PASS_STATE bs);
michael@0 2299 Bfree(PASS_STATE delta);
michael@0 2300 }
michael@0 2301 #ifdef SET_INEXACT
michael@0 2302 if (inexact) {
michael@0 2303 if (!oldinexact) {
michael@0 2304 word0(rv0) = Exp_1 + (70 << Exp_shift);
michael@0 2305 word1(rv0) = 0;
michael@0 2306 dval(rv0) += 1.;
michael@0 2307 }
michael@0 2308 }
michael@0 2309 else if (!oldinexact)
michael@0 2310 clear_inexact();
michael@0 2311 #endif
michael@0 2312 #ifdef Avoid_Underflow
michael@0 2313 if (scale) {
michael@0 2314 word0(rv0) = Exp_1 - 2*P*Exp_msk1;
michael@0 2315 word1(rv0) = 0;
michael@0 2316 dval(rv) *= dval(rv0);
michael@0 2317 #ifndef NO_ERRNO
michael@0 2318 /* try to avoid the bug of testing an 8087 register value */
michael@0 2319 if (word0(rv) == 0 && word1(rv) == 0)
michael@0 2320 errno = ERANGE;
michael@0 2321 #endif
michael@0 2322 }
michael@0 2323 #endif /* Avoid_Underflow */
michael@0 2324 #ifdef SET_INEXACT
michael@0 2325 if (inexact && !(word0(rv) & Exp_mask)) {
michael@0 2326 /* set underflow bit */
michael@0 2327 dval(rv0) = 1e-300;
michael@0 2328 dval(rv0) *= dval(rv0);
michael@0 2329 }
michael@0 2330 #endif
michael@0 2331 retfree:
michael@0 2332 Bfree(PASS_STATE bb);
michael@0 2333 Bfree(PASS_STATE bd);
michael@0 2334 Bfree(PASS_STATE bs);
michael@0 2335 Bfree(PASS_STATE bd0);
michael@0 2336 Bfree(PASS_STATE delta);
michael@0 2337 ret:
michael@0 2338 if (se)
michael@0 2339 *se = (char *)s;
michael@0 2340 return sign ? -dval(rv) : dval(rv);
michael@0 2341 }
michael@0 2342
michael@0 2343 static int
michael@0 2344 quorem
michael@0 2345 #ifdef KR_headers
michael@0 2346 (b, S) Bigint *b, *S;
michael@0 2347 #else
michael@0 2348 (Bigint *b, Bigint *S)
michael@0 2349 #endif
michael@0 2350 {
michael@0 2351 int n;
michael@0 2352 ULong *bx, *bxe, q, *sx, *sxe;
michael@0 2353 #ifdef ULLong
michael@0 2354 ULLong borrow, carry, y, ys;
michael@0 2355 #else
michael@0 2356 ULong borrow, carry, y, ys;
michael@0 2357 #ifdef Pack_32
michael@0 2358 ULong si, z, zs;
michael@0 2359 #endif
michael@0 2360 #endif
michael@0 2361
michael@0 2362 n = S->wds;
michael@0 2363 #ifdef DEBUG
michael@0 2364 /*debug*/ if (b->wds > n)
michael@0 2365 /*debug*/ Bug("oversize b in quorem");
michael@0 2366 #endif
michael@0 2367 if (b->wds < n)
michael@0 2368 return 0;
michael@0 2369 sx = S->x;
michael@0 2370 sxe = sx + --n;
michael@0 2371 bx = b->x;
michael@0 2372 bxe = bx + n;
michael@0 2373 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
michael@0 2374 #ifdef DEBUG
michael@0 2375 /*debug*/ if (q > 9)
michael@0 2376 /*debug*/ Bug("oversized quotient in quorem");
michael@0 2377 #endif
michael@0 2378 if (q) {
michael@0 2379 borrow = 0;
michael@0 2380 carry = 0;
michael@0 2381 do {
michael@0 2382 #ifdef ULLong
michael@0 2383 ys = *sx++ * (ULLong)q + carry;
michael@0 2384 carry = ys >> 32;
michael@0 2385 y = *bx - (ys & FFFFFFFF) - borrow;
michael@0 2386 borrow = y >> 32 & (ULong)1;
michael@0 2387 *bx++ = (ULong) y & FFFFFFFF;
michael@0 2388 #else
michael@0 2389 #ifdef Pack_32
michael@0 2390 si = *sx++;
michael@0 2391 ys = (si & 0xffff) * q + carry;
michael@0 2392 zs = (si >> 16) * q + (ys >> 16);
michael@0 2393 carry = zs >> 16;
michael@0 2394 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
michael@0 2395 borrow = (y & 0x10000) >> 16;
michael@0 2396 z = (*bx >> 16) - (zs & 0xffff) - borrow;
michael@0 2397 borrow = (z & 0x10000) >> 16;
michael@0 2398 Storeinc(bx, z, y);
michael@0 2399 #else
michael@0 2400 ys = *sx++ * q + carry;
michael@0 2401 carry = ys >> 16;
michael@0 2402 y = *bx - (ys & 0xffff) - borrow;
michael@0 2403 borrow = (y & 0x10000) >> 16;
michael@0 2404 *bx++ = y & 0xffff;
michael@0 2405 #endif
michael@0 2406 #endif
michael@0 2407 }
michael@0 2408 while(sx <= sxe);
michael@0 2409 if (!*bxe) {
michael@0 2410 bx = b->x;
michael@0 2411 while(--bxe > bx && !*bxe)
michael@0 2412 --n;
michael@0 2413 b->wds = n;
michael@0 2414 }
michael@0 2415 }
michael@0 2416 if (cmp(b, S) >= 0) {
michael@0 2417 q++;
michael@0 2418 borrow = 0;
michael@0 2419 carry = 0;
michael@0 2420 bx = b->x;
michael@0 2421 sx = S->x;
michael@0 2422 do {
michael@0 2423 #ifdef ULLong
michael@0 2424 ys = *sx++ + carry;
michael@0 2425 carry = ys >> 32;
michael@0 2426 y = *bx - (ys & FFFFFFFF) - borrow;
michael@0 2427 borrow = y >> 32 & (ULong)1;
michael@0 2428 *bx++ = (ULong) y & FFFFFFFF;
michael@0 2429 #else
michael@0 2430 #ifdef Pack_32
michael@0 2431 si = *sx++;
michael@0 2432 ys = (si & 0xffff) + carry;
michael@0 2433 zs = (si >> 16) + (ys >> 16);
michael@0 2434 carry = zs >> 16;
michael@0 2435 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
michael@0 2436 borrow = (y & 0x10000) >> 16;
michael@0 2437 z = (*bx >> 16) - (zs & 0xffff) - borrow;
michael@0 2438 borrow = (z & 0x10000) >> 16;
michael@0 2439 Storeinc(bx, z, y);
michael@0 2440 #else
michael@0 2441 ys = *sx++ + carry;
michael@0 2442 carry = ys >> 16;
michael@0 2443 y = *bx - (ys & 0xffff) - borrow;
michael@0 2444 borrow = (y & 0x10000) >> 16;
michael@0 2445 *bx++ = y & 0xffff;
michael@0 2446 #endif
michael@0 2447 #endif
michael@0 2448 }
michael@0 2449 while(sx <= sxe);
michael@0 2450 bx = b->x;
michael@0 2451 bxe = bx + n;
michael@0 2452 if (!*bxe) {
michael@0 2453 while(--bxe > bx && !*bxe)
michael@0 2454 --n;
michael@0 2455 b->wds = n;
michael@0 2456 }
michael@0 2457 }
michael@0 2458 return q;
michael@0 2459 }
michael@0 2460
michael@0 2461 #if !defined(MULTIPLE_THREADS) && !defined(NO_GLOBAL_STATE)
michael@0 2462 #define USE_DTOA_RESULT 1
michael@0 2463 static char *dtoa_result;
michael@0 2464 #endif
michael@0 2465
michael@0 2466 static char *
michael@0 2467 #ifdef KR_headers
michael@0 2468 rv_alloc(STATE_PARAM i) STATE_PARAM_DECL int i;
michael@0 2469 #else
michael@0 2470 rv_alloc(STATE_PARAM int i)
michael@0 2471 #endif
michael@0 2472 {
michael@0 2473 int j, k, *r;
michael@0 2474
michael@0 2475 j = sizeof(ULong);
michael@0 2476 for(k = 0;
michael@0 2477 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned) i;
michael@0 2478 j <<= 1)
michael@0 2479 k++;
michael@0 2480 r = (int*)Balloc(PASS_STATE k);
michael@0 2481 *r = k;
michael@0 2482 return
michael@0 2483 #ifdef USE_DTOA_RESULT
michael@0 2484 dtoa_result =
michael@0 2485 #endif
michael@0 2486 (char *)(r+1);
michael@0 2487 }
michael@0 2488
michael@0 2489 static char *
michael@0 2490 #ifdef KR_headers
michael@0 2491 nrv_alloc(STATE_PARAM s, rve, n) STATE_PARAM_DECL char *s, **rve; int n;
michael@0 2492 #else
michael@0 2493 nrv_alloc(STATE_PARAM CONST char *s, char **rve, int n)
michael@0 2494 #endif
michael@0 2495 {
michael@0 2496 char *rv, *t;
michael@0 2497
michael@0 2498 t = rv = rv_alloc(PASS_STATE n);
michael@0 2499 while((*t = *s++)) t++;
michael@0 2500 if (rve)
michael@0 2501 *rve = t;
michael@0 2502 return rv;
michael@0 2503 }
michael@0 2504
michael@0 2505 /* freedtoa(s) must be used to free values s returned by dtoa
michael@0 2506 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
michael@0 2507 * but for consistency with earlier versions of dtoa, it is optional
michael@0 2508 * when MULTIPLE_THREADS is not defined.
michael@0 2509 */
michael@0 2510
michael@0 2511 static void
michael@0 2512 #ifdef KR_headers
michael@0 2513 freedtoa(STATE_PARAM s) STATE_PARAM_DECL char *s;
michael@0 2514 #else
michael@0 2515 freedtoa(STATE_PARAM char *s)
michael@0 2516 #endif
michael@0 2517 {
michael@0 2518 Bigint *b = (Bigint *)((int *)s - 1);
michael@0 2519 b->maxwds = 1 << (b->k = *(int*)b);
michael@0 2520 Bfree(PASS_STATE b);
michael@0 2521 #ifdef USE_DTOA_RESULT
michael@0 2522 if (s == dtoa_result)
michael@0 2523 dtoa_result = 0;
michael@0 2524 #endif
michael@0 2525 }
michael@0 2526
michael@0 2527 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
michael@0 2528 *
michael@0 2529 * Inspired by "How to Print Floating-Point Numbers Accurately" by
michael@0 2530 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
michael@0 2531 *
michael@0 2532 * Modifications:
michael@0 2533 * 1. Rather than iterating, we use a simple numeric overestimate
michael@0 2534 * to determine k = floor(log10(d)). We scale relevant
michael@0 2535 * quantities using O(log2(k)) rather than O(k) multiplications.
michael@0 2536 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
michael@0 2537 * try to generate digits strictly left to right. Instead, we
michael@0 2538 * compute with fewer bits and propagate the carry if necessary
michael@0 2539 * when rounding the final digit up. This is often faster.
michael@0 2540 * 3. Under the assumption that input will be rounded nearest,
michael@0 2541 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
michael@0 2542 * That is, we allow equality in stopping tests when the
michael@0 2543 * round-nearest rule will give the same floating-point value
michael@0 2544 * as would satisfaction of the stopping test with strict
michael@0 2545 * inequality.
michael@0 2546 * 4. We remove common factors of powers of 2 from relevant
michael@0 2547 * quantities.
michael@0 2548 * 5. When converting floating-point integers less than 1e16,
michael@0 2549 * we use floating-point arithmetic rather than resorting
michael@0 2550 * to multiple-precision integers.
michael@0 2551 * 6. When asked to produce fewer than 15 digits, we first try
michael@0 2552 * to get by with floating-point arithmetic; we resort to
michael@0 2553 * multiple-precision integer arithmetic only if we cannot
michael@0 2554 * guarantee that the floating-point calculation has given
michael@0 2555 * the correctly rounded result. For k requested digits and
michael@0 2556 * "uniformly" distributed input, the probability is
michael@0 2557 * something like 10^(k-15) that we must resort to the Long
michael@0 2558 * calculation.
michael@0 2559 */
michael@0 2560
michael@0 2561 static char *
michael@0 2562 dtoa
michael@0 2563 #ifdef KR_headers
michael@0 2564 (STATE_PARAM d, mode, ndigits, decpt, sign, rve)
michael@0 2565 STATE_PARAM_DECL U d; int mode, ndigits, *decpt, *sign; char **rve;
michael@0 2566 #else
michael@0 2567 (STATE_PARAM U d, int mode, int ndigits, int *decpt, int *sign, char **rve)
michael@0 2568 #endif
michael@0 2569 {
michael@0 2570 /* Arguments ndigits, decpt, sign are similar to those
michael@0 2571 of ecvt and fcvt; trailing zeros are suppressed from
michael@0 2572 the returned string. If not null, *rve is set to point
michael@0 2573 to the end of the return value. If d is +-Infinity or NaN,
michael@0 2574 then *decpt is set to 9999.
michael@0 2575
michael@0 2576 mode:
michael@0 2577 0 ==> shortest string that yields d when read in
michael@0 2578 and rounded to nearest.
michael@0 2579 1 ==> like 0, but with Steele & White stopping rule;
michael@0 2580 e.g. with IEEE P754 arithmetic , mode 0 gives
michael@0 2581 1e23 whereas mode 1 gives 9.999999999999999e22.
michael@0 2582 2 ==> max(1,ndigits) significant digits. This gives a
michael@0 2583 return value similar to that of ecvt, except
michael@0 2584 that trailing zeros are suppressed.
michael@0 2585 3 ==> through ndigits past the decimal point. This
michael@0 2586 gives a return value similar to that from fcvt,
michael@0 2587 except that trailing zeros are suppressed, and
michael@0 2588 ndigits can be negative.
michael@0 2589 4,5 ==> similar to 2 and 3, respectively, but (in
michael@0 2590 round-nearest mode) with the tests of mode 0 to
michael@0 2591 possibly return a shorter string that rounds to d.
michael@0 2592 With IEEE arithmetic and compilation with
michael@0 2593 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
michael@0 2594 as modes 2 and 3 when FLT_ROUNDS != 1.
michael@0 2595 6-9 ==> Debugging modes similar to mode - 4: don't try
michael@0 2596 fast floating-point estimate (if applicable).
michael@0 2597
michael@0 2598 Values of mode other than 0-9 are treated as mode 0.
michael@0 2599
michael@0 2600 Sufficient space is allocated to the return value
michael@0 2601 to hold the suppressed trailing zeros.
michael@0 2602 */
michael@0 2603
michael@0 2604 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
michael@0 2605 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
michael@0 2606 spec_case, try_quick;
michael@0 2607 Long L;
michael@0 2608 #ifndef Sudden_Underflow
michael@0 2609 int denorm;
michael@0 2610 ULong x;
michael@0 2611 #endif
michael@0 2612 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
michael@0 2613 U d2, eps;
michael@0 2614 double ds;
michael@0 2615 char *s, *s0;
michael@0 2616 #ifdef Honor_FLT_ROUNDS
michael@0 2617 int rounding;
michael@0 2618 #endif
michael@0 2619 #ifdef SET_INEXACT
michael@0 2620 int inexact, oldinexact;
michael@0 2621 #endif
michael@0 2622
michael@0 2623 #ifdef __GNUC__
michael@0 2624 ilim = ilim1 = 0;
michael@0 2625 mlo = NULL;
michael@0 2626 #endif
michael@0 2627
michael@0 2628 #ifdef USE_DTOA_RESULT
michael@0 2629 if (dtoa_result) {
michael@0 2630 freedtoa(PASS_STATE dtoa_result);
michael@0 2631 dtoa_result = 0;
michael@0 2632 }
michael@0 2633 #endif
michael@0 2634
michael@0 2635 if (word0(d) & Sign_bit) {
michael@0 2636 /* set sign for everything, including 0's and NaNs */
michael@0 2637 *sign = 1;
michael@0 2638 word0(d) &= ~Sign_bit; /* clear sign bit */
michael@0 2639 }
michael@0 2640 else
michael@0 2641 *sign = 0;
michael@0 2642
michael@0 2643 #if defined(IEEE_Arith) + defined(VAX)
michael@0 2644 #ifdef IEEE_Arith
michael@0 2645 if ((word0(d) & Exp_mask) == Exp_mask)
michael@0 2646 #else
michael@0 2647 if (word0(d) == 0x8000)
michael@0 2648 #endif
michael@0 2649 {
michael@0 2650 /* Infinity or NaN */
michael@0 2651 *decpt = 9999;
michael@0 2652 #ifdef IEEE_Arith
michael@0 2653 if (!word1(d) && !(word0(d) & 0xfffff))
michael@0 2654 return nrv_alloc(PASS_STATE "Infinity", rve, 8);
michael@0 2655 #endif
michael@0 2656 return nrv_alloc(PASS_STATE "NaN", rve, 3);
michael@0 2657 }
michael@0 2658 #endif
michael@0 2659 #ifdef IBM
michael@0 2660 dval(d) += 0; /* normalize */
michael@0 2661 #endif
michael@0 2662 if (!dval(d)) {
michael@0 2663 *decpt = 1;
michael@0 2664 return nrv_alloc(PASS_STATE "0", rve, 1);
michael@0 2665 }
michael@0 2666
michael@0 2667 #ifdef SET_INEXACT
michael@0 2668 try_quick = oldinexact = get_inexact();
michael@0 2669 inexact = 1;
michael@0 2670 #endif
michael@0 2671 #ifdef Honor_FLT_ROUNDS
michael@0 2672 if ((rounding = Flt_Rounds) >= 2) {
michael@0 2673 if (*sign)
michael@0 2674 rounding = rounding == 2 ? 0 : 2;
michael@0 2675 else
michael@0 2676 if (rounding != 2)
michael@0 2677 rounding = 0;
michael@0 2678 }
michael@0 2679 #endif
michael@0 2680
michael@0 2681 b = d2b(PASS_STATE d, &be, &bbits);
michael@0 2682 #ifdef Sudden_Underflow
michael@0 2683 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
michael@0 2684 #else
michael@0 2685 if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
michael@0 2686 #endif
michael@0 2687 dval(d2) = dval(d);
michael@0 2688 word0(d2) &= Frac_mask1;
michael@0 2689 word0(d2) |= Exp_11;
michael@0 2690 #ifdef IBM
michael@0 2691 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
michael@0 2692 dval(d2) /= 1 << j;
michael@0 2693 #endif
michael@0 2694
michael@0 2695 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
michael@0 2696 * log10(x) = log(x) / log(10)
michael@0 2697 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
michael@0 2698 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
michael@0 2699 *
michael@0 2700 * This suggests computing an approximation k to log10(d) by
michael@0 2701 *
michael@0 2702 * k = (i - Bias)*0.301029995663981
michael@0 2703 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
michael@0 2704 *
michael@0 2705 * We want k to be too large rather than too small.
michael@0 2706 * The error in the first-order Taylor series approximation
michael@0 2707 * is in our favor, so we just round up the constant enough
michael@0 2708 * to compensate for any error in the multiplication of
michael@0 2709 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
michael@0 2710 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
michael@0 2711 * adding 1e-13 to the constant term more than suffices.
michael@0 2712 * Hence we adjust the constant term to 0.1760912590558.
michael@0 2713 * (We could get a more accurate k by invoking log10,
michael@0 2714 * but this is probably not worthwhile.)
michael@0 2715 */
michael@0 2716
michael@0 2717 i -= Bias;
michael@0 2718 #ifdef IBM
michael@0 2719 i <<= 2;
michael@0 2720 i += j;
michael@0 2721 #endif
michael@0 2722 #ifndef Sudden_Underflow
michael@0 2723 denorm = 0;
michael@0 2724 }
michael@0 2725 else {
michael@0 2726 /* d is denormalized */
michael@0 2727
michael@0 2728 i = bbits + be + (Bias + (P-1) - 1);
michael@0 2729 x = i > 32 ? word0(d) << (64 - i) | word1(d) >> (i - 32)
michael@0 2730 : word1(d) << (32 - i);
michael@0 2731 dval(d2) = x;
michael@0 2732 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
michael@0 2733 i -= (Bias + (P-1) - 1) + 1;
michael@0 2734 denorm = 1;
michael@0 2735 }
michael@0 2736 #endif
michael@0 2737 ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
michael@0 2738 k = (int)ds;
michael@0 2739 if (ds < 0. && ds != k)
michael@0 2740 k--; /* want k = floor(ds) */
michael@0 2741 k_check = 1;
michael@0 2742 if (k >= 0 && k <= Ten_pmax) {
michael@0 2743 if (dval(d) < tens[k])
michael@0 2744 k--;
michael@0 2745 k_check = 0;
michael@0 2746 }
michael@0 2747 j = bbits - i - 1;
michael@0 2748 if (j >= 0) {
michael@0 2749 b2 = 0;
michael@0 2750 s2 = j;
michael@0 2751 }
michael@0 2752 else {
michael@0 2753 b2 = -j;
michael@0 2754 s2 = 0;
michael@0 2755 }
michael@0 2756 if (k >= 0) {
michael@0 2757 b5 = 0;
michael@0 2758 s5 = k;
michael@0 2759 s2 += k;
michael@0 2760 }
michael@0 2761 else {
michael@0 2762 b2 -= k;
michael@0 2763 b5 = -k;
michael@0 2764 s5 = 0;
michael@0 2765 }
michael@0 2766 if (mode < 0 || mode > 9)
michael@0 2767 mode = 0;
michael@0 2768
michael@0 2769 #ifndef SET_INEXACT
michael@0 2770 #ifdef Check_FLT_ROUNDS
michael@0 2771 try_quick = Rounding == 1;
michael@0 2772 #else
michael@0 2773 try_quick = 1;
michael@0 2774 #endif
michael@0 2775 #endif /*SET_INEXACT*/
michael@0 2776
michael@0 2777 if (mode > 5) {
michael@0 2778 mode -= 4;
michael@0 2779 try_quick = 0;
michael@0 2780 }
michael@0 2781 leftright = 1;
michael@0 2782 switch(mode) {
michael@0 2783 case 0:
michael@0 2784 case 1:
michael@0 2785 ilim = ilim1 = -1;
michael@0 2786 i = 18;
michael@0 2787 ndigits = 0;
michael@0 2788 break;
michael@0 2789 case 2:
michael@0 2790 leftright = 0;
michael@0 2791 /* no break */
michael@0 2792 case 4:
michael@0 2793 if (ndigits <= 0)
michael@0 2794 ndigits = 1;
michael@0 2795 ilim = ilim1 = i = ndigits;
michael@0 2796 break;
michael@0 2797 case 3:
michael@0 2798 leftright = 0;
michael@0 2799 /* no break */
michael@0 2800 case 5:
michael@0 2801 i = ndigits + k + 1;
michael@0 2802 ilim = i;
michael@0 2803 ilim1 = i - 1;
michael@0 2804 if (i <= 0)
michael@0 2805 i = 1;
michael@0 2806 }
michael@0 2807 s = s0 = rv_alloc(PASS_STATE i);
michael@0 2808
michael@0 2809 #ifdef Honor_FLT_ROUNDS
michael@0 2810 if (mode > 1 && rounding != 1)
michael@0 2811 leftright = 0;
michael@0 2812 #endif
michael@0 2813
michael@0 2814 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
michael@0 2815
michael@0 2816 /* Try to get by with floating-point arithmetic. */
michael@0 2817
michael@0 2818 i = 0;
michael@0 2819 dval(d2) = dval(d);
michael@0 2820 k0 = k;
michael@0 2821 ilim0 = ilim;
michael@0 2822 ieps = 2; /* conservative */
michael@0 2823 if (k > 0) {
michael@0 2824 ds = tens[k&0xf];
michael@0 2825 j = k >> 4;
michael@0 2826 if (j & Bletch) {
michael@0 2827 /* prevent overflows */
michael@0 2828 j &= Bletch - 1;
michael@0 2829 dval(d) /= bigtens[n_bigtens-1];
michael@0 2830 ieps++;
michael@0 2831 }
michael@0 2832 for(; j; j >>= 1, i++)
michael@0 2833 if (j & 1) {
michael@0 2834 ieps++;
michael@0 2835 ds *= bigtens[i];
michael@0 2836 }
michael@0 2837 dval(d) /= ds;
michael@0 2838 }
michael@0 2839 else if ((j1 = -k)) {
michael@0 2840 dval(d) *= tens[j1 & 0xf];
michael@0 2841 for(j = j1 >> 4; j; j >>= 1, i++)
michael@0 2842 if (j & 1) {
michael@0 2843 ieps++;
michael@0 2844 dval(d) *= bigtens[i];
michael@0 2845 }
michael@0 2846 }
michael@0 2847 if (k_check && dval(d) < 1. && ilim > 0) {
michael@0 2848 if (ilim1 <= 0)
michael@0 2849 goto fast_failed;
michael@0 2850 ilim = ilim1;
michael@0 2851 k--;
michael@0 2852 dval(d) *= 10.;
michael@0 2853 ieps++;
michael@0 2854 }
michael@0 2855 dval(eps) = ieps*dval(d) + 7.;
michael@0 2856 word0(eps) -= (P-1)*Exp_msk1;
michael@0 2857 if (ilim == 0) {
michael@0 2858 S = mhi = 0;
michael@0 2859 dval(d) -= 5.;
michael@0 2860 if (dval(d) > dval(eps))
michael@0 2861 goto one_digit;
michael@0 2862 if (dval(d) < -dval(eps))
michael@0 2863 goto no_digits;
michael@0 2864 goto fast_failed;
michael@0 2865 }
michael@0 2866 #ifndef No_leftright
michael@0 2867 if (leftright) {
michael@0 2868 /* Use Steele & White method of only
michael@0 2869 * generating digits needed.
michael@0 2870 */
michael@0 2871 dval(eps) = 0.5/tens[ilim-1] - dval(eps);
michael@0 2872 for(i = 0;;) {
michael@0 2873 L = (ULong) dval(d);
michael@0 2874 dval(d) -= L;
michael@0 2875 *s++ = '0' + (int)L;
michael@0 2876 if (dval(d) < dval(eps))
michael@0 2877 goto ret1;
michael@0 2878 if (1. - dval(d) < dval(eps))
michael@0 2879 goto bump_up;
michael@0 2880 if (++i >= ilim)
michael@0 2881 break;
michael@0 2882 dval(eps) *= 10.;
michael@0 2883 dval(d) *= 10.;
michael@0 2884 }
michael@0 2885 }
michael@0 2886 else {
michael@0 2887 #endif
michael@0 2888 /* Generate ilim digits, then fix them up. */
michael@0 2889 dval(eps) *= tens[ilim-1];
michael@0 2890 for(i = 1;; i++, dval(d) *= 10.) {
michael@0 2891 L = (Long)(dval(d));
michael@0 2892 if (!(dval(d) -= L))
michael@0 2893 ilim = i;
michael@0 2894 *s++ = '0' + (int)L;
michael@0 2895 if (i == ilim) {
michael@0 2896 if (dval(d) > 0.5 + dval(eps))
michael@0 2897 goto bump_up;
michael@0 2898 else if (dval(d) < 0.5 - dval(eps)) {
michael@0 2899 while(*--s == '0');
michael@0 2900 s++;
michael@0 2901 goto ret1;
michael@0 2902 }
michael@0 2903 break;
michael@0 2904 }
michael@0 2905 }
michael@0 2906 #ifndef No_leftright
michael@0 2907 }
michael@0 2908 #endif
michael@0 2909 fast_failed:
michael@0 2910 s = s0;
michael@0 2911 dval(d) = dval(d2);
michael@0 2912 k = k0;
michael@0 2913 ilim = ilim0;
michael@0 2914 }
michael@0 2915
michael@0 2916 /* Do we have a "small" integer? */
michael@0 2917
michael@0 2918 if (be >= 0 && k <= Int_max) {
michael@0 2919 /* Yes. */
michael@0 2920 ds = tens[k];
michael@0 2921 if (ndigits < 0 && ilim <= 0) {
michael@0 2922 S = mhi = 0;
michael@0 2923 if (ilim < 0 || dval(d) < 5*ds)
michael@0 2924 goto no_digits;
michael@0 2925 goto one_digit;
michael@0 2926 }
michael@0 2927 for(i = 1;; i++, dval(d) *= 10.) {
michael@0 2928 L = (Long)(dval(d) / ds);
michael@0 2929 dval(d) -= L*ds;
michael@0 2930 #ifdef Check_FLT_ROUNDS
michael@0 2931 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
michael@0 2932 if (dval(d) < 0) {
michael@0 2933 L--;
michael@0 2934 dval(d) += ds;
michael@0 2935 }
michael@0 2936 #endif
michael@0 2937 *s++ = '0' + (int)L;
michael@0 2938 if (!dval(d)) {
michael@0 2939 #ifdef SET_INEXACT
michael@0 2940 inexact = 0;
michael@0 2941 #endif
michael@0 2942 break;
michael@0 2943 }
michael@0 2944 if (i == ilim) {
michael@0 2945 #ifdef Honor_FLT_ROUNDS
michael@0 2946 if (mode > 1)
michael@0 2947 switch(rounding) {
michael@0 2948 case 0: goto ret1;
michael@0 2949 case 2: goto bump_up;
michael@0 2950 }
michael@0 2951 #endif
michael@0 2952 dval(d) += dval(d);
michael@0 2953 if (dval(d) > ds || (dval(d) == ds && L & 1)) {
michael@0 2954 bump_up:
michael@0 2955 while(*--s == '9')
michael@0 2956 if (s == s0) {
michael@0 2957 k++;
michael@0 2958 *s = '0';
michael@0 2959 break;
michael@0 2960 }
michael@0 2961 ++*s++;
michael@0 2962 }
michael@0 2963 break;
michael@0 2964 }
michael@0 2965 }
michael@0 2966 goto ret1;
michael@0 2967 }
michael@0 2968
michael@0 2969 m2 = b2;
michael@0 2970 m5 = b5;
michael@0 2971 mhi = mlo = 0;
michael@0 2972 if (leftright) {
michael@0 2973 i =
michael@0 2974 #ifndef Sudden_Underflow
michael@0 2975 denorm ? be + (Bias + (P-1) - 1 + 1) :
michael@0 2976 #endif
michael@0 2977 #ifdef IBM
michael@0 2978 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
michael@0 2979 #else
michael@0 2980 1 + P - bbits;
michael@0 2981 #endif
michael@0 2982 b2 += i;
michael@0 2983 s2 += i;
michael@0 2984 mhi = i2b(PASS_STATE 1);
michael@0 2985 }
michael@0 2986 if (m2 > 0 && s2 > 0) {
michael@0 2987 i = m2 < s2 ? m2 : s2;
michael@0 2988 b2 -= i;
michael@0 2989 m2 -= i;
michael@0 2990 s2 -= i;
michael@0 2991 }
michael@0 2992 if (b5 > 0) {
michael@0 2993 if (leftright) {
michael@0 2994 if (m5 > 0) {
michael@0 2995 mhi = pow5mult(PASS_STATE mhi, m5);
michael@0 2996 b1 = mult(PASS_STATE mhi, b);
michael@0 2997 Bfree(PASS_STATE b);
michael@0 2998 b = b1;
michael@0 2999 }
michael@0 3000 if ((j = b5 - m5))
michael@0 3001 b = pow5mult(PASS_STATE b, j);
michael@0 3002 }
michael@0 3003 else
michael@0 3004 b = pow5mult(PASS_STATE b, b5);
michael@0 3005 }
michael@0 3006 S = i2b(PASS_STATE 1);
michael@0 3007 if (s5 > 0)
michael@0 3008 S = pow5mult(PASS_STATE S, s5);
michael@0 3009
michael@0 3010 /* Check for special case that d is a normalized power of 2. */
michael@0 3011
michael@0 3012 spec_case = 0;
michael@0 3013 if ((mode < 2 || leftright)
michael@0 3014 #ifdef Honor_FLT_ROUNDS
michael@0 3015 && rounding == 1
michael@0 3016 #endif
michael@0 3017 ) {
michael@0 3018 if (!word1(d) && !(word0(d) & Bndry_mask)
michael@0 3019 #ifndef Sudden_Underflow
michael@0 3020 && word0(d) & (Exp_mask & ~Exp_msk1)
michael@0 3021 #endif
michael@0 3022 ) {
michael@0 3023 /* The special case */
michael@0 3024 b2 += Log2P;
michael@0 3025 s2 += Log2P;
michael@0 3026 spec_case = 1;
michael@0 3027 }
michael@0 3028 }
michael@0 3029
michael@0 3030 /* Arrange for convenient computation of quotients:
michael@0 3031 * shift left if necessary so divisor has 4 leading 0 bits.
michael@0 3032 *
michael@0 3033 * Perhaps we should just compute leading 28 bits of S once
michael@0 3034 * and for all and pass them and a shift to quorem, so it
michael@0 3035 * can do shifts and ors to compute the numerator for q.
michael@0 3036 */
michael@0 3037 #ifdef Pack_32
michael@0 3038 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
michael@0 3039 i = 32 - i;
michael@0 3040 #else
michael@0 3041 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
michael@0 3042 i = 16 - i;
michael@0 3043 #endif
michael@0 3044 if (i > 4) {
michael@0 3045 i -= 4;
michael@0 3046 b2 += i;
michael@0 3047 m2 += i;
michael@0 3048 s2 += i;
michael@0 3049 }
michael@0 3050 else if (i < 4) {
michael@0 3051 i += 28;
michael@0 3052 b2 += i;
michael@0 3053 m2 += i;
michael@0 3054 s2 += i;
michael@0 3055 }
michael@0 3056 if (b2 > 0)
michael@0 3057 b = lshift(PASS_STATE b, b2);
michael@0 3058 if (s2 > 0)
michael@0 3059 S = lshift(PASS_STATE S, s2);
michael@0 3060 if (k_check) {
michael@0 3061 if (cmp(b,S) < 0) {
michael@0 3062 k--;
michael@0 3063 b = multadd(PASS_STATE b, 10, 0); /* we botched the k estimate */
michael@0 3064 if (leftright)
michael@0 3065 mhi = multadd(PASS_STATE mhi, 10, 0);
michael@0 3066 ilim = ilim1;
michael@0 3067 }
michael@0 3068 }
michael@0 3069 if (ilim <= 0 && (mode == 3 || mode == 5)) {
michael@0 3070 if (ilim < 0 || cmp(b,S = multadd(PASS_STATE S,5,0)) < 0) {
michael@0 3071 /* no digits, fcvt style */
michael@0 3072 no_digits:
michael@0 3073 /* MOZILLA CHANGE: Always return a non-empty string. */
michael@0 3074 *s++ = '0';
michael@0 3075 k = 0;
michael@0 3076 goto ret;
michael@0 3077 }
michael@0 3078 one_digit:
michael@0 3079 *s++ = '1';
michael@0 3080 k++;
michael@0 3081 goto ret;
michael@0 3082 }
michael@0 3083 if (leftright) {
michael@0 3084 if (m2 > 0)
michael@0 3085 mhi = lshift(PASS_STATE mhi, m2);
michael@0 3086
michael@0 3087 /* Compute mlo -- check for special case
michael@0 3088 * that d is a normalized power of 2.
michael@0 3089 */
michael@0 3090
michael@0 3091 mlo = mhi;
michael@0 3092 if (spec_case) {
michael@0 3093 mhi = Balloc(PASS_STATE mhi->k);
michael@0 3094 Bcopy(mhi, mlo);
michael@0 3095 mhi = lshift(PASS_STATE mhi, Log2P);
michael@0 3096 }
michael@0 3097
michael@0 3098 for(i = 1;;i++) {
michael@0 3099 dig = quorem(b,S) + '0';
michael@0 3100 /* Do we yet have the shortest decimal string
michael@0 3101 * that will round to d?
michael@0 3102 */
michael@0 3103 j = cmp(b, mlo);
michael@0 3104 delta = diff(PASS_STATE S, mhi);
michael@0 3105 j1 = delta->sign ? 1 : cmp(b, delta);
michael@0 3106 Bfree(PASS_STATE delta);
michael@0 3107 #ifndef ROUND_BIASED
michael@0 3108 if (j1 == 0 && mode != 1 && !(word1(d) & 1)
michael@0 3109 #ifdef Honor_FLT_ROUNDS
michael@0 3110 && rounding >= 1
michael@0 3111 #endif
michael@0 3112 ) {
michael@0 3113 if (dig == '9')
michael@0 3114 goto round_9_up;
michael@0 3115 if (j > 0)
michael@0 3116 dig++;
michael@0 3117 #ifdef SET_INEXACT
michael@0 3118 else if (!b->x[0] && b->wds <= 1)
michael@0 3119 inexact = 0;
michael@0 3120 #endif
michael@0 3121 *s++ = dig;
michael@0 3122 goto ret;
michael@0 3123 }
michael@0 3124 #endif
michael@0 3125 if (j < 0 || (j == 0 && mode != 1
michael@0 3126 #ifndef ROUND_BIASED
michael@0 3127 && !(word1(d) & 1)
michael@0 3128 #endif
michael@0 3129 )) {
michael@0 3130 if (!b->x[0] && b->wds <= 1) {
michael@0 3131 #ifdef SET_INEXACT
michael@0 3132 inexact = 0;
michael@0 3133 #endif
michael@0 3134 goto accept_dig;
michael@0 3135 }
michael@0 3136 #ifdef Honor_FLT_ROUNDS
michael@0 3137 if (mode > 1)
michael@0 3138 switch(rounding) {
michael@0 3139 case 0: goto accept_dig;
michael@0 3140 case 2: goto keep_dig;
michael@0 3141 }
michael@0 3142 #endif /*Honor_FLT_ROUNDS*/
michael@0 3143 if (j1 > 0) {
michael@0 3144 b = lshift(PASS_STATE b, 1);
michael@0 3145 j1 = cmp(b, S);
michael@0 3146 if ((j1 > 0 || (j1 == 0 && dig & 1))
michael@0 3147 && dig++ == '9')
michael@0 3148 goto round_9_up;
michael@0 3149 }
michael@0 3150 accept_dig:
michael@0 3151 *s++ = dig;
michael@0 3152 goto ret;
michael@0 3153 }
michael@0 3154 if (j1 > 0) {
michael@0 3155 #ifdef Honor_FLT_ROUNDS
michael@0 3156 if (!rounding)
michael@0 3157 goto accept_dig;
michael@0 3158 #endif
michael@0 3159 if (dig == '9') { /* possible if i == 1 */
michael@0 3160 round_9_up:
michael@0 3161 *s++ = '9';
michael@0 3162 goto roundoff;
michael@0 3163 }
michael@0 3164 *s++ = dig + 1;
michael@0 3165 goto ret;
michael@0 3166 }
michael@0 3167 #ifdef Honor_FLT_ROUNDS
michael@0 3168 keep_dig:
michael@0 3169 #endif
michael@0 3170 *s++ = dig;
michael@0 3171 if (i == ilim)
michael@0 3172 break;
michael@0 3173 b = multadd(PASS_STATE b, 10, 0);
michael@0 3174 if (mlo == mhi)
michael@0 3175 mlo = mhi = multadd(PASS_STATE mhi, 10, 0);
michael@0 3176 else {
michael@0 3177 mlo = multadd(PASS_STATE mlo, 10, 0);
michael@0 3178 mhi = multadd(PASS_STATE mhi, 10, 0);
michael@0 3179 }
michael@0 3180 }
michael@0 3181 }
michael@0 3182 else
michael@0 3183 for(i = 1;; i++) {
michael@0 3184 *s++ = dig = quorem(b,S) + '0';
michael@0 3185 if (!b->x[0] && b->wds <= 1) {
michael@0 3186 #ifdef SET_INEXACT
michael@0 3187 inexact = 0;
michael@0 3188 #endif
michael@0 3189 goto ret;
michael@0 3190 }
michael@0 3191 if (i >= ilim)
michael@0 3192 break;
michael@0 3193 b = multadd(PASS_STATE b, 10, 0);
michael@0 3194 }
michael@0 3195
michael@0 3196 /* Round off last digit */
michael@0 3197
michael@0 3198 #ifdef Honor_FLT_ROUNDS
michael@0 3199 switch(rounding) {
michael@0 3200 case 0: goto trimzeros;
michael@0 3201 case 2: goto roundoff;
michael@0 3202 }
michael@0 3203 #endif
michael@0 3204 b = lshift(PASS_STATE b, 1);
michael@0 3205 j = cmp(b, S);
michael@0 3206 if (j >= 0) { /* ECMA compatible rounding needed by Spidermonkey */
michael@0 3207 roundoff:
michael@0 3208 while(*--s == '9')
michael@0 3209 if (s == s0) {
michael@0 3210 k++;
michael@0 3211 *s++ = '1';
michael@0 3212 goto ret;
michael@0 3213 }
michael@0 3214 ++*s++;
michael@0 3215 }
michael@0 3216 else {
michael@0 3217 #ifdef Honor_FLT_ROUNDS
michael@0 3218 trimzeros:
michael@0 3219 #endif
michael@0 3220 while(*--s == '0');
michael@0 3221 s++;
michael@0 3222 }
michael@0 3223 ret:
michael@0 3224 Bfree(PASS_STATE S);
michael@0 3225 if (mhi) {
michael@0 3226 if (mlo && mlo != mhi)
michael@0 3227 Bfree(PASS_STATE mlo);
michael@0 3228 Bfree(PASS_STATE mhi);
michael@0 3229 }
michael@0 3230 ret1:
michael@0 3231 #ifdef SET_INEXACT
michael@0 3232 if (inexact) {
michael@0 3233 if (!oldinexact) {
michael@0 3234 word0(d) = Exp_1 + (70 << Exp_shift);
michael@0 3235 word1(d) = 0;
michael@0 3236 dval(d) += 1.;
michael@0 3237 }
michael@0 3238 }
michael@0 3239 else if (!oldinexact)
michael@0 3240 clear_inexact();
michael@0 3241 #endif
michael@0 3242 Bfree(PASS_STATE b);
michael@0 3243 *s = 0;
michael@0 3244 *decpt = k + 1;
michael@0 3245 if (rve)
michael@0 3246 *rve = s;
michael@0 3247 return s0;
michael@0 3248 }

mercurial