/****************************************************************

 *

 * 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.

 *

 ***************************************************************/

/* Please send bug reports to

    David M. Gay

    Bell Laboratories, Room 2C-463

    600 Mountain Avenue

    Murray Hill, NJ 07974-0636

    U.S.A.

    dmg@bell-labs.com

 */

/* 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).

 */

/*

 * #define IEEE_8087 for IEEE-arithmetic machines where the least

 *  significant byte has the lowest address.

 * #define IEEE_MC68k for IEEE-arithmetic machines where the most

 *  significant byte has the lowest address.

 * #define Long int on machines with 32-bit ints and 64-bit longs.

 * #define IBM for IBM mainframe-style floating-point arithmetic.

 * #define VAX for VAX-style floating-point arithmetic (D_floating).

 * #define No_leftright to omit left-right logic in fast floating-point

 *  computation of dtoa.

 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3

 *  and strtod and dtoa should round accordingly.

 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3

 *  and Honor_FLT_ROUNDS is not #defined.

 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines

 *  that use extended-precision instructions to compute rounded

 *  products and quotients) with IBM.

 * #define ROUND_BIASED for IEEE-format with biased rounding.

 * #define Inaccurate_Divide for IEEE-format with correctly rounded

 *  products but inaccurate quotients, e.g., for Intel i860.

 * #define NO_LONG_LONG on machines that do not have a "long long"

 *  integer type (of >= 64 bits).  On such machines, you can

 *  #define Just_16 to store 16 bits per 32-bit Long when doing

 *  high-precision integer arithmetic.  Whether this speeds things

 *  up or slows things down depends on the machine and the number

 *  being converted.  If long long is available and the name is

 *  something other than "long long", #define Llong to be the name,

 *  and if "unsigned Llong" does not work as an unsigned version of

 *  Llong, #define #ULLong to be the corresponding unsigned type.

 * #define KR_headers for old-style C function headers.

 * #define Bad_float_h if your system lacks a float.h or if it does not

 *  define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,

 *  FLT_RADIX, FLT_ROUNDS, and DBL_MAX.

 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)

 *  if memory is available and otherwise does something you deem

 *  appropriate.  If MALLOC is undefined, malloc will be invoked

 *  directly -- and assumed always to succeed.

 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making

 *  memory allocations from a private pool of memory when possible.

 *  When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,

 *  unless #defined to be a different length.  This default length

 *  suffices to get rid of MALLOC calls except for unusual cases,

 *  such as decimal-to-binary conversion of a very long string of

 *  digits.  The longest string dtoa can return is about 751 bytes

 *  long.  For conversions by strtod of strings of 800 digits and

 *  all dtoa conversions in single-threaded executions with 8-byte

 *  pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte

 *  pointers, PRIVATE_MEM >= 7112 appears adequate.

 * #define INFNAN_CHECK on IEEE systems to cause strtod to check for

 *  Infinity and NaN (case insensitively).  On some systems (e.g.,

 *  some HP systems), it may be necessary to #define NAN_WORD0

 *  appropriately -- to the most significant word of a quiet NaN.

 *  (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)

 *  When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,

 *  strtod also accepts (case insensitively) strings of the form

 *  NaN(x), where x is a string of hexadecimal digits and spaces;

 *  if there is only one string of hexadecimal digits, it is taken

 *  for the 52 fraction bits of the resulting NaN; if there are two

 *  or more strings of hex digits, the first is for the high 20 bits,

 *  the second and subsequent for the low 32 bits, with intervening

 *  white space ignored; but if this results in none of the 52

 *  fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0

 *  and NAN_WORD1 are used instead.

 * #define MULTIPLE_THREADS if the system offers preemptively scheduled

 *  multiple threads.  In this case, you must provide (or suitably

 *  #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed

 *  by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed

 *  in pow5mult, ensures lazy evaluation of only one copy of high

 *  powers of 5; omitting this lock would introduce a small

 *  probability of wasting memory, but would otherwise be harmless.)

 *  You must also invoke freedtoa(s) to free the value s returned by

 *  dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.

 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that

 *  avoids underflows on inputs whose result does not underflow.

 *  If you #define NO_IEEE_Scale on a machine that uses IEEE-format

 *  floating-point numbers and flushes underflows to zero rather

 *  than implementing gradual underflow, then you must also #define

 *  Sudden_Underflow.

 * #define YES_ALIAS to permit aliasing certain double values with

 *  arrays of ULongs.  This leads to slightly better code with

 *  some compilers and was always used prior to 19990916, but it

 *  is not strictly legal and can cause trouble with aggressively

 *  optimizing compilers (e.g., gcc 2.95.1 under -O2).

 * #define USE_LOCALE to use the current locale's decimal_point value.

 * #define SET_INEXACT if IEEE arithmetic is being used and extra

 *  computation should be done to set the inexact flag when the

 *  result is inexact and avoid setting inexact when the result

 *  is exact.  In this case, dtoa.c must be compiled in

 *  an environment, perhaps provided by #include "dtoa.c" in a

 *  suitable wrapper, that defines two functions,

 *      int get_inexact(void);

 *      void clear_inexact(void);

 *  such that get_inexact() returns a nonzero value if the

 *  inexact bit is already set, and clear_inexact() sets the

 *  inexact bit to 0.  When SET_INEXACT is #defined, strtod

 *  also does extra computations to set the underflow and overflow

 *  flags when appropriate (i.e., when the result is tiny and

 *  inexact or when it is a numeric value rounded to +-infinity).

 * #define NO_ERRNO if strtod should not assign errno = ERANGE when

 *  the result overflows to +-Infinity or underflows to 0.

 */

