Changeset 31288 in webkit for trunk/JavaScriptCore/kjs/dtoa.cpp


Ignore:
Timestamp:
Mar 25, 2008, 12:43:15 PM (17 years ago)
Author:
[email protected]
Message:

Rubber-stamped by Darin.

Cleanup dtoa.cpp style.

  • kjs/dtoa.cpp: (Bigint::Balloc): (Bigint::Bfree): (Bigint::multadd): (Bigint::s2b): (Bigint::hi0bits): (Bigint::lo0bits): (Bigint::i2b): (Bigint::mult): (Bigint::pow5mult): (Bigint::lshift): (Bigint::cmp): (Bigint::diff): (Bigint::ulp): (Bigint::b2d): (Bigint::d2b): (Bigint::ratio): (Bigint::): (Bigint::match): (Bigint::hexnan): (Bigint::kjs_strtod): (Bigint::quorem): (Bigint::rv_alloc): (Bigint::nrv_alloc): (Bigint::kjs_freedtoa): (Bigint::kjs_dtoa):
  • kjs/dtoa.h:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/JavaScriptCore/kjs/dtoa.cpp

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