diff options
Diffstat (limited to 'media/sphinxbase/src/libsphinxbase/util/dtoa.c')
| -rw-r--r-- | media/sphinxbase/src/libsphinxbase/util/dtoa.c | 2979 |
1 files changed, 2979 insertions, 0 deletions
diff --git a/media/sphinxbase/src/libsphinxbase/util/dtoa.c b/media/sphinxbase/src/libsphinxbase/util/dtoa.c new file mode 100644 index 000000000..4673ae003 --- /dev/null +++ b/media/sphinxbase/src/libsphinxbase/util/dtoa.c @@ -0,0 +1,2979 @@ +/**************************************************************** + * + * The author of this software is David M. Gay. + * + * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. + * + * Permission to use, copy, modify, and distribute this software for any + * purpose without fee is hereby granted, provided that this entire notice + * is included in all copies of any software which is or includes a copy + * or modification of this software and in all copies of the supporting + * documentation for such software. + * + * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED + * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY + * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY + * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. + * + ***************************************************************/ + +/**************************************************************** + * This is dtoa.c by David M. Gay, downloaded from + * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for + * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith. + * It was taken from Python distribution then and imported into sphinxbase. + * Python version is preferred due to cleanups, though original + * version at netlib is still maintained. + * + * Please remember to check http://www.netlib.org/fp regularly for bugfixes and updates. + * + * The major modifications from Gay's original code are as follows: + * + * 0. The original code has been specialized to Sphinxbase's needs by removing + * many of the #ifdef'd sections. In particular, code to support VAX and + * IBM floating-point formats, hex NaNs, hex floats, locale-aware + * treatment of the decimal point, and setting of the inexact flag have + * been removed. + * + * 1. We use cdk_calloc and ckd_free in place of malloc and free. + * + * 2. The public functions strtod, dtoa and freedtoa all now have + * a sb_ prefix. + * + * 3. Instead of assuming that malloc always succeeds, we thread + * malloc failures through the code. The functions + * + * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b + * + * of return type *Bigint all return NULL to indicate a malloc failure. + * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on + * failure. bigcomp now has return type int (it used to be void) and + * returns -1 on failure and 0 otherwise. sb_dtoa returns NULL + * on failure. sb_strtod indicates failure due to malloc failure + * by returning -1.0, setting errno=ENOMEM and *se to s00. + * + * 4. The static variable dtoa_result has been removed. Callers of + * sb_dtoa are expected to call sb_freedtoa to free the memory allocated + * by sb_dtoa. + * + * 5. The code has been reformatted to better fit with C style. + * + * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory + * that hasn't been MALLOC'ed, private_mem should only be used when k <= + * Kmax. + * + * 7. sb_strtod has been modified so that it doesn't accept strings with + * leading whitespace. + * + ***************************************************************/ + +/* Please send bug reports for the original dtoa.c code to David M. Gay (dmg + * at acm dot org, with " at " changed at "@" and " dot " changed to "."). + */ + +/* On a machine with IEEE extended-precision registers, it is + * necessary to specify double-precision (53-bit) rounding precision + * before invoking strtod or dtoa. If the machine uses (the equivalent + * of) Intel 80x87 arithmetic, the call + * _control87(PC_53, MCW_PC); + * does this with many compilers. Whether this or another call is + * appropriate depends on the compiler; for this to work, it may be + * necessary to #include "float.h" or another system-dependent header + * file. + */ + +/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. + * + * This strtod returns a nearest machine number to the input decimal + * string (or sets errno to ERANGE). With IEEE arithmetic, ties are + * broken by the IEEE round-even rule. Otherwise ties are broken by + * biased rounding (add half and chop). + * + * Inspired loosely by William D. Clinger's paper "How to Read Floating + * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. + * + * Modifications: + * + * 1. We only require IEEE, IBM, or VAX double-precision + * arithmetic (not IEEE double-extended). + * 2. We get by with floating-point arithmetic in a case that + * Clinger missed -- when we're computing d * 10^n + * for a small integer d and the integer n is not too + * much larger than 22 (the maximum integer k for which + * we can represent 10^k exactly), we may be able to + * compute (d*10^k) * 10^(e-k) with just one roundoff. + * 3. Rather than a bit-at-a-time adjustment of the binary + * result in the hard case, we use floating-point + * arithmetic to determine the adjustment to within + * one bit; only in really hard cases do we need to + * compute a second residual. + * 4. Because of 3., we don't need a large table of powers of 10 + * for ten-to-e (just some small tables, e.g. of 10^k + * for 0 <= k <= 22). + */ + +/* Linking of sphinxbase's #defines to Gay's #defines starts here. */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include <errno.h> +#include <string.h> +#include <assert.h> +#include <stdio.h> + +#include <sphinxbase/ckd_alloc.h> +#include <sphinxbase/prim_type.h> + +#ifdef WORDS_BIGENDIAN +#define IEEE_MC68k +#else +#define IEEE_8087 +#endif + +#define Long int32 /* ZOMG */ +#define ULong uint32 /* WTF */ +#ifdef HAVE_LONG_LONG +#define ULLong uint64 +#endif + +#define MALLOC ckd_malloc +#define FREE ckd_free + +#define DBL_DIG 15 +#define DBL_MAX_10_EXP 308 +#define DBL_MAX_EXP 1024 +#define FLT_RADIX 2 + +/* maximum permitted exponent value for strtod; exponents larger than + MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP + should fit into an int. */ +#ifndef MAX_ABS_EXP +#define MAX_ABS_EXP 1100000000U +#endif +/* Bound on length of pieces of input strings in sb_strtod; specifically, + this is used to bound the total number of digits ignoring leading zeros and + the number of digits that follow the decimal point. Ideally, MAX_DIGITS + should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the + exponent clipping in sb_strtod can't affect the value of the output. */ +#ifndef MAX_DIGITS +#define MAX_DIGITS 1000000000U +#endif + +/* End sphinxbase #define linking */ + +#ifdef DEBUG +#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} +#endif + +#ifndef PRIVATE_MEM +#define PRIVATE_MEM 2304 +#endif +#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) +static double private_mem[PRIVATE_mem], *pmem_next = private_mem; + +#ifdef __cplusplus +extern "C" { +#endif + +typedef union { double d; ULong L[2]; } U; + +#ifdef IEEE_8087 +#define word0(x) (x)->L[1] +#define word1(x) (x)->L[0] +#else +#define word0(x) (x)->L[0] +#define word1(x) (x)->L[1] +#endif +#define dval(x) (x)->d + +#ifndef STRTOD_DIGLIM +#define STRTOD_DIGLIM 40 +#endif + +/* maximum permitted exponent value for strtod; exponents larger than + MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP + should fit into an int. */ +#ifndef MAX_ABS_EXP +#define MAX_ABS_EXP 1100000000U +#endif +/* Bound on length of pieces of input strings in sb_strtod; specifically, + this is used to bound the total number of digits ignoring leading zeros and + the number of digits that follow the decimal point. Ideally, MAX_DIGITS + should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the + exponent clipping in sb_strtod can't affect the value of the output. */ +#ifndef MAX_DIGITS +#define MAX_DIGITS 1000000000U +#endif + +/* Guard against trying to use the above values on unusual platforms with ints + * of width less than 32 bits. */ +#if MAX_ABS_EXP > 0x7fffffff +#error "MAX_ABS_EXP should fit in an int" +#endif +#if MAX_DIGITS > 0x7fffffff +#error "MAX_DIGITS should fit in an int" +#endif + +/* The following definition of Storeinc is appropriate for MIPS processors. + * An alternative that might be better on some machines is + * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) + */ +#if defined(IEEE_8087) +#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ + ((unsigned short *)a)[0] = (unsigned short)c, a++) +#else +#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ + ((unsigned short *)a)[1] = (unsigned short)c, a++) +#endif + +/* #define P DBL_MANT_DIG */ +/* Ten_pmax = floor(P*log(2)/log(5)) */ +/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ +/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ +/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ + +#define Exp_shift 20 +#define Exp_shift1 20 +#define Exp_msk1 0x100000 +#define Exp_msk11 0x100000 +#define Exp_mask 0x7ff00000 +#define P 53 +#define Nbits 53 +#define Bias 1023 +#define Emax 1023 +#define Emin (-1022) +#define Etiny (-1074) /* smallest denormal is 2**Etiny */ +#define Exp_1 0x3ff00000 +#define Exp_11 0x3ff00000 +#define Ebits 11 +#define Frac_mask 0xfffff +#define Frac_mask1 0xfffff +#define Ten_pmax 22 +#define Bletch 0x10 +#define Bndry_mask 0xfffff +#define Bndry_mask1 0xfffff +#define Sign_bit 0x80000000 +#define Log2P 1 +#define Tiny0 0 +#define Tiny1 1 +#define Quick_max 14 +#define Int_max 14 + +#ifndef Flt_Rounds +#ifdef FLT_ROUNDS +#define Flt_Rounds FLT_ROUNDS +#else +#define Flt_Rounds 1 +#endif +#endif /*Flt_Rounds*/ + +#define Rounding Flt_Rounds + +#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) +#define Big1 0xffffffff + +/* Standard NaN used by sb_stdnan. */ + +#define NAN_WORD0 0x7ff80000 +#define NAN_WORD1 0 + +/* Bits of the representation of positive infinity. */ + +#define POSINF_WORD0 0x7ff00000 +#define POSINF_WORD1 0 + +/* struct BCinfo is used to pass information from sb_strtod to bigcomp */ + +typedef struct BCinfo BCinfo; +struct +BCinfo { + int e0, nd, nd0, scale; +}; + +#define FFFFFFFF 0xffffffffUL + +#define Kmax 7 + +/* struct Bigint is used to represent arbitrary-precision integers. These + integers are stored in sign-magnitude format, with the magnitude stored as + an array of base 2**32 digits. Bigints are always normalized: if x is a + Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero. + + The Bigint fields are as follows: + + - next is a header used by Balloc and Bfree to keep track of lists + of freed Bigints; it's also used for the linked list of + powers of 5 of the form 5**2**i used by pow5mult. + - k indicates which pool this Bigint was allocated from + - maxwds is the maximum number of words space was allocated for + (usually maxwds == 2**k) + - sign is 1 for negative Bigints, 0 for positive. The sign is unused + (ignored on inputs, set to 0 on outputs) in almost all operations + involving Bigints: a notable exception is the diff function, which + ignores signs on inputs but sets the sign of the output correctly. + - wds is the actual number of significant words + - x contains the vector of words (digits) for this Bigint, from least + significant (x[0]) to most significant (x[wds-1]). +*/ + +struct +Bigint { + struct Bigint *next; + int k, maxwds, sign, wds; + ULong x[1]; +}; + +typedef struct Bigint Bigint; + +#ifndef Py_USING_MEMORY_DEBUGGER + +/* Memory management: memory is allocated from, and returned to, Kmax+1 pools + of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds == + 1 << k. These pools are maintained as linked lists, with freelist[k] + pointing to the head of the list for pool k. + + On allocation, if there's no free slot in the appropriate pool, MALLOC is + called to get more memory. This memory is not returned to the system until + Python quits. There's also a private memory pool that's allocated from + in preference to using MALLOC. + + For Bigints with more than (1 << Kmax) digits (which implies at least 1233 + decimal digits), memory is directly allocated using MALLOC, and freed using + FREE. + + XXX: it would be easy to bypass this memory-management system and + translate each call to Balloc into a call to PyMem_Malloc, and each + Bfree to PyMem_Free. Investigate whether this has any significant + performance on impact. */ + +static Bigint *freelist[Kmax+1]; + +/* Allocate space for a Bigint with up to 1<<k digits */ + +static Bigint * +Balloc(int k) +{ + int x; + Bigint *rv; + unsigned int len; + + if (k <= Kmax && (rv = freelist[k])) + freelist[k] = rv->next; + else { + x = 1 << k; + len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) + /sizeof(double); + if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { + rv = (Bigint*)pmem_next; + pmem_next += len; + } + else { + rv = (Bigint*)MALLOC(len*sizeof(double)); + if (rv == NULL) + return NULL; + } + rv->k = k; + rv->maxwds = x; + } + rv->sign = rv->wds = 0; + return rv; +} + +/* Free a Bigint allocated with Balloc */ + +static void +Bfree(Bigint *v) +{ + if (v) { + if (v->k > Kmax) + FREE((void*)v); + else { + v->next = freelist[v->k]; + freelist[v->k] = v; + } + } +} + +#else + +/* Alternative versions of Balloc and Bfree that use PyMem_Malloc and + PyMem_Free directly in place of the custom memory allocation scheme above. + These are provided for the benefit of memory debugging tools like + Valgrind. */ + +/* Allocate space for a Bigint with up to 1<<k digits */ + +static Bigint * +Balloc(int k) +{ + int x; + Bigint *rv; + unsigned int len; + + x = 1 << k; + len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) + /sizeof(double); + + rv = (Bigint*)MALLOC(len*sizeof(double)); + if (rv == NULL) + return NULL; + + rv->k = k; + rv->maxwds = x; + rv->sign = rv->wds = 0; + return rv; +} + +/* Free a Bigint allocated with Balloc */ + +static void +Bfree(Bigint *v) +{ + if (v) { + FREE((void*)v); + } +} + +#endif /* Py_USING_MEMORY_DEBUGGER */ + +#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ + y->wds*sizeof(Long) + 2*sizeof(int)) + +/* Multiply a Bigint b by m and add a. Either modifies b in place and returns + a pointer to the modified b, or Bfrees b and returns a pointer to a copy. + On failure, return NULL. In this case, b will have been already freed. */ + +static Bigint * +multadd(Bigint *b, int m, int a) /* multiply by m and add a */ +{ + int i, wds; +#ifdef ULLong + ULong *x; + ULLong carry, y; +#else + ULong carry, *x, y; + ULong xi, z; +#endif + Bigint *b1; + + wds = b->wds; + x = b->x; + i = 0; + carry = a; + do { +#ifdef ULLong + y = *x * (ULLong)m + carry; + carry = y >> 32; + *x++ = (ULong)(y & FFFFFFFF); +#else + xi = *x; + y = (xi & 0xffff) * m + carry; + z = (xi >> 16) * m + (y >> 16); + carry = z >> 16; + *x++ = (z << 16) + (y & 0xffff); +#endif + } + while(++i < wds); + if (carry) { + if (wds >= b->maxwds) { + b1 = Balloc(b->k+1); + if (b1 == NULL){ + Bfree(b); + return NULL; + } + Bcopy(b1, b); + Bfree(b); + b = b1; + } + b->x[wds++] = (ULong)carry; + b->wds = wds; + } + return b; +} + +/* convert a string s containing nd decimal digits (possibly containing a + decimal separator at position nd0, which is ignored) to a Bigint. This + function carries on where the parsing code in sb_strtod leaves off: on + entry, y9 contains the result of converting the first 9 digits. Returns + NULL on failure. */ + +static Bigint * +s2b(const char *s, int nd0, int nd, ULong y9) +{ + Bigint *b; + int i, k; + Long x, y; + + x = (nd + 8) / 9; + for(k = 0, y = 1; x > y; y <<= 1, k++) ; + b = Balloc(k); + if (b == NULL) + return NULL; + b->x[0] = y9; + b->wds = 1; + + if (nd <= 9) + return b; + + s += 9; + for (i = 9; i < nd0; i++) { + b = multadd(b, 10, *s++ - '0'); + if (b == NULL) + return NULL; + } + s++; + for(; i < nd; i++) { + b = multadd(b, 10, *s++ - '0'); + if (b == NULL) + return NULL; + } + return b; +} + +/* count leading 0 bits in the 32-bit integer x. */ + +static int +hi0bits(ULong x) +{ + int k = 0; + + if (!(x & 0xffff0000)) { + k = 16; + x <<= 16; + } + if (!(x & 0xff000000)) { + k += 8; + x <<= 8; + } + if (!(x & 0xf0000000)) { + k += 4; + x <<= 4; + } + if (!(x & 0xc0000000)) { + k += 2; + x <<= 2; + } + if (!(x & 0x80000000)) { + k++; + if (!(x & 0x40000000)) + return 32; + } + return k; +} + +/* count trailing 0 bits in the 32-bit integer y, and shift y right by that + number of bits. */ + +static int +lo0bits(ULong *y) +{ + int k; + ULong x = *y; + + if (x & 7) { + if (x & 1) + return 0; + if (x & 2) { + *y = x >> 1; + return 1; + } + *y = x >> 2; + return 2; + } + k = 0; + if (!(x & 0xffff)) { + k = 16; + x >>= 16; + } + if (!(x & 0xff)) { + k += 8; + x >>= 8; + } + if (!(x & 0xf)) { + k += 4; + x >>= 4; + } + if (!(x & 0x3)) { + k += 2; + x >>= 2; + } + if (!(x & 1)) { + k++; + x >>= 1; + if (!x) + return 32; + } + *y = x; + return k; +} + +/* convert a small nonnegative integer to a Bigint */ + +static Bigint * +i2b(int i) +{ + Bigint *b; + + b = Balloc(1); + if (b == NULL) + return NULL; + b->x[0] = i; + b->wds = 1; + return b; +} + +/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores + the signs of a and b. */ + +static Bigint * +mult(Bigint *a, Bigint *b) +{ + Bigint *c; + int k, wa, wb, wc; + ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; + ULong y; +#ifdef ULLong + ULLong carry, z; +#else + ULong carry, z; + ULong z2; +#endif + + if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) { + c = Balloc(0); + if (c == NULL) + return NULL; + c->wds = 1; + c->x[0] = 0; + return c; + } + + if (a->wds < b->wds) { + c = a; + a = b; + b = c; + } + k = a->k; + wa = a->wds; + wb = b->wds; + wc = wa + wb; + if (wc > a->maxwds) + k++; + c = Balloc(k); + if (c == NULL) + return NULL; + for(x = c->x, xa = x + wc; x < xa; x++) + *x = 0; + xa = a->x; + xae = xa + wa; + xb = b->x; + xbe = xb + wb; + xc0 = c->x; +#ifdef ULLong + for(; xb < xbe; xc0++) { + if ((y = *xb++)) { + x = xa; + xc = xc0; + carry = 0; + do { + z = *x++ * (ULLong)y + *xc + carry; + carry = z >> 32; + *xc++ = (ULong)(z & FFFFFFFF); + } + while(x < xae); + *xc = (ULong)carry; + } + } +#else + for(; xb < xbe; xb++, xc0++) { + if (y = *xb & 0xffff) { + x = xa; + xc = xc0; + carry = 0; + do { + z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; + carry = z >> 16; + z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; + carry = z2 >> 16; + Storeinc(xc, z2, z); + } + while(x < xae); + *xc = carry; + } + if (y = *xb >> 16) { + x = xa; + xc = xc0; + carry = 0; + z2 = *xc; + do { + z = (*x & 0xffff) * y + (*xc >> 16) + carry; + carry = z >> 16; + Storeinc(xc, z, z2); + z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; + carry = z2 >> 16; + } + while(x < xae); + *xc = z2; + } + } +#endif + for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; + c->wds = wc; + return c; +} + +#ifndef Py_USING_MEMORY_DEBUGGER + +/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */ + +static Bigint *p5s; + +/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on + failure; if the returned pointer is distinct from b then the original + Bigint b will have been Bfree'd. Ignores the sign of b. */ + +static Bigint * +pow5mult(Bigint *b, int k) +{ + Bigint *b1, *p5, *p51; + int i; + static int p05[3] = { 5, 25, 125 }; + + if ((i = k & 3)) { + b = multadd(b, p05[i-1], 0); + if (b == NULL) + return NULL; + } + + if (!(k >>= 2)) + return b; + p5 = p5s; + if (!p5) { + /* first time */ + p5 = i2b(625); + if (p5 == NULL) { + Bfree(b); + return NULL; + } + p5s = p5; + p5->next = 0; + } + for(;;) { + if (k & 1) { + b1 = mult(b, p5); + Bfree(b); + b = b1; + if (b == NULL) + return NULL; + } + if (!(k >>= 1)) + break; + p51 = p5->next; + if (!p51) { + p51 = mult(p5,p5); + if (p51 == NULL) { + Bfree(b); + return NULL; + } + p51->next = 0; + p5->next = p51; + } + p5 = p51; + } + return b; +} + +#else + +/* Version of pow5mult that doesn't cache powers of 5. Provided for + the benefit of memory debugging tools like Valgrind. */ + +static Bigint * +pow5mult(Bigint *b, int k) +{ + Bigint *b1, *p5, *p51; + int i; + static int p05[3] = { 5, 25, 125 }; + + if ((i = k & 3)) { + b = multadd(b, p05[i-1], 0); + if (b == NULL) + return NULL; + } + + if (!(k >>= 2)) + return b; + p5 = i2b(625); + if (p5 == NULL) { + Bfree(b); + return NULL; + } + + for(;;) { + if (k & 1) { + b1 = mult(b, p5); + Bfree(b); + b = b1; + if (b == NULL) { + Bfree(p5); + return NULL; + } + } + if (!(k >>= 1)) + break; + p51 = mult(p5, p5); + Bfree(p5); + p5 = p51; + if (p5 == NULL) { + Bfree(b); + return NULL; + } + } + Bfree(p5); + return b; +} + +#endif /* Py_USING_MEMORY_DEBUGGER */ + +/* shift a Bigint b left by k bits. Return a pointer to the shifted result, + or NULL on failure. If the returned pointer is distinct from b then the + original b will have been Bfree'd. Ignores the sign of b. */ + +static Bigint * +lshift(Bigint *b, int k) +{ + int i, k1, n, n1; + Bigint *b1; + ULong *x, *x1, *xe, z; + + if (!k || (!b->x[0] && b->wds == 1)) + return b; + + n = k >> 5; + k1 = b->k; + n1 = n + b->wds + 1; + for(i = b->maxwds; n1 > i; i <<= 1) + k1++; + b1 = Balloc(k1); + if (b1 == NULL) { + Bfree(b); + return NULL; + } + x1 = b1->x; + for(i = 0; i < n; i++) + *x1++ = 0; + x = b->x; + xe = x + b->wds; + if (k &= 0x1f) { + k1 = 32 - k; + z = 0; + do { + *x1++ = *x << k | z; + z = *x++ >> k1; + } + while(x < xe); + if ((*x1 = z)) + ++n1; + } + else do + *x1++ = *x++; + while(x < xe); + b1->wds = n1 - 1; + Bfree(b); + return b1; +} + +/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and + 1 if a > b. Ignores signs of a and b. */ + +static int +cmp(Bigint *a, Bigint *b) +{ + ULong *xa, *xa0, *xb, *xb0; + int i, j; + + i = a->wds; + j = b->wds; +#ifdef DEBUG + if (i > 1 && !a->x[i-1]) + Bug("cmp called with a->x[a->wds-1] == 0"); + if (j > 1 && !b->x[j-1]) + Bug("cmp called with b->x[b->wds-1] == 0"); +#endif + if (i -= j) + return i; + xa0 = a->x; + xa = xa0 + j; + xb0 = b->x; + xb = xb0 + j; + for(;;) { + if (*--xa != *--xb) + return *xa < *xb ? -1 : 1; + if (xa <= xa0) + break; + } + return 0; +} + +/* Take the difference of Bigints a and b, returning a new Bigint. Returns + NULL on failure. The signs of a and b are ignored, but the sign of the + result is set appropriately. */ + +static Bigint * +diff(Bigint *a, Bigint *b) +{ + Bigint *c; + int i, wa, wb; + ULong *xa, *xae, *xb, *xbe, *xc; +#ifdef ULLong + ULLong borrow, y; +#else + ULong borrow, y; + ULong z; +#endif + + i = cmp(a,b); + if (!i) { + c = Balloc(0); + if (c == NULL) + return NULL; + c->wds = 1; + c->x[0] = 0; + return c; + } + if (i < 0) { + c = a; + a = b; + b = c; + i = 1; + } + else + i = 0; + c = Balloc(a->k); + if (c == NULL) + return NULL; + c->sign = i; + wa = a->wds; + xa = a->x; + xae = xa + wa; + wb = b->wds; + xb = b->x; + xbe = xb + wb; + xc = c->x; + borrow = 0; +#ifdef ULLong + do { + y = (ULLong)*xa++ - *xb++ - borrow; + borrow = y >> 32 & (ULong)1; + *xc++ = (ULong)(y & FFFFFFFF); + } + while(xb < xbe); + while(xa < xae) { + y = *xa++ - borrow; + borrow = y >> 32 & (ULong)1; + *xc++ = (ULong)(y & FFFFFFFF); + } +#else + do { + y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(xc, z, y); + } + while(xb < xbe); + while(xa < xae) { + y = (*xa & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*xa++ >> 16) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(xc, z, y); + } +#endif + while(!*--xc) + wa--; + c->wds = wa; + return c; +} + +/* Given a positive normal double x, return the difference between x and the + next double up. Doesn't give correct results for subnormals. */ + +static double +ulp(U *x) +{ + Long L; + U u; + + L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; + word0(&u) = L; + word1(&u) = 0; + return dval(&u); +} + +/* Convert a Bigint to a double plus an exponent */ + +static double +b2d(Bigint *a, int *e) +{ + ULong *xa, *xa0, w, y, z; + int k; + U d; + + xa0 = a->x; + xa = xa0 + a->wds; + y = *--xa; +#ifdef DEBUG + if (!y) Bug("zero y in b2d"); +#endif + k = hi0bits(y); + *e = 32 - k; + if (k < Ebits) { + word0(&d) = Exp_1 | y >> (Ebits - k); + w = xa > xa0 ? *--xa : 0; + word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k); + goto ret_d; + } + z = xa > xa0 ? *--xa : 0; + if (k -= Ebits) { + word0(&d) = Exp_1 | y << k | z >> (32 - k); + y = xa > xa0 ? *--xa : 0; + word1(&d) = z << k | y >> (32 - k); + } + else { + word0(&d) = Exp_1 | y; + word1(&d) = z; + } + ret_d: + return dval(&d); +} + +/* Convert a scaled double to a Bigint plus an exponent. Similar to d2b, + except that it accepts the scale parameter used in sb_strtod (which + should be either 0 or 2*P), and the normalization for the return value is + different (see below). On input, d should be finite and nonnegative, and d + / 2**scale should be exactly representable as an IEEE 754 double. + + Returns a Bigint b and an integer e such that + + dval(d) / 2**scale = b * 2**e. + + Unlike d2b, b is not necessarily odd: b and e are normalized so + that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P + and e == Etiny. This applies equally to an input of 0.0: in that + case the return values are b = 0 and e = Etiny. + + The above normalization ensures that for all possible inputs d, + 2**e gives ulp(d/2**scale). + + Returns NULL on failure. +*/ + +static Bigint * +sd2b(U *d, int scale, int *e) +{ + Bigint *b; + + b = Balloc(1); + if (b == NULL) + return NULL; + + /* First construct b and e assuming that scale == 0. */ + b->wds = 2; + b->x[0] = word1(d); + b->x[1] = word0(d) & Frac_mask; + *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift); + if (*e < Etiny) + *e = Etiny; + else + b->x[1] |= Exp_msk1; + + /* Now adjust for scale, provided that b != 0. */ + if (scale && (b->x[0] || b->x[1])) { + *e -= scale; + if (*e < Etiny) { + scale = Etiny - *e; + *e = Etiny; + /* We can't shift more than P-1 bits without shifting out a 1. */ + assert(0 < scale && scale <= P - 1); + if (scale >= 32) { + /* The bits shifted out should all be zero. */ + assert(b->x[0] == 0); + b->x[0] = b->x[1]; + b->x[1] = 0; + scale -= 32; + } + if (scale) { + /* The bits shifted out should all be zero. */ + assert(b->x[0] << (32 - scale) == 0); + b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale)); + b->x[1] >>= scale; + } + } + } + /* Ensure b is normalized. */ + if (!b->x[1]) + b->wds = 1; + + return b; +} + +/* Convert a double to a Bigint plus an exponent. Return NULL on failure. + + Given a finite nonzero double d, return an odd Bigint b and exponent *e + such that fabs(d) = b * 2**e. On return, *bbits gives the number of + significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits). + + If d is zero, then b == 0, *e == -1010, *bbits = 0. + */ + +static Bigint * +d2b(U *d, int *e, int *bits) +{ + Bigint *b; + int de, k; + ULong *x, y, z; + int i; + + b = Balloc(1); + if (b == NULL) + return NULL; + x = b->x; + + z = word0(d) & Frac_mask; + word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */ + if ((de = (int)(word0(d) >> Exp_shift))) + z |= Exp_msk1; + if ((y = word1(d))) { + if ((k = lo0bits(&y))) { + x[0] = y | z << (32 - k); + z >>= k; + } + else + x[0] = y; + i = + b->wds = (x[1] = z) ? 2 : 1; + } + else { + k = lo0bits(&z); + x[0] = z; + i = + b->wds = 1; + k += 32; + } + if (de) { + *e = de - Bias - (P-1) + k; + *bits = P - k; + } + else { + *e = de - Bias - (P-1) + 1 + k; + *bits = 32*i - hi0bits(x[i-1]); + } + return b; +} + +/* Compute the ratio of two Bigints, as a double. The result may have an + error of up to 2.5 ulps. */ + +static double +ratio(Bigint *a, Bigint *b) +{ + U da, db; + int k, ka, kb; + + dval(&da) = b2d(a, &ka); + dval(&db) = b2d(b, &kb); + k = ka - kb + 32*(a->wds - b->wds); + if (k > 0) + word0(&da) += k*Exp_msk1; + else { + k = -k; + word0(&db) += k*Exp_msk1; + } + return dval(&da) / dval(&db); +} + +static const double +tens[] = { + 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, + 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, + 1e20, 1e21, 1e22 +}; + +static const double +bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; +static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, + 9007199254740992.*9007199254740992.e-256 + /* = 2^106 * 1e-256 */ +}; +/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ +/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ +#define Scale_Bit 0x10 +#define n_bigtens 5 + +#define ULbits 32 +#define kshift 5 +#define kmask 31 + + +static int +dshift(Bigint *b, int p2) +{ + int rv = hi0bits(b->x[b->wds-1]) - 4; + if (p2 > 0) + rv -= p2; + return rv & kmask; +} + +/* special case of Bigint division. The quotient is always in the range 0 <= + quotient < 10, and on entry the divisor S is normalized so that its top 4 + bits (28--31) are zero and bit 27 is set. */ + +static int +quorem(Bigint *b, Bigint *S) +{ + int n; + ULong *bx, *bxe, q, *sx, *sxe; +#ifdef ULLong + ULLong borrow, carry, y, ys; +#else + ULong borrow, carry, y, ys; + ULong si, z, zs; +#endif + + n = S->wds; +#ifdef DEBUG + /*debug*/ if (b->wds > n) + /*debug*/ Bug("oversize b in quorem"); +#endif + if (b->wds < n) + return 0; + sx = S->x; + sxe = sx + --n; + bx = b->x; + bxe = bx + n; + q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ +#ifdef DEBUG + /*debug*/ if (q > 9) + /*debug*/ Bug("oversized quotient in quorem"); +#endif + if (q) { + borrow = 0; + carry = 0; + do { +#ifdef ULLong + ys = *sx++ * (ULLong)q + carry; + carry = ys >> 32; + y = *bx - (ys & FFFFFFFF) - borrow; + borrow = y >> 32 & (ULong)1; + *bx++ = (ULong)(y & FFFFFFFF); +#else + si = *sx++; + ys = (si & 0xffff) * q + carry; + zs = (si >> 16) * q + (ys >> 16); + carry = zs >> 16; + y = (*bx & 0xffff) - (ys & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*bx >> 16) - (zs & 0xffff) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(bx, z, y); +#endif + } + while(sx <= sxe); + if (!*bxe) { + bx = b->x; + while(--bxe > bx && !*bxe) + --n; + b->wds = n; + } + } + if (cmp(b, S) >= 0) { + q++; + borrow = 0; + carry = 0; + bx = b->x; + sx = S->x; + do { +#ifdef ULLong + ys = *sx++ + carry; + carry = ys >> 32; + y = *bx - (ys & FFFFFFFF) - borrow; + borrow = y >> 32 & (ULong)1; + *bx++ = (ULong)(y & FFFFFFFF); +#else + si = *sx++; + ys = (si & 0xffff) + carry; + zs = (si >> 16) + (ys >> 16); + carry = zs >> 16; + y = (*bx & 0xffff) - (ys & 0xffff) - borrow; + borrow = (y & 0x10000) >> 16; + z = (*bx >> 16) - (zs & 0xffff) - borrow; + borrow = (z & 0x10000) >> 16; + Storeinc(bx, z, y); +#endif + } + while(sx <= sxe); + bx = b->x; + bxe = bx + n; + if (!*bxe) { + while(--bxe > bx && !*bxe) + --n; + b->wds = n; + } + } + return q; +} + +/* sulp(x) is a version of ulp(x) that takes bc.scale into account. + + Assuming that x is finite and nonnegative (positive zero is fine + here) and x / 2^bc.scale is exactly representable as a double, + sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */ + +static double +sulp(U *x, BCinfo *bc) +{ + U u; + + if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) { + /* rv/2^bc->scale is subnormal */ + word0(&u) = (P+2)*Exp_msk1; + word1(&u) = 0; + return u.d; + } + else { + assert(word0(x) || word1(x)); /* x != 0.0 */ + return ulp(x); + } +} + +/* The bigcomp function handles some hard cases for strtod, for inputs + with more than STRTOD_DIGLIM digits. It's called once an initial + estimate for the double corresponding to the input string has + already been obtained by the code in sb_strtod. + + The bigcomp function is only called after sb_strtod has found a + double value rv such that either rv or rv + 1ulp represents the + correctly rounded value corresponding to the original string. It + determines which of these two values is the correct one by + computing the decimal digits of rv + 0.5ulp and comparing them with + the corresponding digits of s0. + + In the following, write dv for the absolute value of the number represented + by the input string. + + Inputs: + + s0 points to the first significant digit of the input string. + + rv is a (possibly scaled) estimate for the closest double value to the + value represented by the original input to sb_strtod. If + bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to + the input value. + + bc is a struct containing information gathered during the parsing and + estimation steps of sb_strtod. Description of fields follows: + + bc->e0 gives the exponent of the input value, such that dv = (integer + given by the bd->nd digits of s0) * 10**e0 + + bc->nd gives the total number of significant digits of s0. It will + be at least 1. + + bc->nd0 gives the number of significant digits of s0 before the + decimal separator. If there's no decimal separator, bc->nd0 == + bc->nd. + + bc->scale is the value used to scale rv to avoid doing arithmetic with + subnormal values. It's either 0 or 2*P (=106). + + Outputs: + + On successful exit, rv/2^(bc->scale) is the closest double to dv. + + Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */ + +static int +bigcomp(U *rv, const char *s0, BCinfo *bc) +{ + Bigint *b, *d; + int b2, d2, dd, i, nd, nd0, odd, p2, p5; + + nd = bc->nd; + nd0 = bc->nd0; + p5 = nd + bc->e0; + b = sd2b(rv, bc->scale, &p2); + if (b == NULL) + return -1; + + /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway + case, this is used for round to even. */ + odd = b->x[0] & 1; + + /* left shift b by 1 bit and or a 1 into the least significant bit; + this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */ + b = lshift(b, 1); + if (b == NULL) + return -1; + b->x[0] |= 1; + p2--; + + p2 -= p5; + d = i2b(1); + if (d == NULL) { + Bfree(b); + return -1; + } + /* Arrange for convenient computation of quotients: + * shift left if necessary so divisor has 4 leading 0 bits. + */ + if (p5 > 0) { + d = pow5mult(d, p5); + if (d == NULL) { + Bfree(b); + return -1; + } + } + else if (p5 < 0) { + b = pow5mult(b, -p5); + if (b == NULL) { + Bfree(d); + return -1; + } + } + if (p2 > 0) { + b2 = p2; + d2 = 0; + } + else { + b2 = 0; + d2 = -p2; + } + i = dshift(d, d2); + if ((b2 += i) > 0) { + b = lshift(b, b2); + if (b == NULL) { + Bfree(d); + return -1; + } + } + if ((d2 += i) > 0) { + d = lshift(d, d2); + if (d == NULL) { + Bfree(b); + return -1; + } + } + + /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 == + * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing + * a number in the range [0.1, 1). */ + if (cmp(b, d) >= 0) + /* b/d >= 1 */ + dd = -1; + else { + i = 0; + for(;;) { + b = multadd(b, 10, 0); + if (b == NULL) { + Bfree(d); + return -1; + } + dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d); + i++; + + if (dd) + break; + if (!b->x[0] && b->wds == 1) { + /* b/d == 0 */ + dd = i < nd; + break; + } + if (!(i < nd)) { + /* b/d != 0, but digits of s0 exhausted */ + dd = -1; + break; + } + } + } + Bfree(b); + Bfree(d); + if (dd > 0 || (dd == 0 && odd)) + dval(rv) += sulp(rv, bc); + return 0; +} + +/* Return a 'standard' NaN value. + + There are exactly two quiet NaNs that don't arise by 'quieting' signaling + NaNs (see IEEE 754-2008, section 6.2.1). If sign == 0, return the one whose + sign bit is cleared. Otherwise, return the one whose sign bit is set. +*/ + +double +sb_stdnan(int sign) +{ + U rv; + word0(&rv) = NAN_WORD0; + word1(&rv) = NAN_WORD1; + if (sign) + word0(&rv) |= Sign_bit; + return dval(&rv); +} + +/* Return positive or negative infinity, according to the given sign (0 for + * positive infinity, 1 for negative infinity). */ + +double +sb_infinity(int sign) +{ + U rv; + word0(&rv) = POSINF_WORD0; + word1(&rv) = POSINF_WORD1; + return sign ? -dval(&rv) : dval(&rv); +} + +double +sb_strtod(const char *s00, char **se) +{ + int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error; + int esign, i, j, k, lz, nd, nd0, odd, sign; + const char *s, *s0, *s1; + double aadj, aadj1; + U aadj2, adj, rv, rv0; + ULong y, z, abs_exp; + Long L; + BCinfo bc; + Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; + size_t ndigits, fraclen; + + dval(&rv) = 0.; + + /* Start parsing. */ + c = *(s = s00); + + /* Parse optional sign, if present. */ + sign = 0; + switch (c) { + case '-': + sign = 1; + /* no break */ + case '+': + c = *++s; + } + + /* Skip leading zeros: lz is true iff there were leading zeros. */ + s1 = s; + while (c == '0') + c = *++s; + lz = s != s1; + + /* Point s0 at the first nonzero digit (if any). fraclen will be the + number of digits between the decimal point and the end of the + digit string. ndigits will be the total number of digits ignoring + leading zeros. */ + s0 = s1 = s; + while ('0' <= c && c <= '9') + c = *++s; + ndigits = s - s1; + fraclen = 0; + + /* Parse decimal point and following digits. */ + if (c == '.') { + c = *++s; + if (!ndigits) { + s1 = s; + while (c == '0') + c = *++s; + lz = lz || s != s1; + fraclen += (s - s1); + s0 = s; + } + s1 = s; + while ('0' <= c && c <= '9') + c = *++s; + ndigits += s - s1; + fraclen += s - s1; + } + + /* Now lz is true if and only if there were leading zero digits, and + ndigits gives the total number of digits ignoring leading zeros. A + valid input must have at least one digit. */ + if (!ndigits && !lz) { + if (se) + *se = (char *)s00; + goto parse_error; + } + + /* Range check ndigits and fraclen to make sure that they, and values + computed with them, can safely fit in an int. */ + if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) { + if (se) + *se = (char *)s00; + goto parse_error; + } + nd = (int)ndigits; + nd0 = (int)ndigits - (int)fraclen; + + /* Parse exponent. */ + e = 0; + if (c == 'e' || c == 'E') { + s00 = s; + c = *++s; + + /* Exponent sign. */ + esign = 0; + switch (c) { + case '-': + esign = 1; + /* no break */ + case '+': + c = *++s; + } + + /* Skip zeros. lz is true iff there are leading zeros. */ + s1 = s; + while (c == '0') + c = *++s; + lz = s != s1; + + /* Get absolute value of the exponent. */ + s1 = s; + abs_exp = 0; + while ('0' <= c && c <= '9') { + abs_exp = 10*abs_exp + (c - '0'); + c = *++s; + } + + /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if + there are at most 9 significant exponent digits then overflow is + impossible. */ + if (s - s1 > 9 || abs_exp > MAX_ABS_EXP) + e = (int)MAX_ABS_EXP; + else + e = (int)abs_exp; + if (esign) + e = -e; + + /* A valid exponent must have at least one digit. */ + if (s == s1 && !lz) + s = s00; + } + + /* Adjust exponent to take into account position of the point. */ + e -= nd - nd0; + if (nd0 <= 0) + nd0 = nd; + + /* Finished parsing. Set se to indicate how far we parsed */ + if (se) + *se = (char *)s; + + /* If all digits were zero, exit with return value +-0.0. Otherwise, + strip trailing zeros: scan back until we hit a nonzero digit. */ + if (!nd) + goto ret; + for (i = nd; i > 0; ) { + --i; + if (s0[i < nd0 ? i : i+1] != '0') { + ++i; + break; + } + } + e += nd - i; + nd = i; + if (nd0 > nd) + nd0 = nd; + + /* Summary of parsing results. After parsing, and dealing with zero + * inputs, we have values s0, nd0, nd, e, sign, where: + * + * - s0 points to the first significant digit of the input string + * + * - nd is the total number of significant digits (here, and + * below, 'significant digits' means the set of digits of the + * significand of the input that remain after ignoring leading + * and trailing zeros). + * + * - nd0 indicates the position of the decimal point, if present; it + * satisfies 1 <= nd0 <= nd. The nd significant digits are in + * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice + * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if + * nd0 == nd, then s0[nd0] could be any non-digit character.) + * + * - e is the adjusted exponent: the absolute value of the number + * represented by the original input string is n * 10**e, where + * n is the integer represented by the concatenation of + * s0[0:nd0] and s0[nd0+1:nd+1] + * + * - sign gives the sign of the input: 1 for negative, 0 for positive + * + * - the first and last significant digits are nonzero + */ + + /* put first DBL_DIG+1 digits into integer y and z. + * + * - y contains the value represented by the first min(9, nd) + * significant digits + * + * - if nd > 9, z contains the value represented by significant digits + * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z + * gives the value represented by the first min(16, nd) sig. digits. + */ + + bc.e0 = e1 = e; + y = z = 0; + for (i = 0; i < nd; i++) { + if (i < 9) + y = 10*y + s0[i < nd0 ? i : i+1] - '0'; + else if (i < DBL_DIG+1) + z = 10*z + s0[i < nd0 ? i : i+1] - '0'; + else + break; + } + + k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; + dval(&rv) = y; + if (k > 9) { + dval(&rv) = tens[k - 9] * dval(&rv) + z; + } + bd0 = 0; + if (nd <= DBL_DIG + && Flt_Rounds == 1 + ) { + if (!e) + goto ret; + if (e > 0) { + if (e <= Ten_pmax) { + dval(&rv) *= tens[e]; + goto ret; + } + i = DBL_DIG - nd; + if (e <= Ten_pmax + i) { + /* A fancier test would sometimes let us do + * this for larger i values. + */ + e -= i; + dval(&rv) *= tens[i]; + dval(&rv) *= tens[e]; + goto ret; + } + } + else if (e >= -Ten_pmax) { + dval(&rv) /= tens[-e]; + goto ret; + } + } + e1 += nd - k; + + bc.scale = 0; + + /* Get starting approximation = rv * 10**e1 */ + + if (e1 > 0) { + if ((i = e1 & 15)) + dval(&rv) *= tens[i]; + if (e1 &= ~15) { + if (e1 > DBL_MAX_10_EXP) + goto ovfl; + e1 >>= 4; + for(j = 0; e1 > 1; j++, e1 >>= 1) + if (e1 & 1) + dval(&rv) *= bigtens[j]; + /* The last multiplication could overflow. */ + word0(&rv) -= P*Exp_msk1; + dval(&rv) *= bigtens[j]; + if ((z = word0(&rv) & Exp_mask) + > Exp_msk1*(DBL_MAX_EXP+Bias-P)) + goto ovfl; + if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { + /* set to largest number */ + /* (Can't trust DBL_MAX) */ + word0(&rv) = Big0; + word1(&rv) = Big1; + } + else + word0(&rv) += P*Exp_msk1; + } + } + else if (e1 < 0) { + /* The input decimal value lies in [10**e1, 10**(e1+16)). + + If e1 <= -512, underflow immediately. + If e1 <= -256, set bc.scale to 2*P. + + So for input value < 1e-256, bc.scale is always set; + for input value >= 1e-240, bc.scale is never set. + For input values in [1e-256, 1e-240), bc.scale may or may + not be set. */ + + e1 = -e1; + if ((i = e1 & 15)) + dval(&rv) /= tens[i]; + if (e1 >>= 4) { + if (e1 >= 1 << n_bigtens) + goto undfl; + if (e1 & Scale_Bit) + bc.scale = 2*P; + for(j = 0; e1 > 0; j++, e1 >>= 1) + if (e1 & 1) + dval(&rv) *= tinytens[j]; + if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) + >> Exp_shift)) > 0) { + /* scaled rv is denormal; clear j low bits */ + if (j >= 32) { + word1(&rv) = 0; + if (j >= 53) + word0(&rv) = (P+2)*Exp_msk1; + else + word0(&rv) &= 0xffffffff << (j-32); + } + else + word1(&rv) &= 0xffffffff << j; + } + if (!dval(&rv)) + goto undfl; + } + } + + /* Now the hard part -- adjusting rv to the correct value.*/ + + /* Put digits into bd: true value = bd * 10^e */ + + bc.nd = nd; + bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */ + /* to silence an erroneous warning about bc.nd0 */ + /* possibly not being initialized. */ + if (nd > STRTOD_DIGLIM) { + /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */ + /* minimum number of decimal digits to distinguish double values */ + /* in IEEE arithmetic. */ + + /* Truncate input to 18 significant digits, then discard any trailing + zeros on the result by updating nd, nd0, e and y suitably. (There's + no need to update z; it's not reused beyond this point.) */ + for (i = 18; i > 0; ) { + /* scan back until we hit a nonzero digit. significant digit 'i' + is s0[i] if i < nd0, s0[i+1] if i >= nd0. */ + --i; + if (s0[i < nd0 ? i : i+1] != '0') { + ++i; + break; + } + } + e += nd - i; + nd = i; + if (nd0 > nd) + nd0 = nd; + if (nd < 9) { /* must recompute y */ + y = 0; + for(i = 0; i < nd0; ++i) + y = 10*y + s0[i] - '0'; + for(; i < nd; ++i) + y = 10*y + s0[i+1] - '0'; + } + } + bd0 = s2b(s0, nd0, nd, y); + if (bd0 == NULL) + goto failed_malloc; + + /* Notation for the comments below. Write: + + - dv for the absolute value of the number represented by the original + decimal input string. + + - if we've truncated dv, write tdv for the truncated value. + Otherwise, set tdv == dv. + + - srv for the quantity rv/2^bc.scale; so srv is the current binary + approximation to tdv (and dv). It should be exactly representable + in an IEEE 754 double. + */ + + for(;;) { + + /* This is the main correction loop for sb_strtod. + + We've got a decimal value tdv, and a floating-point approximation + srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is + close enough (i.e., within 0.5 ulps) to tdv, and to compute a new + approximation if not. + + To determine whether srv is close enough to tdv, compute integers + bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv) + respectively, and then use integer arithmetic to determine whether + |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv). + */ + + bd = Balloc(bd0->k); + if (bd == NULL) { + Bfree(bd0); + goto failed_malloc; + } + Bcopy(bd, bd0); + bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */ + if (bb == NULL) { + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + /* Record whether lsb of bb is odd, in case we need this + for the round-to-even step later. */ + odd = bb->x[0] & 1; + + /* tdv = bd * 10**e; srv = bb * 2**bbe */ + bs = i2b(1); + if (bs == NULL) { + Bfree(bb); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + + if (e >= 0) { + bb2 = bb5 = 0; + bd2 = bd5 = e; + } + else { + bb2 = bb5 = -e; + bd2 = bd5 = 0; + } + if (bbe >= 0) + bb2 += bbe; + else + bd2 -= bbe; + bs2 = bb2; + bb2++; + bd2++; + + /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1, + and bs == 1, so: + + tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5) + srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2) + 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2) + + It follows that: + + M * tdv = bd * 2**bd2 * 5**bd5 + M * srv = bb * 2**bb2 * 5**bb5 + M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5 + + for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but + this fact is not needed below.) + */ + + /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */ + i = bb2 < bd2 ? bb2 : bd2; + if (i > bs2) + i = bs2; + if (i > 0) { + bb2 -= i; + bd2 -= i; + bs2 -= i; + } + + /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */ + if (bb5 > 0) { + bs = pow5mult(bs, bb5); + if (bs == NULL) { + Bfree(bb); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + bb1 = mult(bs, bb); + Bfree(bb); + bb = bb1; + if (bb == NULL) { + Bfree(bs); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + } + if (bb2 > 0) { + bb = lshift(bb, bb2); + if (bb == NULL) { + Bfree(bs); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + } + if (bd5 > 0) { + bd = pow5mult(bd, bd5); + if (bd == NULL) { + Bfree(bb); + Bfree(bs); + Bfree(bd0); + goto failed_malloc; + } + } + if (bd2 > 0) { + bd = lshift(bd, bd2); + if (bd == NULL) { + Bfree(bb); + Bfree(bs); + Bfree(bd0); + goto failed_malloc; + } + } + if (bs2 > 0) { + bs = lshift(bs, bs2); + if (bs == NULL) { + Bfree(bb); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + } + + /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv), + respectively. Compute the difference |tdv - srv|, and compare + with 0.5 ulp(srv). */ + + delta = diff(bb, bd); + if (delta == NULL) { + Bfree(bb); + Bfree(bs); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + dsign = delta->sign; + delta->sign = 0; + i = cmp(delta, bs); + if (bc.nd > nd && i <= 0) { + if (dsign) + break; /* Must use bigcomp(). */ + + /* Here rv overestimates the truncated decimal value by at most + 0.5 ulp(rv). Hence rv either overestimates the true decimal + value by <= 0.5 ulp(rv), or underestimates it by some small + amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of + the true decimal value, so it's possible to exit. + + Exception: if scaled rv is a normal exact power of 2, but not + DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the + next double, so the correctly rounded result is either rv - 0.5 + ulp(rv) or rv; in this case, use bigcomp to distinguish. */ + + if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) { + /* rv can't be 0, since it's an overestimate for some + nonzero value. So rv is a normal power of 2. */ + j = (int)(word0(&rv) & Exp_mask) >> Exp_shift; + /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if + rv / 2^bc.scale >= 2^-1021. */ + if (j - bc.scale >= 2) { + dval(&rv) -= 0.5 * sulp(&rv, &bc); + break; /* Use bigcomp. */ + } + } + + { + bc.nd = nd; + i = -1; /* Discarded digits make delta smaller. */ + } + } + + if (i < 0) { + /* Error is less than half an ulp -- check for + * special case of mantissa a power of two. + */ + if (dsign || word1(&rv) || word0(&rv) & Bndry_mask + || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 + ) { + break; + } + if (!delta->x[0] && delta->wds <= 1) { + /* exact result */ + break; + } + delta = lshift(delta,Log2P); + if (delta == NULL) { + Bfree(bb); + Bfree(bs); + Bfree(bd); + Bfree(bd0); + goto failed_malloc; + } + if (cmp(delta, bs) > 0) + goto drop_down; + break; + } + if (i == 0) { + /* exactly half-way between */ + if (dsign) { + if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 + && word1(&rv) == ( + (bc.scale && + (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ? + (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : + 0xffffffff)) { + /*boundary case -- increment exponent*/ + word0(&rv) = (word0(&rv) & Exp_mask) + + Exp_msk1 + ; + word1(&rv) = 0; + /* dsign = 0; */ + break; + } + } + else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { + drop_down: + /* boundary case -- decrement exponent */ + if (bc.scale) { + L = word0(&rv) & Exp_mask; + if (L <= (2*P+1)*Exp_msk1) { + if (L > (P+2)*Exp_msk1) + /* round even ==> */ + /* accept rv */ + break; + /* rv = smallest denormal */ + if (bc.nd > nd) + break; + goto undfl; + } + } + L = (word0(&rv) & Exp_mask) - Exp_msk1; + word0(&rv) = L | Bndry_mask1; + word1(&rv) = 0xffffffff; + break; + } + if (!odd) + break; + if (dsign) + dval(&rv) += sulp(&rv, &bc); + else { + dval(&rv) -= sulp(&rv, &bc); + if (!dval(&rv)) { + if (bc.nd >nd) + break; + goto undfl; + } + } + /* dsign = 1 - dsign; */ + break; + } + if ((aadj = ratio(delta, bs)) <= 2.) { + if (dsign) + aadj = aadj1 = 1.; + else if (word1(&rv) || word0(&rv) & Bndry_mask) { + if (word1(&rv) == Tiny1 && !word0(&rv)) { + if (bc.nd >nd) + break; + goto undfl; + } + aadj = 1.; + aadj1 = -1.; + } + else { + /* special case -- power of FLT_RADIX to be */ + /* rounded down... */ + + if (aadj < 2./FLT_RADIX) + aadj = 1./FLT_RADIX; + else + aadj *= 0.5; + aadj1 = -aadj; + } + } + else { + aadj *= 0.5; + aadj1 = dsign ? aadj : -aadj; + if (Flt_Rounds == 0) + aadj1 += 0.5; + } + y = word0(&rv) & Exp_mask; + + /* Check for overflow */ + + if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { + dval(&rv0) = dval(&rv); + word0(&rv) -= P*Exp_msk1; + adj.d = aadj1 * ulp(&rv); + dval(&rv) += adj.d; + if ((word0(&rv) & Exp_mask) >= + Exp_msk1*(DBL_MAX_EXP+Bias-P)) { + if (word0(&rv0) == Big0 && word1(&rv0) == Big1) { + Bfree(bb); + Bfree(bd); + Bfree(bs); + Bfree(bd0); + Bfree(delta); + goto ovfl; + } + word0(&rv) = Big0; + word1(&rv) = Big1; + goto cont; + } + else + word0(&rv) += P*Exp_msk1; + } + else { + if (bc.scale && y <= 2*P*Exp_msk1) { + if (aadj <= 0x7fffffff) { + if ((z = (ULong)aadj) <= 0) + z = 1; + aadj = z; + aadj1 = dsign ? aadj : -aadj; + } + dval(&aadj2) = aadj1; + word0(&aadj2) += (2*P+1)*Exp_msk1 - y; + aadj1 = dval(&aadj2); + } + adj.d = aadj1 * ulp(&rv); + dval(&rv) += adj.d; + } + z = word0(&rv) & Exp_mask; + if (bc.nd == nd) { + if (!bc.scale) + if (y == z) { + /* Can we stop now? */ + L = (Long)aadj; + aadj -= L; + /* The tolerances below are conservative. */ + if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { + if (aadj < .4999999 || aadj > .5000001) + break; + } + else if (aadj < .4999999/FLT_RADIX) + break; + } + } + cont: + Bfree(bb); + Bfree(bd); + Bfree(bs); + Bfree(delta); + } + Bfree(bb); + Bfree(bd); + Bfree(bs); + Bfree(bd0); + Bfree(delta); + if (bc.nd > nd) { + error = bigcomp(&rv, s0, &bc); + if (error) + goto failed_malloc; + } + + if (bc.scale) { + word0(&rv0) = Exp_1 - 2*P*Exp_msk1; + word1(&rv0) = 0; + dval(&rv) *= dval(&rv0); + } + + ret: + return sign ? -dval(&rv) : dval(&rv); + + parse_error: + return 0.0; + + failed_malloc: + errno = ENOMEM; + return -1.0; + + undfl: + return sign ? -0.0 : 0.0; + + ovfl: + errno = ERANGE; + /* Can't trust HUGE_VAL */ + word0(&rv) = Exp_mask; + word1(&rv) = 0; + return sign ? -dval(&rv) : dval(&rv); + +} + +static char * +rv_alloc(int i) +{ + int j, k, *r; + + j = sizeof(ULong); + for(k = 0; + sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; + j <<= 1) + k++; + r = (int*)Balloc(k); + if (r == NULL) + return NULL; + *r = k; + return (char *)(r+1); +} + +static char * +nrv_alloc(char *s, char **rve, int n) +{ + char *rv, *t; + + rv = rv_alloc(n); + if (rv == NULL) + return NULL; + t = rv; + while((*t = *s++)) t++; + if (rve) + *rve = t; + return rv; +} + +/* freedtoa(s) must be used to free values s returned by dtoa + * when MULTIPLE_THREADS is #defined. It should be used in all cases, + * but for consistency with earlier versions of dtoa, it is optional + * when MULTIPLE_THREADS is not defined. + */ + +void +sb_freedtoa(char *s) +{ + Bigint *b = (Bigint *)((int *)s - 1); + b->maxwds = 1 << (b->k = *(int*)b); + Bfree(b); +} + +/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. + * + * Inspired by "How to Print Floating-Point Numbers Accurately" by + * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. + * + * Modifications: + * 1. Rather than iterating, we use a simple numeric overestimate + * to determine k = floor(log10(d)). We scale relevant + * quantities using O(log2(k)) rather than O(k) multiplications. + * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't + * try to generate digits strictly left to right. Instead, we + * compute with fewer bits and propagate the carry if necessary + * when rounding the final digit up. This is often faster. + * 3. Under the assumption that input will be rounded nearest, + * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. + * That is, we allow equality in stopping tests when the + * round-nearest rule will give the same floating-point value + * as would satisfaction of the stopping test with strict + * inequality. + * 4. We remove common factors of powers of 2 from relevant + * quantities. + * 5. When converting floating-point integers less than 1e16, + * we use floating-point arithmetic rather than resorting + * to multiple-precision integers. + * 6. When asked to produce fewer than 15 digits, we first try + * to get by with floating-point arithmetic; we resort to + * multiple-precision integer arithmetic only if we cannot + * guarantee that the floating-point calculation has given + * the correctly rounded result. For k requested digits and + * "uniformly" distributed input, the probability is + * something like 10^(k-15) that we must resort to the Long + * calculation. + */ + +/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory + leakage, a successful call to sb_dtoa should always be matched by a + call to sb_freedtoa. */ + +char * +sb_dtoa(double dd, int mode, int ndigits, + int *decpt, int *sign, char **rve) +{ + /* Arguments ndigits, decpt, sign are similar to those + of ecvt and fcvt; trailing zeros are suppressed from + the returned string. If not null, *rve is set to point + to the end of the return value. If d is +-Infinity or NaN, + then *decpt is set to 9999. + + mode: + 0 ==> shortest string that yields d when read in + and rounded to nearest. + 1 ==> like 0, but with Steele & White stopping rule; + e.g. with IEEE P754 arithmetic , mode 0 gives + 1e23 whereas mode 1 gives 9.999999999999999e22. + 2 ==> max(1,ndigits) significant digits. This gives a + return value similar to that of ecvt, except + that trailing zeros are suppressed. + 3 ==> through ndigits past the decimal point. This + gives a return value similar to that from fcvt, + except that trailing zeros are suppressed, and + ndigits can be negative. + 4,5 ==> similar to 2 and 3, respectively, but (in + round-nearest mode) with the tests of mode 0 to + possibly return a shorter string that rounds to d. + With IEEE arithmetic and compilation with + -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same + as modes 2 and 3 when FLT_ROUNDS != 1. + 6-9 ==> Debugging modes similar to mode - 4: don't try + fast floating-point estimate (if applicable). + + Values of mode other than 0-9 are treated as mode 0. + + Sufficient space is allocated to the return value + to hold the suppressed trailing zeros. + */ + + int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, + j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, + spec_case, try_quick; + Long L; + int denorm; + ULong x; + Bigint *b, *b1, *delta, *mlo, *mhi, *S; + U d2, eps, u; + double ds; + char *s, *s0; + + /* set pointers to NULL, to silence gcc compiler warnings and make + cleanup easier on error */ + mlo = mhi = S = 0; + s0 = 0; + + u.d = dd; + if (word0(&u) & Sign_bit) { + /* set sign for everything, including 0's and NaNs */ + *sign = 1; + word0(&u) &= ~Sign_bit; /* clear sign bit */ + } + else + *sign = 0; + + /* quick return for Infinities, NaNs and zeros */ + if ((word0(&u) & Exp_mask) == Exp_mask) + { + /* Infinity or NaN */ + *decpt = 9999; + if (!word1(&u) && !(word0(&u) & 0xfffff)) + return nrv_alloc("Infinity", rve, 8); + return nrv_alloc("NaN", rve, 3); + } + if (!dval(&u)) { + *decpt = 1; + return nrv_alloc("0", rve, 1); + } + + /* compute k = floor(log10(d)). The computation may leave k + one too large, but should never leave k too small. */ + b = d2b(&u, &be, &bbits); + if (b == NULL) + goto failed_malloc; + if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { + dval(&d2) = dval(&u); + word0(&d2) &= Frac_mask1; + word0(&d2) |= Exp_11; + + /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 + * log10(x) = log(x) / log(10) + * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) + * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) + * + * This suggests computing an approximation k to log10(d) by + * + * k = (i - Bias)*0.301029995663981 + * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); + * + * We want k to be too large rather than too small. + * The error in the first-order Taylor series approximation + * is in our favor, so we just round up the constant enough + * to compensate for any error in the multiplication of + * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, + * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, + * adding 1e-13 to the constant term more than suffices. + * Hence we adjust the constant term to 0.1760912590558. + * (We could get a more accurate k by invoking log10, + * but this is probably not worthwhile.) + */ + + i -= Bias; + denorm = 0; + } + else { + /* d is denormalized */ + + i = bbits + be + (Bias + (P-1) - 1); + x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) + : word1(&u) << (32 - i); + dval(&d2) = x; + word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ + i -= (Bias + (P-1) - 1) + 1; + denorm = 1; + } + ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + + i*0.301029995663981; + k = (int)ds; + if (ds < 0. && ds != k) + k--; /* want k = floor(ds) */ + k_check = 1; + if (k >= 0 && k <= Ten_pmax) { + if (dval(&u) < tens[k]) + k--; + k_check = 0; + } + j = bbits - i - 1; + if (j >= 0) { + b2 = 0; + s2 = j; + } + else { + b2 = -j; + s2 = 0; + } + if (k >= 0) { + b5 = 0; + s5 = k; + s2 += k; + } + else { + b2 -= k; + b5 = -k; + s5 = 0; + } + if (mode < 0 || mode > 9) + mode = 0; + + try_quick = 1; + + if (mode > 5) { + mode -= 4; + try_quick = 0; + } + leftright = 1; + ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ + /* silence erroneous "gcc -Wall" warning. */ + switch(mode) { + case 0: + case 1: + i = 18; + ndigits = 0; + break; + case 2: + leftright = 0; + /* no break */ + case 4: + if (ndigits <= 0) + ndigits = 1; + ilim = ilim1 = i = ndigits; + break; + case 3: + leftright = 0; + /* no break */ + case 5: + i = ndigits + k + 1; + ilim = i; + ilim1 = i - 1; + if (i <= 0) + i = 1; + } + s0 = rv_alloc(i); + if (s0 == NULL) + goto failed_malloc; + s = s0; + + + if (ilim >= 0 && ilim <= Quick_max && try_quick) { + + /* Try to get by with floating-point arithmetic. */ + + i = 0; + dval(&d2) = dval(&u); + k0 = k; + ilim0 = ilim; + ieps = 2; /* conservative */ + if (k > 0) { + ds = tens[k&0xf]; + j = k >> 4; + if (j & Bletch) { + /* prevent overflows */ + j &= Bletch - 1; + dval(&u) /= bigtens[n_bigtens-1]; + ieps++; + } + for(; j; j >>= 1, i++) + if (j & 1) { + ieps++; + ds *= bigtens[i]; + } + dval(&u) /= ds; + } + else if ((j1 = -k)) { + dval(&u) *= tens[j1 & 0xf]; + for(j = j1 >> 4; j; j >>= 1, i++) + if (j & 1) { + ieps++; + dval(&u) *= bigtens[i]; + } + } + if (k_check && dval(&u) < 1. && ilim > 0) { + if (ilim1 <= 0) + goto fast_failed; + ilim = ilim1; + k--; + dval(&u) *= 10.; + ieps++; + } + dval(&eps) = ieps*dval(&u) + 7.; + word0(&eps) -= (P-1)*Exp_msk1; + if (ilim == 0) { + S = mhi = 0; + dval(&u) -= 5.; + if (dval(&u) > dval(&eps)) + goto one_digit; + if (dval(&u) < -dval(&eps)) + goto no_digits; + goto fast_failed; + } + if (leftright) { + /* Use Steele & White method of only + * generating digits needed. + */ + dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); + for(i = 0;;) { + L = (Long)dval(&u); + dval(&u) -= L; + *s++ = '0' + (int)L; + if (dval(&u) < dval(&eps)) + goto ret1; + if (1. - dval(&u) < dval(&eps)) + goto bump_up; + if (++i >= ilim) + break; + dval(&eps) *= 10.; + dval(&u) *= 10.; + } + } + else { + /* Generate ilim digits, then fix them up. */ + dval(&eps) *= tens[ilim-1]; + for(i = 1;; i++, dval(&u) *= 10.) { + L = (Long)(dval(&u)); + if (!(dval(&u) -= L)) + ilim = i; + *s++ = '0' + (int)L; + if (i == ilim) { + if (dval(&u) > 0.5 + dval(&eps)) + goto bump_up; + else if (dval(&u) < 0.5 - dval(&eps)) { + while(*--s == '0'); + s++; + goto ret1; + } + break; + } + } + } + fast_failed: + s = s0; + dval(&u) = dval(&d2); + k = k0; + ilim = ilim0; + } + + /* Do we have a "small" integer? */ + + if (be >= 0 && k <= Int_max) { + /* Yes. */ + ds = tens[k]; + if (ndigits < 0 && ilim <= 0) { + S = mhi = 0; + if (ilim < 0 || dval(&u) <= 5*ds) + goto no_digits; + goto one_digit; + } + for(i = 1;; i++, dval(&u) *= 10.) { + L = (Long)(dval(&u) / ds); + dval(&u) -= L*ds; + *s++ = '0' + (int)L; + if (!dval(&u)) { + break; + } + if (i == ilim) { + dval(&u) += dval(&u); + if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { + bump_up: + while(*--s == '9') + if (s == s0) { + k++; + *s = '0'; + break; + } + ++*s++; + } + break; + } + } + goto ret1; + } + + m2 = b2; + m5 = b5; + if (leftright) { + i = + denorm ? be + (Bias + (P-1) - 1 + 1) : + 1 + P - bbits; + b2 += i; + s2 += i; + mhi = i2b(1); + if (mhi == NULL) + goto failed_malloc; + } + if (m2 > 0 && s2 > 0) { + i = m2 < s2 ? m2 : s2; + b2 -= i; + m2 -= i; + s2 -= i; + } + if (b5 > 0) { + if (leftright) { + if (m5 > 0) { + mhi = pow5mult(mhi, m5); + if (mhi == NULL) + goto failed_malloc; + b1 = mult(mhi, b); + Bfree(b); + b = b1; + if (b == NULL) + goto failed_malloc; + } + if ((j = b5 - m5)) { + b = pow5mult(b, j); + if (b == NULL) + goto failed_malloc; + } + } + else { + b = pow5mult(b, b5); + if (b == NULL) + goto failed_malloc; + } + } + S = i2b(1); + if (S == NULL) + goto failed_malloc; + if (s5 > 0) { + S = pow5mult(S, s5); + if (S == NULL) + goto failed_malloc; + } + + /* Check for special case that d is a normalized power of 2. */ + + spec_case = 0; + if ((mode < 2 || leftright) + ) { + if (!word1(&u) && !(word0(&u) & Bndry_mask) + && word0(&u) & (Exp_mask & ~Exp_msk1) + ) { + /* The special case */ + b2 += Log2P; + s2 += Log2P; + spec_case = 1; + } + } + + /* Arrange for convenient computation of quotients: + * shift left if necessary so divisor has 4 leading 0 bits. + * + * Perhaps we should just compute leading 28 bits of S once + * and for all and pass them and a shift to quorem, so it + * can do shifts and ors to compute the numerator for q. + */ +#define iInc 28 + i = dshift(S, s2); + b2 += i; + m2 += i; + s2 += i; + if (b2 > 0) { + b = lshift(b, b2); + if (b == NULL) + goto failed_malloc; + } + if (s2 > 0) { + S = lshift(S, s2); + if (S == NULL) + goto failed_malloc; + } + if (k_check) { + if (cmp(b,S) < 0) { + k--; + b = multadd(b, 10, 0); /* we botched the k estimate */ + if (b == NULL) + goto failed_malloc; + if (leftright) { + mhi = multadd(mhi, 10, 0); + if (mhi == NULL) + goto failed_malloc; + } + ilim = ilim1; + } + } + if (ilim <= 0 && (mode == 3 || mode == 5)) { + if (ilim < 0) { + /* no digits, fcvt style */ + no_digits: + k = -1 - ndigits; + goto ret; + } + else { + S = multadd(S, 5, 0); + if (S == NULL) + goto failed_malloc; + if (cmp(b, S) <= 0) + goto no_digits; + } + one_digit: + *s++ = '1'; + k++; + goto ret; + } + if (leftright) { + if (m2 > 0) { + mhi = lshift(mhi, m2); + if (mhi == NULL) + goto failed_malloc; + } + + /* Compute mlo -- check for special case + * that d is a normalized power of 2. + */ + + mlo = mhi; + if (spec_case) { + mhi = Balloc(mhi->k); + if (mhi == NULL) + goto failed_malloc; + Bcopy(mhi, mlo); + mhi = lshift(mhi, Log2P); + if (mhi == NULL) + goto failed_malloc; + } + + for(i = 1;;i++) { + dig = quorem(b,S) + '0'; + /* Do we yet have the shortest decimal string + * that will round to d? + */ + j = cmp(b, mlo); + delta = diff(S, mhi); + if (delta == NULL) + goto failed_malloc; + j1 = delta->sign ? 1 : cmp(b, delta); + Bfree(delta); + if (j1 == 0 && mode != 1 && !(word1(&u) & 1) + ) { + if (dig == '9') + goto round_9_up; + if (j > 0) + dig++; + *s++ = dig; + goto ret; + } + if (j < 0 || (j == 0 && mode != 1 + && !(word1(&u) & 1) + )) { + if (!b->x[0] && b->wds <= 1) { + goto accept_dig; + } + if (j1 > 0) { + b = lshift(b, 1); + if (b == NULL) + goto failed_malloc; + j1 = cmp(b, S); + if ((j1 > 0 || (j1 == 0 && dig & 1)) + && dig++ == '9') + goto round_9_up; + } + accept_dig: + *s++ = dig; + goto ret; + } + if (j1 > 0) { + if (dig == '9') { /* possible if i == 1 */ + round_9_up: + *s++ = '9'; + goto roundoff; + } + *s++ = dig + 1; + goto ret; + } + *s++ = dig; + if (i == ilim) + break; + b = multadd(b, 10, 0); + if (b == NULL) + goto failed_malloc; + if (mlo == mhi) { + mlo = mhi = multadd(mhi, 10, 0); + if (mlo == NULL) + goto failed_malloc; + } + else { + mlo = multadd(mlo, 10, 0); + if (mlo == NULL) + goto failed_malloc; + mhi = multadd(mhi, 10, 0); + if (mhi == NULL) + goto failed_malloc; + } + } + } + else + for(i = 1;; i++) { + *s++ = dig = quorem(b,S) + '0'; + if (!b->x[0] && b->wds <= 1) { + goto ret; + } + if (i >= ilim) + break; + b = multadd(b, 10, 0); + if (b == NULL) + goto failed_malloc; + } + + /* Round off last digit */ + + b = lshift(b, 1); + if (b == NULL) + goto failed_malloc; + j = cmp(b, S); + if (j > 0 || (j == 0 && dig & 1)) { + roundoff: + while(*--s == '9') + if (s == s0) { + k++; + *s++ = '1'; + goto ret; + } + ++*s++; + } + else { + while(*--s == '0'); + s++; + } + ret: + Bfree(S); + if (mhi) { + if (mlo && mlo != mhi) + Bfree(mlo); + Bfree(mhi); + } + ret1: + Bfree(b); + *s = 0; + *decpt = k + 1; + if (rve) + *rve = s; + return s0; + failed_malloc: + if (S) + Bfree(S); + if (mlo && mlo != mhi) + Bfree(mlo); + if (mhi) + Bfree(mhi); + if (b) + Bfree(b); + if (s0) + sb_freedtoa(s0); + return NULL; +} +#ifdef __cplusplus +} +#endif |