#include "config.h"

#include "dtoa.h"

#if COMPILER(MSVC)

#pragma warning(disable: 4244)

#pragma warning(disable: 4245)

#pragma warning(disable: 4554)

#endif

#if PLATFORM(BIG_ENDIAN)

#define IEEE_MC68k

#else

#define IEEE_8087

#endif

#define INFNAN_CHECK

#ifndef Long

#define Long int

#endif

#ifndef ULong

typedef unsigned Long ULong;

#endif

#ifdef DEBUG

#include <stdio.h>

#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}

#endif

#include <stdlib.h>

#include <string.h>

#ifdef USE_LOCALE

#include <locale.h>

#endif

#ifdef MALLOC

#ifdef KR_headers

extern char *MALLOC();

#else

extern void *MALLOC(size_t);

#endif

#else

#define MALLOC malloc

#endif

#ifndef Omit_Private_Memory

#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;

#endif

#undef IEEE_Arith

#undef Avoid_Underflow

#ifdef IEEE_MC68k

#define IEEE_Arith

#endif

#ifdef IEEE_8087

#define IEEE_Arith

#endif

#include <errno.h>

#ifdef Bad_float_h

#ifdef IEEE_Arith

#define DBL_DIG 15

#define DBL_MAX_10_EXP 308

#define DBL_MAX_EXP 1024

#define FLT_RADIX 2

#endif /*IEEE_Arith*/

#ifdef IBM

#define DBL_DIG 16

#define DBL_MAX_10_EXP 75

#define DBL_MAX_EXP 63

#define FLT_RADIX 16

#define DBL_MAX 7.2370055773322621e+75

#endif

#ifdef VAX

#define DBL_DIG 16

#define DBL_MAX_10_EXP 38

#define DBL_MAX_EXP 127

#define FLT_RADIX 2

#define DBL_MAX 1.7014118346046923e+38

#endif

#ifndef LONG_MAX

#define LONG_MAX 2147483647

#endif

#else /* ifndef Bad_float_h */

#include <float.h>

#endif /* Bad_float_h */

#ifndef __MATH_H__

#include <math.h>

#endif

#define strtod kjs_strtod

#define dtoa kjs_dtoa

#define freedtoa kjs_freedtoa

#ifdef __cplusplus

extern "C" {

#endif

#ifndef CONST_

#ifdef KR_headers

#define CONST_ /* blank */

#else

#define CONST_ const

#endif

#endif

#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1

Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.

#endif

typedef union { double d; ULong L[2]; } U;

#ifdef YES_ALIAS

#define dval(x) x

#ifdef IEEE_8087

#define word0(x) ((ULong *)&x)[1]

#define word1(x) ((ULong *)&x)[0]

#else

#define word0(x) ((ULong *)&x)[0]

#define word1(x) ((ULong *)&x)[1]

#endif

#else

#ifdef IEEE_8087

#define word0(x) ((U*)&x)->L[1]

#define word1(x) ((U*)&x)->L[0]

#else

#define word0(x) ((U*)&x)->L[0]

#define word1(x) ((U*)&x)->L[1]

#endif

#define dval(x) ((U*)&x)->d

#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) + defined(VAX)

#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) */

#ifdef IEEE_Arith

#define Exp_shift  20

#define Exp_shift1 20

#define Exp_msk1    0x100000

#define Exp_msk11   0x100000

#define Exp_mask  0x7ff00000

#define P 53

#define Bias 1023

#define Emin (-1022)

#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 LSB 1

#define Sign_bit 0x80000000

#define Log2P 1

#define Tiny0 0

#define Tiny1 1

#define Quick_max 14

#define Int_max 14

#ifndef NO_IEEE_Scale

#define Avoid_Underflow

#ifdef Flush_Denorm /* debugging option */

#undef Sudden_Underflow

#endif

#endif

#ifndef Flt_Rounds

#ifdef FLT_ROUNDS

#define Flt_Rounds FLT_ROUNDS

#else

#define Flt_Rounds 1

#endif

#endif /*Flt_Rounds*/

#ifdef Honor_FLT_ROUNDS

#define Rounding rounding

#undef Check_FLT_ROUNDS

#define Check_FLT_ROUNDS

#else

#define Rounding Flt_Rounds

#endif

#else /* ifndef IEEE_Arith */

#undef Check_FLT_ROUNDS

#undef Honor_FLT_ROUNDS

#undef SET_INEXACT

#undef  Sudden_Underflow

#define Sudden_Underflow

#ifdef IBM

#undef Flt_Rounds

#define Flt_Rounds 0

#define Exp_shift  24

#define Exp_shift1 24

#define Exp_msk1   0x1000000

#define Exp_msk11  0x1000000

#define Exp_mask  0x7f000000

#define P 14

#define Bias 65

#define Exp_1  0x41000000

#define Exp_11 0x41000000

#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */

#define Frac_mask  0xffffff

#define Frac_mask1 0xffffff

#define Bletch 4

#define Ten_pmax 22

#define Bndry_mask  0xefffff

#define Bndry_mask1 0xffffff

#define LSB 1

#define Sign_bit 0x80000000

#define Log2P 4

#define Tiny0 0x100000

#define Tiny1 0

#define Quick_max 14

#define Int_max 15

#else /* VAX */

#undef Flt_Rounds

#define Flt_Rounds 1

#define Exp_shift  23

#define Exp_shift1 7

#define Exp_msk1    0x80

#define Exp_msk11   0x800000

#define Exp_mask  0x7f80

#define P 56

#define Bias 129

#define Exp_1  0x40800000

#define Exp_11 0x4080

#define Ebits 8

#define Frac_mask  0x7fffff

#define Frac_mask1 0xffff007f

#define Ten_pmax 24

#define Bletch 2

#define Bndry_mask  0xffff007f

#define Bndry_mask1 0xffff007f

#define LSB 0x10000

#define Sign_bit 0x8000

#define Log2P 1

#define Tiny0 0x80

#define Tiny1 0

#define Quick_max 15

#define Int_max 15

#endif /* IBM, VAX */

#endif /* IEEE_Arith */

#ifndef IEEE_Arith

#define ROUND_BIASED

#endif

#ifdef RND_PRODQUOT

#define rounded_product(a,b) a = rnd_prod(a, b)

#define rounded_quotient(a,b) a = rnd_quot(a, b)

#ifdef KR_headers

extern double rnd_prod(), rnd_quot();

#else

extern double rnd_prod(double, double), rnd_quot(double, double);

#endif

#else

#define rounded_product(a,b) a *= b

#define rounded_quotient(a,b) a /= b

#endif

#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))

#define Big1 0xffffffff

#ifndef Pack_32

#define Pack_32

#endif

#ifdef KR_headers

#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)

#else

#define FFFFFFFF 0xffffffffUL

#endif

#ifdef NO_LONG_LONG

#undef ULLong

#ifdef Just_16

#undef Pack_32

/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.

 * This makes some inner loops simpler and sometimes saves work

 * during multiplications, but it often seems to make things slightly

 * slower.  Hence the default is now to store 32 bits per Long.

 */

#endif

#else   /* long long available */

#ifndef Llong

#define Llong long long

#endif

#ifndef ULLong

#define ULLong unsigned Llong

#endif

#endif /* NO_LONG_LONG */

#ifndef MULTIPLE_THREADS

#define ACQUIRE_DTOA_LOCK(n)    /*nothing*/

#define FREE_DTOA_LOCK(n)   /*nothing*/

#endif

#define Kmax 15

 struct

Bigint {

    struct Bigint *next;

    int k, maxwds, sign, wds;

    ULong x[1];

    };

 typedef struct Bigint Bigint;

 static Bigint *freelist[Kmax+1];

 static Bigint *

Balloc

#ifdef KR_headers

    (k) int k;

#else

    (int k)

#endif

{

    int x;

    Bigint *rv;

#ifndef Omit_Private_Memory

    unsigned int len;

#endif

    ACQUIRE_DTOA_LOCK(0);

    if ((rv = freelist[k])) {

        freelist[k] = rv->next;

        }