Back to index

lightning-sunbird  0.9+nobinonly
prdtoa.c
Go to the documentation of this file.
00001 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
00002 /* ***** BEGIN LICENSE BLOCK *****
00003  * Version: MPL 1.1/GPL 2.0/LGPL 2.1
00004  *
00005  * The contents of this file are subject to the Mozilla Public License Version
00006  * 1.1 (the "License"); you may not use this file except in compliance with
00007  * the License. You may obtain a copy of the License at
00008  * http://www.mozilla.org/MPL/
00009  *
00010  * Software distributed under the License is distributed on an "AS IS" basis,
00011  * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
00012  * for the specific language governing rights and limitations under the
00013  * License.
00014  *
00015  * The Original Code is the Netscape Portable Runtime (NSPR).
00016  *
00017  * The Initial Developer of the Original Code is
00018  * Netscape Communications Corporation.
00019  * Portions created by the Initial Developer are Copyright (C) 1998-2000
00020  * the Initial Developer. All Rights Reserved.
00021  *
00022  * Contributor(s):
00023  *
00024  * Alternatively, the contents of this file may be used under the terms of
00025  * either the GNU General Public License Version 2 or later (the "GPL"), or
00026  * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
00027  * in which case the provisions of the GPL or the LGPL are applicable instead
00028  * of those above. If you wish to allow use of your version of this file only
00029  * under the terms of either the GPL or the LGPL, and not to allow others to
00030  * use your version of this file under the terms of the MPL, indicate your
00031  * decision by deleting the provisions above and replace them with the notice
00032  * and other provisions required by the GPL or the LGPL. If you do not delete
00033  * the provisions above, a recipient may use your version of this file under
00034  * the terms of any one of the MPL, the GPL or the LGPL.
00035  *
00036  * ***** END LICENSE BLOCK ***** */
00037 
00038 #include "primpl.h"
00039 
00040 #define MULTIPLE_THREADS
00041 #define ACQUIRE_DTOA_LOCK(n)       PR_Lock(dtoa_lock[n])
00042 #define FREE_DTOA_LOCK(n)   PR_Unlock(dtoa_lock[n])
00043 
00044 static PRLock *dtoa_lock[2];
00045 
00046 void _PR_InitDtoa(void)
00047 {
00048     dtoa_lock[0] = PR_NewLock();
00049     dtoa_lock[1] = PR_NewLock();
00050 }
00051 
00052 void _PR_CleanupDtoa(void)
00053 {
00054     PR_DestroyLock(dtoa_lock[0]);
00055     dtoa_lock[0] = NULL;
00056     PR_DestroyLock(dtoa_lock[1]);
00057     dtoa_lock[1] = NULL;
00058 
00059     /* FIXME: deal with freelist and p5s. */
00060 }
00061 
00062 #if defined(__arm) || defined(__arm__) || defined(__arm26__) \
00063     || defined(__arm32__)
00064 #define IEEE_ARM
00065 #elif defined(IS_LITTLE_ENDIAN)
00066 #define IEEE_8087
00067 #else
00068 #define IEEE_MC68k
00069 #endif
00070 
00071 #define Long PRInt32
00072 #define ULong PRUint32
00073 #define NO_LONG_LONG
00074 
00075 /****************************************************************
00076  *
00077  * The author of this software is David M. Gay.
00078  *
00079  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
00080  *
00081  * Permission to use, copy, modify, and distribute this software for any
00082  * purpose without fee is hereby granted, provided that this entire notice
00083  * is included in all copies of any software which is or includes a copy
00084  * or modification of this software and in all copies of the supporting
00085  * documentation for such software.
00086  *
00087  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00088  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
00089  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00090  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00091  *
00092  ***************************************************************/
00093 
00094 /* Please send bug reports to David M. Gay (dmg at acm dot org,
00095  * with " at " changed at "@" and " dot " changed to ".").     */
00096 
00097 /* On a machine with IEEE extended-precision registers, it is
00098  * necessary to specify double-precision (53-bit) rounding precision
00099  * before invoking strtod or dtoa.  If the machine uses (the equivalent
00100  * of) Intel 80x87 arithmetic, the call
00101  *     _control87(PC_53, MCW_PC);
00102  * does this with many compilers.  Whether this or another call is
00103  * appropriate depends on the compiler; for this to work, it may be
00104  * necessary to #include "float.h" or another system-dependent header
00105  * file.
00106  */
00107 
00108 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
00109  *
00110  * This strtod returns a nearest machine number to the input decimal
00111  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
00112  * broken by the IEEE round-even rule.  Otherwise ties are broken by
00113  * biased rounding (add half and chop).
00114  *
00115  * Inspired loosely by William D. Clinger's paper "How to Read Floating
00116  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
00117  *
00118  * Modifications:
00119  *
00120  *     1. We only require IEEE, IBM, or VAX double-precision
00121  *            arithmetic (not IEEE double-extended).
00122  *     2. We get by with floating-point arithmetic in a case that
00123  *            Clinger missed -- when we're computing d * 10^n
00124  *            for a small integer d and the integer n is not too
00125  *            much larger than 22 (the maximum integer k for which
00126  *            we can represent 10^k exactly), we may be able to
00127  *            compute (d*10^k) * 10^(e-k) with just one roundoff.
00128  *     3. Rather than a bit-at-a-time adjustment of the binary
00129  *            result in the hard case, we use floating-point
00130  *            arithmetic to determine the adjustment to within
00131  *            one bit; only in really hard cases do we need to
00132  *            compute a second residual.
00133  *     4. Because of 3., we don't need a large table of powers of 10
00134  *            for ten-to-e (just some small tables, e.g. of 10^k
00135  *            for 0 <= k <= 22).
00136  */
00137 
00138 /*
00139  * #define IEEE_8087 for IEEE-arithmetic machines where the least
00140  *     significant byte has the lowest address.
00141  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
00142  *     significant byte has the lowest address.
00143  * #define IEEE_ARM for IEEE-arithmetic machines where the two words
00144  *     in a double are stored in big endian order but the two shorts
00145  *     in a word are still stored in little endian order.
00146  * #define Long int on machines with 32-bit ints and 64-bit longs.
00147  * #define IBM for IBM mainframe-style floating-point arithmetic.
00148  * #define VAX for VAX-style floating-point arithmetic (D_floating).
00149  * #define No_leftright to omit left-right logic in fast floating-point
00150  *     computation of dtoa.
00151  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00152  *     and strtod and dtoa should round accordingly.
00153  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
00154  *     and Honor_FLT_ROUNDS is not #defined.
00155  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
00156  *     that use extended-precision instructions to compute rounded
00157  *     products and quotients) with IBM.
00158  * #define ROUND_BIASED for IEEE-format with biased rounding.
00159  * #define Inaccurate_Divide for IEEE-format with correctly rounded
00160  *     products but inaccurate quotients, e.g., for Intel i860.
00161  * #define NO_LONG_LONG on machines that do not have a "long long"
00162  *     integer type (of >= 64 bits).  On such machines, you can
00163  *     #define Just_16 to store 16 bits per 32-bit Long when doing
00164  *     high-precision integer arithmetic.  Whether this speeds things
00165  *     up or slows things down depends on the machine and the number
00166  *     being converted.  If long long is available and the name is
00167  *     something other than "long long", #define Llong to be the name,
00168  *     and if "unsigned Llong" does not work as an unsigned version of
00169  *     Llong, #define #ULLong to be the corresponding unsigned type.
00170  * #define KR_headers for old-style C function headers.
00171  * #define Bad_float_h if your system lacks a float.h or if it does not
00172  *     define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
00173  *     FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
00174  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
00175  *     if memory is available and otherwise does something you deem
00176  *     appropriate.  If MALLOC is undefined, malloc will be invoked
00177  *     directly -- and assumed always to succeed.
00178  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
00179  *     memory allocations from a private pool of memory when possible.
00180  *     When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
00181  *     unless #defined to be a different length.  This default length
00182  *     suffices to get rid of MALLOC calls except for unusual cases,
00183  *     such as decimal-to-binary conversion of a very long string of
00184  *     digits.  The longest string dtoa can return is about 751 bytes
00185  *     long.  For conversions by strtod of strings of 800 digits and
00186  *     all dtoa conversions in single-threaded executions with 8-byte
00187  *     pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
00188  *     pointers, PRIVATE_MEM >= 7112 appears adequate.
00189  * #define INFNAN_CHECK on IEEE systems to cause strtod to check for
00190  *     Infinity and NaN (case insensitively).  On some systems (e.g.,
00191  *     some HP systems), it may be necessary to #define NAN_WORD0
00192  *     appropriately -- to the most significant word of a quiet NaN.
00193  *     (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
00194  *     When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
00195  *     strtod also accepts (case insensitively) strings of the form
00196  *     NaN(x), where x is a string of hexadecimal digits and spaces;
00197  *     if there is only one string of hexadecimal digits, it is taken
00198  *     for the 52 fraction bits of the resulting NaN; if there are two
00199  *     or more strings of hex digits, the first is for the high 20 bits,
00200  *     the second and subsequent for the low 32 bits, with intervening
00201  *     white space ignored; but if this results in none of the 52
00202  *     fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
00203  *     and NAN_WORD1 are used instead.
00204  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
00205  *     multiple threads.  In this case, you must provide (or suitably
00206  *     #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
00207  *     by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
00208  *     in pow5mult, ensures lazy evaluation of only one copy of high
00209  *     powers of 5; omitting this lock would introduce a small
00210  *     probability of wasting memory, but would otherwise be harmless.)
00211  *     You must also invoke freedtoa(s) to free the value s returned by
00212  *     dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
00213  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
00214  *     avoids underflows on inputs whose result does not underflow.
00215  *     If you #define NO_IEEE_Scale on a machine that uses IEEE-format
00216  *     floating-point numbers and flushes underflows to zero rather
00217  *     than implementing gradual underflow, then you must also #define
00218  *     Sudden_Underflow.
00219  * #define YES_ALIAS to permit aliasing certain double values with
00220  *     arrays of ULongs.  This leads to slightly better code with
00221  *     some compilers and was always used prior to 19990916, but it
00222  *     is not strictly legal and can cause trouble with aggressively
00223  *     optimizing compilers (e.g., gcc 2.95.1 under -O2).
00224  * #define USE_LOCALE to use the current locale's decimal_point value.
00225  * #define SET_INEXACT if IEEE arithmetic is being used and extra
00226  *     computation should be done to set the inexact flag when the
00227  *     result is inexact and avoid setting inexact when the result
00228  *     is exact.  In this case, dtoa.c must be compiled in
00229  *     an environment, perhaps provided by #include "dtoa.c" in a
00230  *     suitable wrapper, that defines two functions,
00231  *            int get_inexact(void);
00232  *            void clear_inexact(void);
00233  *     such that get_inexact() returns a nonzero value if the
00234  *     inexact bit is already set, and clear_inexact() sets the
00235  *     inexact bit to 0.  When SET_INEXACT is #defined, strtod
00236  *     also does extra computations to set the underflow and overflow
00237  *     flags when appropriate (i.e., when the result is tiny and
00238  *     inexact or when it is a numeric value rounded to +-infinity).
00239  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
00240  *     the result overflows to +-Infinity or underflows to 0.
00241  */
00242 
00243 #ifndef Long
00244 #define Long long
00245 #endif
00246 #ifndef ULong
00247 typedef unsigned Long ULong;
00248 #endif
00249 
00250 #ifdef DEBUG
00251 #include "stdio.h"
00252 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
00253 #endif
00254 
00255 #include "stdlib.h"
00256 #include "string.h"
00257 
00258 #ifdef USE_LOCALE
00259 #include "locale.h"
00260 #endif
00261 
00262 #ifdef MALLOC
00263 #ifdef KR_headers
00264 extern char *MALLOC();
00265 #else
00266 extern void *MALLOC(size_t);
00267 #endif
00268 #else
00269 #define MALLOC malloc
00270 #endif
00271 
00272 #ifndef Omit_Private_Memory
00273 #ifndef PRIVATE_MEM
00274 #define PRIVATE_MEM 2304
00275 #endif
00276 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
00277 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
00278 #endif
00279 
00280 #undef IEEE_Arith
00281 #undef Avoid_Underflow
00282 #ifdef IEEE_MC68k
00283 #define IEEE_Arith
00284 #endif
00285 #ifdef IEEE_8087
00286 #define IEEE_Arith
00287 #endif
00288 #ifdef IEEE_ARM
00289 #define IEEE_Arith
00290 #endif
00291 
00292 #include "errno.h"
00293 
00294 #ifdef Bad_float_h
00295 
00296 #ifdef IEEE_Arith
00297 #define DBL_DIG 15
00298 #define DBL_MAX_10_EXP 308
00299 #define DBL_MAX_EXP 1024
00300 #define FLT_RADIX 2
00301 #endif /*IEEE_Arith*/
00302 
00303 #ifdef IBM
00304 #define DBL_DIG 16
00305 #define DBL_MAX_10_EXP 75
00306 #define DBL_MAX_EXP 63
00307 #define FLT_RADIX 16
00308 #define DBL_MAX 7.2370055773322621e+75
00309 #endif
00310 
00311 #ifdef VAX
00312 #define DBL_DIG 16
00313 #define DBL_MAX_10_EXP 38
00314 #define DBL_MAX_EXP 127
00315 #define FLT_RADIX 2
00316 #define DBL_MAX 1.7014118346046923e+38
00317 #endif
00318 
00319 #ifndef LONG_MAX
00320 #define LONG_MAX 2147483647
00321 #endif
00322 
00323 #else /* ifndef Bad_float_h */
00324 #include "float.h"
00325 /*
00326  * MacOS 10.2 defines the macro FLT_ROUNDS to an internal function
00327  * which does not exist on 10.1.  We can safely #define it to 1 here
00328  * to allow 10.2 builds to run on 10.1, since we can't use fesetround()
00329  * (which does not exist on 10.1 either).
00330  */
00331 #if defined(XP_MACOSX) && (!defined(MAC_OS_X_VERSION_10_2) || \
00332     MAC_OS_X_VERSION_MIN_REQUIRED < MAC_OS_X_VERSION_10_2)
00333 #undef FLT_ROUNDS
00334 #define FLT_ROUNDS 1
00335 #endif /* DT < 10.2 */
00336 #endif /* Bad_float_h */
00337 
00338 #ifndef __MATH_H__
00339 #include "math.h"
00340 #endif
00341 
00342 #ifdef __cplusplus
00343 extern "C" {
00344 #endif
00345 
00346 #ifndef CONST
00347 #ifdef KR_headers
00348 #define CONST /* blank */
00349 #else
00350 #define CONST const
00351 #endif
00352 #endif
00353 
00354 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_ARM) + defined(VAX) + defined(IBM) != 1
00355 Exactly one of IEEE_8087, IEEE_MC68k, IEEE_ARM, VAX, or IBM should be defined.
00356 #endif
00357 
00358 typedef union { double d; ULong L[2]; } U;
00359 
00360 #ifdef YES_ALIAS
00361 #define dval(x) x
00362 #ifdef IEEE_8087
00363 #define word0(x) ((ULong *)&x)[1]
00364 #define word1(x) ((ULong *)&x)[0]
00365 #else
00366 #define word0(x) ((ULong *)&x)[0]
00367 #define word1(x) ((ULong *)&x)[1]
00368 #endif
00369 #else
00370 #ifdef IEEE_8087
00371 #define word0(x) ((U*)&x)->L[1]
00372 #define word1(x) ((U*)&x)->L[0]
00373 #else
00374 #define word0(x) ((U*)&x)->L[0]
00375 #define word1(x) ((U*)&x)->L[1]
00376 #endif
00377 #define dval(x) ((U*)&x)->d
00378 #endif
00379 
00380 /* The following definition of Storeinc is appropriate for MIPS processors.
00381  * An alternative that might be better on some machines is
00382  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
00383  */
00384 #if defined(IEEE_8087) + defined(IEEE_ARM) + defined(VAX)
00385 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
00386 ((unsigned short *)a)[0] = (unsigned short)c, a++)
00387 #else
00388 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
00389 ((unsigned short *)a)[1] = (unsigned short)c, a++)
00390 #endif
00391 
00392 /* #define P DBL_MANT_DIG */
00393 /* Ten_pmax = floor(P*log(2)/log(5)) */
00394 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
00395 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
00396 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
00397 
00398 #ifdef IEEE_Arith
00399 #define Exp_shift  20
00400 #define Exp_shift1 20
00401 #define Exp_msk1    0x100000
00402 #define Exp_msk11   0x100000
00403 #define Exp_mask  0x7ff00000
00404 #define P 53
00405 #define Bias 1023
00406 #define Emin (-1022)
00407 #define Exp_1  0x3ff00000
00408 #define Exp_11 0x3ff00000
00409 #define Ebits 11
00410 #define Frac_mask  0xfffff
00411 #define Frac_mask1 0xfffff
00412 #define Ten_pmax 22
00413 #define Bletch 0x10
00414 #define Bndry_mask  0xfffff
00415 #define Bndry_mask1 0xfffff
00416 #define LSB 1
00417 #define Sign_bit 0x80000000
00418 #define Log2P 1
00419 #define Tiny0 0
00420 #define Tiny1 1
00421 #define Quick_max 14
00422 #define Int_max 14
00423 #ifndef NO_IEEE_Scale
00424 #define Avoid_Underflow
00425 #ifdef Flush_Denorm  /* debugging option */
00426 #undef Sudden_Underflow
00427 #endif
00428 #endif
00429 
00430 #ifndef Flt_Rounds
00431 #ifdef FLT_ROUNDS
00432 #define Flt_Rounds FLT_ROUNDS
00433 #else
00434 #define Flt_Rounds 1
00435 #endif
00436 #endif /*Flt_Rounds*/
00437 
00438 #ifdef Honor_FLT_ROUNDS
00439 #define Rounding rounding
00440 #undef Check_FLT_ROUNDS
00441 #define Check_FLT_ROUNDS
00442 #else
00443 #define Rounding Flt_Rounds
00444 #endif
00445 
00446 #else /* ifndef IEEE_Arith */
00447 #undef Check_FLT_ROUNDS
00448 #undef Honor_FLT_ROUNDS
00449 #undef SET_INEXACT
00450 #undef  Sudden_Underflow
00451 #define Sudden_Underflow
00452 #ifdef IBM
00453 #undef Flt_Rounds
00454 #define Flt_Rounds 0
00455 #define Exp_shift  24
00456 #define Exp_shift1 24
00457 #define Exp_msk1   0x1000000
00458 #define Exp_msk11  0x1000000
00459 #define Exp_mask  0x7f000000
00460 #define P 14
00461 #define Bias 65
00462 #define Exp_1  0x41000000
00463 #define Exp_11 0x41000000
00464 #define Ebits 8      /* exponent has 7 bits, but 8 is the right value in b2d */
00465 #define Frac_mask  0xffffff
00466 #define Frac_mask1 0xffffff
00467 #define Bletch 4
00468 #define Ten_pmax 22
00469 #define Bndry_mask  0xefffff
00470 #define Bndry_mask1 0xffffff
00471 #define LSB 1
00472 #define Sign_bit 0x80000000
00473 #define Log2P 4
00474 #define Tiny0 0x100000
00475 #define Tiny1 0
00476 #define Quick_max 14
00477 #define Int_max 15
00478 #else /* VAX */
00479 #undef Flt_Rounds
00480 #define Flt_Rounds 1
00481 #define Exp_shift  23
00482 #define Exp_shift1 7
00483 #define Exp_msk1    0x80
00484 #define Exp_msk11   0x800000
00485 #define Exp_mask  0x7f80
00486 #define P 56
00487 #define Bias 129
00488 #define Exp_1  0x40800000
00489 #define Exp_11 0x4080
00490 #define Ebits 8
00491 #define Frac_mask  0x7fffff
00492 #define Frac_mask1 0xffff007f
00493 #define Ten_pmax 24
00494 #define Bletch 2
00495 #define Bndry_mask  0xffff007f
00496 #define Bndry_mask1 0xffff007f
00497 #define LSB 0x10000
00498 #define Sign_bit 0x8000
00499 #define Log2P 1
00500 #define Tiny0 0x80
00501 #define Tiny1 0
00502 #define Quick_max 15
00503 #define Int_max 15
00504 #endif /* IBM, VAX */
00505 #endif /* IEEE_Arith */
00506 
00507 #ifndef IEEE_Arith
00508 #define ROUND_BIASED
00509 #endif
00510 
00511 #ifdef RND_PRODQUOT
00512 #define rounded_product(a,b) a = rnd_prod(a, b)
00513 #define rounded_quotient(a,b) a = rnd_quot(a, b)
00514 #ifdef KR_headers
00515 extern double rnd_prod(), rnd_quot();
00516 #else
00517 extern double rnd_prod(double, double), rnd_quot(double, double);
00518 #endif
00519 #else
00520 #define rounded_product(a,b) a *= b
00521 #define rounded_quotient(a,b) a /= b
00522 #endif
00523 
00524 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
00525 #define Big1 0xffffffff
00526 
00527 #ifndef Pack_32
00528 #define Pack_32
00529 #endif
00530 
00531 #ifdef KR_headers
00532 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
00533 #else
00534 #define FFFFFFFF 0xffffffffUL
00535 #endif
00536 
00537 #ifdef NO_LONG_LONG
00538 #undef ULLong
00539 #ifdef Just_16
00540 #undef Pack_32
00541 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
00542  * This makes some inner loops simpler and sometimes saves work
00543  * during multiplications, but it often seems to make things slightly
00544  * slower.  Hence the default is now to store 32 bits per Long.
00545  */
00546 #endif
00547 #else  /* long long available */
00548 #ifndef Llong
00549 #define Llong long long
00550 #endif
00551 #ifndef ULLong
00552 #define ULLong unsigned Llong
00553 #endif
00554 #endif /* NO_LONG_LONG */
00555 
00556 #ifndef MULTIPLE_THREADS
00557 #define ACQUIRE_DTOA_LOCK(n)       /*nothing*/
00558 #define FREE_DTOA_LOCK(n)   /*nothing*/
00559 #endif
00560 
00561 #define Kmax 15
00562 
00563  struct
00564 Bigint {
00565        struct Bigint *next;
00566        int k, maxwds, sign, wds;
00567        ULong x[1];
00568        };
00569 
00570  typedef struct Bigint Bigint;
00571 
00572  static Bigint *freelist[Kmax+1];
00573 
00574  static Bigint *
00575 Balloc
00576 #ifdef KR_headers
00577        (k) int k;
00578 #else
00579        (int k)
00580 #endif
00581 {
00582        int x;
00583        Bigint *rv;
00584 #ifndef Omit_Private_Memory
00585        unsigned int len;
00586 #endif
00587 
00588        ACQUIRE_DTOA_LOCK(0);
00589        if (rv = freelist[k]) {
00590               freelist[k] = rv->next;
00591               }
00592        else {
00593               x = 1 << k;
00594 #ifdef Omit_Private_Memory
00595               rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
00596 #else
00597               len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
00598                      /sizeof(double);
00599               if (pmem_next - private_mem + len <= PRIVATE_mem) {
00600                      rv = (Bigint*)pmem_next;
00601                      pmem_next += len;
00602                      }
00603               else
00604                      rv = (Bigint*)MALLOC(len*sizeof(double));
00605 #endif
00606               rv->k = k;
00607               rv->maxwds = x;
00608               }
00609        FREE_DTOA_LOCK(0);
00610        rv->sign = rv->wds = 0;
00611        return rv;
00612        }
00613 
00614  static void
00615 Bfree
00616 #ifdef KR_headers
00617        (v) Bigint *v;
00618 #else
00619        (Bigint *v)
00620 #endif
00621 {
00622        if (v) {
00623               ACQUIRE_DTOA_LOCK(0);
00624               v->next = freelist[v->k];
00625               freelist[v->k] = v;
00626               FREE_DTOA_LOCK(0);
00627               }
00628        }
00629 
00630 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
00631 y->wds*sizeof(Long) + 2*sizeof(int))
00632 
00633  static Bigint *
00634 multadd
00635 #ifdef KR_headers
00636        (b, m, a) Bigint *b; int m, a;
00637 #else
00638        (Bigint *b, int m, int a)   /* multiply by m and add a */
00639 #endif
00640 {
00641        int i, wds;
00642 #ifdef ULLong
00643        ULong *x;
00644        ULLong carry, y;
00645 #else
00646        ULong carry, *x, y;
00647 #ifdef Pack_32
00648        ULong xi, z;
00649 #endif
00650 #endif
00651        Bigint *b1;
00652 
00653        wds = b->wds;
00654        x = b->x;
00655        i = 0;
00656        carry = a;
00657        do {
00658 #ifdef ULLong
00659               y = *x * (ULLong)m + carry;
00660               carry = y >> 32;
00661               *x++ = y & FFFFFFFF;
00662 #else
00663 #ifdef Pack_32
00664               xi = *x;
00665               y = (xi & 0xffff) * m + carry;
00666               z = (xi >> 16) * m + (y >> 16);
00667               carry = z >> 16;
00668               *x++ = (z << 16) + (y & 0xffff);
00669 #else
00670               y = *x * m + carry;
00671               carry = y >> 16;
00672               *x++ = y & 0xffff;
00673 #endif
00674 #endif
00675               }
00676               while(++i < wds);
00677        if (carry) {
00678               if (wds >= b->maxwds) {
00679                      b1 = Balloc(b->k+1);
00680                      Bcopy(b1, b);
00681                      Bfree(b);
00682                      b = b1;
00683                      }
00684               b->x[wds++] = carry;
00685               b->wds = wds;
00686               }
00687        return b;
00688        }
00689 
00690  static Bigint *
00691 s2b
00692 #ifdef KR_headers
00693        (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
00694 #else
00695        (CONST char *s, int nd0, int nd, ULong y9)
00696 #endif
00697 {
00698        Bigint *b;
00699        int i, k;
00700        Long x, y;
00701 
00702        x = (nd + 8) / 9;
00703        for(k = 0, y = 1; x > y; y <<= 1, k++) ;
00704 #ifdef Pack_32
00705        b = Balloc(k);
00706        b->x[0] = y9;
00707        b->wds = 1;
00708 #else
00709        b = Balloc(k+1);
00710        b->x[0] = y9 & 0xffff;
00711        b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
00712 #endif
00713 
00714        i = 9;
00715        if (9 < nd0) {
00716               s += 9;
00717               do b = multadd(b, 10, *s++ - '0');
00718                      while(++i < nd0);
00719               s++;
00720               }
00721        else
00722               s += 10;
00723        for(; i < nd; i++)
00724               b = multadd(b, 10, *s++ - '0');
00725        return b;
00726        }
00727 
00728  static int
00729 hi0bits
00730 #ifdef KR_headers
00731        (x) register ULong x;
00732 #else
00733        (register ULong x)
00734 #endif
00735 {
00736        register int k = 0;
00737 
00738        if (!(x & 0xffff0000)) {
00739               k = 16;
00740               x <<= 16;
00741               }
00742        if (!(x & 0xff000000)) {
00743               k += 8;
00744               x <<= 8;
00745               }
00746        if (!(x & 0xf0000000)) {
00747               k += 4;
00748               x <<= 4;
00749               }
00750        if (!(x & 0xc0000000)) {
00751               k += 2;
00752               x <<= 2;
00753               }
00754        if (!(x & 0x80000000)) {
00755               k++;
00756               if (!(x & 0x40000000))
00757                      return 32;
00758               }
00759        return k;
00760        }
00761 
00762  static int
00763 lo0bits
00764 #ifdef KR_headers
00765        (y) ULong *y;
00766 #else
00767        (ULong *y)
00768 #endif
00769 {
00770        register int k;
00771        register ULong x = *y;
00772 
00773        if (x & 7) {
00774               if (x & 1)
00775                      return 0;
00776               if (x & 2) {
00777                      *y = x >> 1;
00778                      return 1;
00779                      }
00780               *y = x >> 2;
00781               return 2;
00782               }
00783        k = 0;
00784        if (!(x & 0xffff)) {
00785               k = 16;
00786               x >>= 16;
00787               }
00788        if (!(x & 0xff)) {
00789               k += 8;
00790               x >>= 8;
00791               }
00792        if (!(x & 0xf)) {
00793               k += 4;
00794               x >>= 4;
00795               }
00796        if (!(x & 0x3)) {
00797               k += 2;
00798               x >>= 2;
00799               }
00800        if (!(x & 1)) {
00801               k++;
00802               x >>= 1;
00803               if (!x)
00804                      return 32;
00805               }
00806        *y = x;
00807        return k;
00808        }
00809 
00810  static Bigint *
00811 i2b
00812 #ifdef KR_headers
00813        (i) int i;
00814 #else
00815        (int i)
00816 #endif
00817 {
00818        Bigint *b;
00819 
00820        b = Balloc(1);
00821        b->x[0] = i;
00822        b->wds = 1;
00823        return b;
00824        }
00825 
00826  static Bigint *
00827 mult
00828 #ifdef KR_headers
00829        (a, b) Bigint *a, *b;
00830 #else
00831        (Bigint *a, Bigint *b)
00832 #endif
00833 {
00834        Bigint *c;
00835        int k, wa, wb, wc;
00836        ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
00837        ULong y;
00838 #ifdef ULLong
00839        ULLong carry, z;
00840 #else
00841        ULong carry, z;
00842 #ifdef Pack_32
00843        ULong z2;
00844 #endif
00845 #endif
00846 
00847        if (a->wds < b->wds) {
00848               c = a;
00849               a = b;
00850               b = c;
00851               }
00852        k = a->k;
00853        wa = a->wds;
00854        wb = b->wds;
00855        wc = wa + wb;
00856        if (wc > a->maxwds)
00857               k++;
00858        c = Balloc(k);
00859        for(x = c->x, xa = x + wc; x < xa; x++)
00860               *x = 0;
00861        xa = a->x;
00862        xae = xa + wa;
00863        xb = b->x;
00864        xbe = xb + wb;
00865        xc0 = c->x;
00866 #ifdef ULLong
00867        for(; xb < xbe; xc0++) {
00868               if (y = *xb++) {
00869                      x = xa;
00870                      xc = xc0;
00871                      carry = 0;
00872                      do {
00873                             z = *x++ * (ULLong)y + *xc + carry;
00874                             carry = z >> 32;
00875                             *xc++ = z & FFFFFFFF;
00876                             }
00877                             while(x < xae);
00878                      *xc = carry;
00879                      }
00880               }
00881 #else
00882 #ifdef Pack_32
00883        for(; xb < xbe; xb++, xc0++) {
00884               if (y = *xb & 0xffff) {
00885                      x = xa;
00886                      xc = xc0;
00887                      carry = 0;
00888                      do {
00889                             z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
00890                             carry = z >> 16;
00891                             z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
00892                             carry = z2 >> 16;
00893                             Storeinc(xc, z2, z);
00894                             }
00895                             while(x < xae);
00896                      *xc = carry;
00897                      }
00898               if (y = *xb >> 16) {
00899                      x = xa;
00900                      xc = xc0;
00901                      carry = 0;
00902                      z2 = *xc;
00903                      do {
00904                             z = (*x & 0xffff) * y + (*xc >> 16) + carry;
00905                             carry = z >> 16;
00906                             Storeinc(xc, z, z2);
00907                             z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
00908                             carry = z2 >> 16;
00909                             }
00910                             while(x < xae);
00911                      *xc = z2;
00912                      }
00913               }
00914 #else
00915        for(; xb < xbe; xc0++) {
00916               if (y = *xb++) {
00917                      x = xa;
00918                      xc = xc0;
00919                      carry = 0;
00920                      do {
00921                             z = *x++ * y + *xc + carry;
00922                             carry = z >> 16;
00923                             *xc++ = z & 0xffff;
00924                             }
00925                             while(x < xae);
00926                      *xc = carry;
00927                      }
00928               }
00929 #endif
00930 #endif
00931        for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
00932        c->wds = wc;
00933        return c;
00934        }
00935 
00936  static Bigint *p5s;
00937 
00938  static Bigint *
00939 pow5mult
00940 #ifdef KR_headers
00941        (b, k) Bigint *b; int k;
00942 #else
00943        (Bigint *b, int k)
00944 #endif
00945 {
00946        Bigint *b1, *p5, *p51;
00947        int i;
00948        static int p05[3] = { 5, 25, 125 };
00949 
00950        if (i = k & 3)
00951               b = multadd(b, p05[i-1], 0);
00952 
00953        if (!(k >>= 2))
00954               return b;
00955        if (!(p5 = p5s)) {
00956               /* first time */
00957 #ifdef MULTIPLE_THREADS
00958               ACQUIRE_DTOA_LOCK(1);
00959               if (!(p5 = p5s)) {
00960                      p5 = p5s = i2b(625);
00961                      p5->next = 0;
00962                      }
00963               FREE_DTOA_LOCK(1);
00964 #else
00965               p5 = p5s = i2b(625);
00966               p5->next = 0;
00967 #endif
00968               }
00969        for(;;) {
00970               if (k & 1) {
00971                      b1 = mult(b, p5);
00972                      Bfree(b);
00973                      b = b1;
00974                      }
00975               if (!(k >>= 1))
00976                      break;
00977               if (!(p51 = p5->next)) {
00978 #ifdef MULTIPLE_THREADS
00979                      ACQUIRE_DTOA_LOCK(1);
00980                      if (!(p51 = p5->next)) {
00981                             p51 = p5->next = mult(p5,p5);
00982                             p51->next = 0;
00983                             }
00984                      FREE_DTOA_LOCK(1);
00985 #else
00986                      p51 = p5->next = mult(p5,p5);
00987                      p51->next = 0;
00988 #endif
00989                      }
00990               p5 = p51;
00991               }
00992        return b;
00993        }
00994 
00995  static Bigint *
00996 lshift
00997 #ifdef KR_headers
00998        (b, k) Bigint *b; int k;
00999 #else
01000        (Bigint *b, int k)
01001 #endif
01002 {
01003        int i, k1, n, n1;
01004        Bigint *b1;
01005        ULong *x, *x1, *xe, z;
01006 
01007 #ifdef Pack_32
01008        n = k >> 5;
01009 #else
01010        n = k >> 4;
01011 #endif
01012        k1 = b->k;
01013        n1 = n + b->wds + 1;
01014        for(i = b->maxwds; n1 > i; i <<= 1)
01015               k1++;
01016        b1 = Balloc(k1);
01017        x1 = b1->x;
01018        for(i = 0; i < n; i++)
01019               *x1++ = 0;
01020        x = b->x;
01021        xe = x + b->wds;
01022 #ifdef Pack_32
01023        if (k &= 0x1f) {
01024               k1 = 32 - k;
01025               z = 0;
01026               do {
01027                      *x1++ = *x << k | z;
01028                      z = *x++ >> k1;
01029                      }
01030                      while(x < xe);
01031               if (*x1 = z)
01032                      ++n1;
01033               }
01034 #else
01035        if (k &= 0xf) {
01036               k1 = 16 - k;
01037               z = 0;
01038               do {
01039                      *x1++ = *x << k  & 0xffff | z;
01040                      z = *x++ >> k1;
01041                      }
01042                      while(x < xe);
01043               if (*x1 = z)
01044                      ++n1;
01045               }
01046 #endif
01047        else do
01048               *x1++ = *x++;
01049               while(x < xe);
01050        b1->wds = n1 - 1;
01051        Bfree(b);
01052        return b1;
01053        }
01054 
01055  static int
01056 cmp
01057 #ifdef KR_headers
01058        (a, b) Bigint *a, *b;
01059 #else
01060        (Bigint *a, Bigint *b)
01061 #endif
01062 {
01063        ULong *xa, *xa0, *xb, *xb0;
01064        int i, j;
01065 
01066        i = a->wds;
01067        j = b->wds;
01068 #ifdef DEBUG
01069        if (i > 1 && !a->x[i-1])
01070               Bug("cmp called with a->x[a->wds-1] == 0");
01071        if (j > 1 && !b->x[j-1])
01072               Bug("cmp called with b->x[b->wds-1] == 0");
01073 #endif
01074        if (i -= j)
01075               return i;
01076        xa0 = a->x;
01077        xa = xa0 + j;
01078        xb0 = b->x;
01079        xb = xb0 + j;
01080        for(;;) {
01081               if (*--xa != *--xb)
01082                      return *xa < *xb ? -1 : 1;
01083               if (xa <= xa0)
01084                      break;
01085               }
01086        return 0;
01087        }
01088 
01089  static Bigint *
01090 diff
01091 #ifdef KR_headers
01092        (a, b) Bigint *a, *b;
01093 #else
01094        (Bigint *a, Bigint *b)
01095 #endif
01096 {
01097        Bigint *c;
01098        int i, wa, wb;
01099        ULong *xa, *xae, *xb, *xbe, *xc;
01100 #ifdef ULLong
01101        ULLong borrow, y;
01102 #else
01103        ULong borrow, y;
01104 #ifdef Pack_32
01105        ULong z;
01106 #endif
01107 #endif
01108 
01109        i = cmp(a,b);
01110        if (!i) {
01111               c = Balloc(0);
01112               c->wds = 1;
01113               c->x[0] = 0;
01114               return c;
01115               }
01116        if (i < 0) {
01117               c = a;
01118               a = b;
01119               b = c;
01120               i = 1;
01121               }
01122        else
01123               i = 0;
01124        c = Balloc(a->k);
01125        c->sign = i;
01126        wa = a->wds;
01127        xa = a->x;
01128        xae = xa + wa;
01129        wb = b->wds;
01130        xb = b->x;
01131        xbe = xb + wb;
01132        xc = c->x;
01133        borrow = 0;
01134 #ifdef ULLong
01135        do {
01136               y = (ULLong)*xa++ - *xb++ - borrow;
01137               borrow = y >> 32 & (ULong)1;
01138               *xc++ = y & FFFFFFFF;
01139               }
01140               while(xb < xbe);
01141        while(xa < xae) {
01142               y = *xa++ - borrow;
01143               borrow = y >> 32 & (ULong)1;
01144               *xc++ = y & FFFFFFFF;
01145               }
01146 #else
01147 #ifdef Pack_32
01148        do {
01149               y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
01150               borrow = (y & 0x10000) >> 16;
01151               z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
01152               borrow = (z & 0x10000) >> 16;
01153               Storeinc(xc, z, y);
01154               }
01155               while(xb < xbe);
01156        while(xa < xae) {
01157               y = (*xa & 0xffff) - borrow;
01158               borrow = (y & 0x10000) >> 16;
01159               z = (*xa++ >> 16) - borrow;
01160               borrow = (z & 0x10000) >> 16;
01161               Storeinc(xc, z, y);
01162               }
01163 #else
01164        do {
01165               y = *xa++ - *xb++ - borrow;
01166               borrow = (y & 0x10000) >> 16;
01167               *xc++ = y & 0xffff;
01168               }
01169               while(xb < xbe);
01170        while(xa < xae) {
01171               y = *xa++ - borrow;
01172               borrow = (y & 0x10000) >> 16;
01173               *xc++ = y & 0xffff;
01174               }
01175 #endif
01176 #endif
01177        while(!*--xc)
01178               wa--;
01179        c->wds = wa;
01180        return c;
01181        }
01182 
01183  static double
01184 ulp
01185 #ifdef KR_headers
01186        (x) double x;
01187 #else
01188        (double x)
01189 #endif
01190 {
01191        register Long L;
01192        double a;
01193 
01194        L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
01195 #ifndef Avoid_Underflow
01196 #ifndef Sudden_Underflow
01197        if (L > 0) {
01198 #endif
01199 #endif
01200 #ifdef IBM
01201               L |= Exp_msk1 >> 4;
01202 #endif
01203               word0(a) = L;
01204               word1(a) = 0;
01205 #ifndef Avoid_Underflow
01206 #ifndef Sudden_Underflow
01207               }
01208        else {
01209               L = -L >> Exp_shift;
01210               if (L < Exp_shift) {
01211                      word0(a) = 0x80000 >> L;
01212                      word1(a) = 0;
01213                      }
01214               else {
01215                      word0(a) = 0;
01216                      L -= Exp_shift;
01217                      word1(a) = L >= 31 ? 1 : 1 << 31 - L;
01218                      }
01219               }
01220 #endif
01221 #endif
01222        return dval(a);
01223        }
01224 
01225  static double
01226 b2d
01227 #ifdef KR_headers
01228        (a, e) Bigint *a; int *e;
01229 #else
01230        (Bigint *a, int *e)
01231 #endif
01232 {
01233        ULong *xa, *xa0, w, y, z;
01234        int k;
01235        double d;
01236 #ifdef VAX
01237        ULong d0, d1;
01238 #else
01239 #define d0 word0(d)
01240 #define d1 word1(d)
01241 #endif
01242 
01243        xa0 = a->x;
01244        xa = xa0 + a->wds;
01245        y = *--xa;
01246 #ifdef DEBUG
01247        if (!y) Bug("zero y in b2d");
01248 #endif
01249        k = hi0bits(y);
01250        *e = 32 - k;
01251 #ifdef Pack_32
01252        if (k < Ebits) {
01253               d0 = Exp_1 | y >> Ebits - k;
01254               w = xa > xa0 ? *--xa : 0;
01255               d1 = y << (32-Ebits) + k | w >> Ebits - k;
01256               goto ret_d;
01257               }
01258        z = xa > xa0 ? *--xa : 0;
01259        if (k -= Ebits) {
01260               d0 = Exp_1 | y << k | z >> 32 - k;
01261               y = xa > xa0 ? *--xa : 0;
01262               d1 = z << k | y >> 32 - k;
01263               }
01264        else {
01265               d0 = Exp_1 | y;
01266               d1 = z;
01267               }
01268 #else
01269        if (k < Ebits + 16) {
01270               z = xa > xa0 ? *--xa : 0;
01271               d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
01272               w = xa > xa0 ? *--xa : 0;
01273               y = xa > xa0 ? *--xa : 0;
01274               d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
01275               goto ret_d;
01276               }
01277        z = xa > xa0 ? *--xa : 0;
01278        w = xa > xa0 ? *--xa : 0;
01279        k -= Ebits + 16;
01280        d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
01281        y = xa > xa0 ? *--xa : 0;
01282        d1 = w << k + 16 | y << k;
01283 #endif
01284  ret_d:
01285 #ifdef VAX
01286        word0(d) = d0 >> 16 | d0 << 16;
01287        word1(d) = d1 >> 16 | d1 << 16;
01288 #else
01289 #undef d0
01290 #undef d1
01291 #endif
01292        return dval(d);
01293        }
01294 
01295  static Bigint *
01296 d2b
01297 #ifdef KR_headers
01298        (d, e, bits) double d; int *e, *bits;
01299 #else
01300        (double d, int *e, int *bits)
01301 #endif
01302 {
01303        Bigint *b;
01304        int de, k;
01305        ULong *x, y, z;
01306 #ifndef Sudden_Underflow
01307        int i;
01308 #endif
01309 #ifdef VAX
01310        ULong d0, d1;
01311        d0 = word0(d) >> 16 | word0(d) << 16;
01312        d1 = word1(d) >> 16 | word1(d) << 16;
01313 #else
01314 #define d0 word0(d)
01315 #define d1 word1(d)
01316 #endif
01317 
01318 #ifdef Pack_32
01319        b = Balloc(1);
01320 #else
01321        b = Balloc(2);
01322 #endif
01323        x = b->x;
01324 
01325        z = d0 & Frac_mask;
01326        d0 &= 0x7fffffff;    /* clear sign bit, which we ignore */
01327 #ifdef Sudden_Underflow
01328        de = (int)(d0 >> Exp_shift);
01329 #ifndef IBM
01330        z |= Exp_msk11;
01331 #endif
01332 #else
01333        if (de = (int)(d0 >> Exp_shift))
01334               z |= Exp_msk1;
01335 #endif
01336 #ifdef Pack_32
01337        if (y = d1) {
01338               if (k = lo0bits(&y)) {
01339                      x[0] = y | z << 32 - k;
01340                      z >>= k;
01341                      }
01342               else
01343                      x[0] = y;
01344 #ifndef Sudden_Underflow
01345               i =
01346 #endif
01347                   b->wds = (x[1] = z) ? 2 : 1;
01348               }
01349        else {
01350 #ifdef DEBUG
01351               if (!z)
01352                      Bug("Zero passed to d2b");
01353 #endif
01354               k = lo0bits(&z);
01355               x[0] = z;
01356 #ifndef Sudden_Underflow
01357               i =
01358 #endif
01359                   b->wds = 1;
01360               k += 32;
01361               }
01362 #else
01363        if (y = d1) {
01364               if (k = lo0bits(&y))
01365                      if (k >= 16) {
01366                             x[0] = y | z << 32 - k & 0xffff;
01367                             x[1] = z >> k - 16 & 0xffff;
01368                             x[2] = z >> k;
01369                             i = 2;
01370                             }
01371                      else {
01372                             x[0] = y & 0xffff;
01373                             x[1] = y >> 16 | z << 16 - k & 0xffff;
01374                             x[2] = z >> k & 0xffff;
01375                             x[3] = z >> k+16;
01376                             i = 3;
01377                             }
01378               else {
01379                      x[0] = y & 0xffff;
01380                      x[1] = y >> 16;
01381                      x[2] = z & 0xffff;
01382                      x[3] = z >> 16;
01383                      i = 3;
01384                      }
01385               }
01386        else {
01387 #ifdef DEBUG
01388               if (!z)
01389                      Bug("Zero passed to d2b");
01390 #endif
01391               k = lo0bits(&z);
01392               if (k >= 16) {
01393                      x[0] = z;
01394                      i = 0;
01395                      }
01396               else {
01397                      x[0] = z & 0xffff;
01398                      x[1] = z >> 16;
01399                      i = 1;
01400                      }
01401               k += 32;
01402               }
01403        while(!x[i])
01404               --i;
01405        b->wds = i + 1;
01406 #endif
01407 #ifndef Sudden_Underflow
01408        if (de) {
01409 #endif
01410 #ifdef IBM
01411               *e = (de - Bias - (P-1) << 2) + k;
01412               *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
01413 #else
01414               *e = de - Bias - (P-1) + k;
01415               *bits = P - k;
01416 #endif
01417 #ifndef Sudden_Underflow
01418               }
01419        else {
01420               *e = de - Bias - (P-1) + 1 + k;
01421 #ifdef Pack_32
01422               *bits = 32*i - hi0bits(x[i-1]);
01423 #else
01424               *bits = (i+2)*16 - hi0bits(x[i]);
01425 #endif
01426               }
01427 #endif
01428        return b;
01429        }
01430 #undef d0
01431 #undef d1
01432 
01433  static double
01434 ratio
01435 #ifdef KR_headers
01436        (a, b) Bigint *a, *b;
01437 #else
01438        (Bigint *a, Bigint *b)
01439 #endif
01440 {
01441        double da, db;
01442        int k, ka, kb;
01443 
01444        dval(da) = b2d(a, &ka);
01445        dval(db) = b2d(b, &kb);
01446 #ifdef Pack_32
01447        k = ka - kb + 32*(a->wds - b->wds);
01448 #else
01449        k = ka - kb + 16*(a->wds - b->wds);
01450 #endif
01451 #ifdef IBM
01452        if (k > 0) {
01453               word0(da) += (k >> 2)*Exp_msk1;
01454               if (k &= 3)
01455                      dval(da) *= 1 << k;
01456               }
01457        else {
01458               k = -k;
01459               word0(db) += (k >> 2)*Exp_msk1;
01460               if (k &= 3)
01461                      dval(db) *= 1 << k;
01462               }
01463 #else
01464        if (k > 0)
01465               word0(da) += k*Exp_msk1;
01466        else {
01467               k = -k;
01468               word0(db) += k*Exp_msk1;
01469               }
01470 #endif
01471        return dval(da) / dval(db);
01472        }
01473 
01474  static CONST double
01475 tens[] = {
01476               1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
01477               1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
01478               1e20, 1e21, 1e22
01479 #ifdef VAX
01480               , 1e23, 1e24
01481 #endif
01482               };
01483 
01484  static CONST double
01485 #ifdef IEEE_Arith
01486 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
01487 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
01488 #ifdef Avoid_Underflow
01489               9007199254740992.*9007199254740992.e-256
01490               /* = 2^106 * 1e-53 */
01491 #else
01492               1e-256
01493 #endif
01494               };
01495 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
01496 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
01497 #define Scale_Bit 0x10
01498 #define n_bigtens 5
01499 #else
01500 #ifdef IBM
01501 bigtens[] = { 1e16, 1e32, 1e64 };
01502 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
01503 #define n_bigtens 3
01504 #else
01505 bigtens[] = { 1e16, 1e32 };
01506 static CONST double tinytens[] = { 1e-16, 1e-32 };
01507 #define n_bigtens 2
01508 #endif
01509 #endif
01510 
01511 #ifndef IEEE_Arith
01512 #undef INFNAN_CHECK
01513 #endif
01514 
01515 #ifdef INFNAN_CHECK
01516 
01517 #ifndef NAN_WORD0
01518 #define NAN_WORD0 0x7ff80000
01519 #endif
01520 
01521 #ifndef NAN_WORD1
01522 #define NAN_WORD1 0
01523 #endif
01524 
01525  static int
01526 match
01527 #ifdef KR_headers
01528        (sp, t) char **sp, *t;
01529 #else
01530        (CONST char **sp, char *t)
01531 #endif
01532 {
01533        int c, d;
01534        CONST char *s = *sp;
01535 
01536        while(d = *t++) {
01537               if ((c = *++s) >= 'A' && c <= 'Z')
01538                      c += 'a' - 'A';
01539               if (c != d)
01540                      return 0;
01541               }
01542        *sp = s + 1;
01543        return 1;
01544        }
01545 
01546 #ifndef No_Hex_NaN
01547  static void
01548 hexnan
01549 #ifdef KR_headers
01550        (rvp, sp) double *rvp; CONST char **sp;
01551 #else
01552        (double *rvp, CONST char **sp)
01553 #endif
01554 {
01555        ULong c, x[2];
01556        CONST char *s;
01557        int havedig, udx0, xshift;
01558 
01559        x[0] = x[1] = 0;
01560        havedig = xshift = 0;
01561        udx0 = 1;
01562        s = *sp;
01563        while(c = *(CONST unsigned char*)++s) {
01564               if (c >= '0' && c <= '9')
01565                      c -= '0';
01566               else if (c >= 'a' && c <= 'f')
01567                      c += 10 - 'a';
01568               else if (c >= 'A' && c <= 'F')
01569                      c += 10 - 'A';
01570               else if (c <= ' ') {
01571                      if (udx0 && havedig) {
01572                             udx0 = 0;
01573                             xshift = 1;
01574                             }
01575                      continue;
01576                      }
01577               else if (/*(*/ c == ')' && havedig) {
01578                      *sp = s + 1;
01579                      break;
01580                      }
01581               else
01582                      return;       /* invalid form: don't change *sp */
01583               havedig = 1;
01584               if (xshift) {
01585                      xshift = 0;
01586                      x[0] = x[1];
01587                      x[1] = 0;
01588                      }
01589               if (udx0)
01590                      x[0] = (x[0] << 4) | (x[1] >> 28);
01591               x[1] = (x[1] << 4) | c;
01592               }
01593        if ((x[0] &= 0xfffff) || x[1]) {
01594               word0(*rvp) = Exp_mask | x[0];
01595               word1(*rvp) = x[1];
01596               }
01597        }
01598 #endif /*No_Hex_NaN*/
01599 #endif /* INFNAN_CHECK */
01600 
01601  PR_IMPLEMENT(double)
01602 PR_strtod
01603 #ifdef KR_headers
01604        (s00, se) CONST char *s00; char **se;
01605 #else
01606        (CONST char *s00, char **se)
01607 #endif
01608 {
01609 #ifdef Avoid_Underflow
01610        int scale;
01611 #endif
01612        int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
01613                e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
01614        CONST char *s, *s0, *s1;
01615        double aadj, aadj1, adj, rv, rv0;
01616        Long L;
01617        ULong y, z;
01618        Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
01619 #ifdef SET_INEXACT
01620        int inexact, oldinexact;
01621 #endif
01622 #ifdef Honor_FLT_ROUNDS
01623        int rounding;
01624 #endif
01625 #ifdef USE_LOCALE
01626        CONST char *s2;
01627 #endif
01628 
01629        if (!_pr_initialized) _PR_ImplicitInitialization();
01630 
01631        sign = nz0 = nz = 0;
01632        dval(rv) = 0.;
01633        for(s = s00;;s++) switch(*s) {
01634               case '-':
01635                      sign = 1;
01636                      /* no break */
01637               case '+':
01638                      if (*++s)
01639                             goto break2;
01640                      /* no break */
01641               case 0:
01642                      goto ret0;
01643               case '\t':
01644               case '\n':
01645               case '\v':
01646               case '\f':
01647               case '\r':
01648               case ' ':
01649                      continue;
01650               default:
01651                      goto break2;
01652               }
01653  break2:
01654        if (*s == '0') {
01655               nz0 = 1;
01656               while(*++s == '0') ;
01657               if (!*s)
01658                      goto ret;
01659               }
01660        s0 = s;
01661        y = z = 0;
01662        for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
01663               if (nd < 9)
01664                      y = 10*y + c - '0';
01665               else if (nd < 16)
01666                      z = 10*z + c - '0';
01667        nd0 = nd;
01668 #ifdef USE_LOCALE
01669        s1 = localeconv()->decimal_point;
01670        if (c == *s1) {
01671               c = '.';
01672               if (*++s1) {
01673                      s2 = s;
01674                      for(;;) {
01675                             if (*++s2 != *s1) {
01676                                    c = 0;
01677                                    break;
01678                                    }
01679                             if (!*++s1) {
01680                                    s = s2;
01681                                    break;
01682                                    }
01683                             }
01684                      }
01685               }
01686 #endif
01687        if (c == '.') {
01688               c = *++s;
01689               if (!nd) {
01690                      for(; c == '0'; c = *++s)
01691                             nz++;
01692                      if (c > '0' && c <= '9') {
01693                             s0 = s;
01694                             nf += nz;
01695                             nz = 0;
01696                             goto have_dig;
01697                             }
01698                      goto dig_done;
01699                      }
01700               for(; c >= '0' && c <= '9'; c = *++s) {
01701  have_dig:
01702                      nz++;
01703                      if (c -= '0') {
01704                             nf += nz;
01705                             for(i = 1; i < nz; i++)
01706                                    if (nd++ < 9)
01707                                           y *= 10;
01708                                    else if (nd <= DBL_DIG + 1)
01709                                           z *= 10;
01710                             if (nd++ < 9)
01711                                    y = 10*y + c;
01712                             else if (nd <= DBL_DIG + 1)
01713                                    z = 10*z + c;
01714                             nz = 0;
01715                             }
01716                      }
01717               }
01718  dig_done:
01719        e = 0;
01720        if (c == 'e' || c == 'E') {
01721               if (!nd && !nz && !nz0) {
01722                      goto ret0;
01723                      }
01724               s00 = s;
01725               esign = 0;
01726               switch(c = *++s) {
01727                      case '-':
01728                             esign = 1;
01729                      case '+':
01730                             c = *++s;
01731                      }
01732               if (c >= '0' && c <= '9') {
01733                      while(c == '0')
01734                             c = *++s;
01735                      if (c > '0' && c <= '9') {
01736                             L = c - '0';
01737                             s1 = s;
01738                             while((c = *++s) >= '0' && c <= '9')
01739                                    L = 10*L + c - '0';
01740                             if (s - s1 > 8 || L > 19999)
01741                                    /* Avoid confusion from exponents
01742                                     * so large that e might overflow.
01743                                     */
01744                                    e = 19999; /* safe for 16 bit ints */
01745                             else
01746                                    e = (int)L;
01747                             if (esign)
01748                                    e = -e;
01749                             }
01750                      else
01751                             e = 0;
01752                      }
01753               else
01754                      s = s00;
01755               }
01756        if (!nd) {
01757               if (!nz && !nz0) {
01758 #ifdef INFNAN_CHECK
01759                      /* Check for Nan and Infinity */
01760                      switch(c) {
01761                        case 'i':
01762                        case 'I':
01763                             if (match(&s,"nf")) {
01764                                    --s;
01765                                    if (!match(&s,"inity"))
01766                                           ++s;
01767                                    word0(rv) = 0x7ff00000;
01768                                    word1(rv) = 0;
01769                                    goto ret;
01770                                    }
01771                             break;
01772                        case 'n':
01773                        case 'N':
01774                             if (match(&s, "an")) {
01775                                    word0(rv) = NAN_WORD0;
01776                                    word1(rv) = NAN_WORD1;
01777 #ifndef No_Hex_NaN
01778                                    if (*s == '(') /*)*/
01779                                           hexnan(&rv, &s);
01780 #endif
01781                                    goto ret;
01782                                    }
01783                        }
01784 #endif /* INFNAN_CHECK */
01785  ret0:
01786                      s = s00;
01787                      sign = 0;
01788                      }
01789               goto ret;
01790               }
01791        e1 = e -= nf;
01792 
01793        /* Now we have nd0 digits, starting at s0, followed by a
01794         * decimal point, followed by nd-nd0 digits.  The number we're
01795         * after is the integer represented by those digits times
01796         * 10**e */
01797 
01798        if (!nd0)
01799               nd0 = nd;
01800        k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
01801        dval(rv) = y;
01802        if (k > 9) {
01803 #ifdef SET_INEXACT
01804               if (k > DBL_DIG)
01805                      oldinexact = get_inexact();
01806 #endif
01807               dval(rv) = tens[k - 9] * dval(rv) + z;
01808               }
01809        bd0 = 0;
01810        if (nd <= DBL_DIG
01811 #ifndef RND_PRODQUOT
01812 #ifndef Honor_FLT_ROUNDS
01813               && Flt_Rounds == 1
01814 #endif
01815 #endif
01816                      ) {
01817               if (!e)
01818                      goto ret;
01819               if (e > 0) {
01820                      if (e <= Ten_pmax) {
01821 #ifdef VAX
01822                             goto vax_ovfl_check;
01823 #else
01824 #ifdef Honor_FLT_ROUNDS
01825                             /* round correctly FLT_ROUNDS = 2 or 3 */
01826                             if (sign) {
01827                                    rv = -rv;
01828                                    sign = 0;
01829                                    }
01830 #endif
01831                             /* rv = */ rounded_product(dval(rv), tens[e]);
01832                             goto ret;
01833 #endif
01834                             }
01835                      i = DBL_DIG - nd;
01836                      if (e <= Ten_pmax + i) {
01837                             /* A fancier test would sometimes let us do
01838                              * this for larger i values.
01839                              */
01840 #ifdef Honor_FLT_ROUNDS
01841                             /* round correctly FLT_ROUNDS = 2 or 3 */
01842                             if (sign) {
01843                                    rv = -rv;
01844                                    sign = 0;
01845                                    }
01846 #endif
01847                             e -= i;
01848                             dval(rv) *= tens[i];
01849 #ifdef VAX
01850                             /* VAX exponent range is so narrow we must
01851                              * worry about overflow here...
01852                              */
01853  vax_ovfl_check:
01854                             word0(rv) -= P*Exp_msk1;
01855                             /* rv = */ rounded_product(dval(rv), tens[e]);
01856                             if ((word0(rv) & Exp_mask)
01857                              > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
01858                                    goto ovfl;
01859                             word0(rv) += P*Exp_msk1;
01860 #else
01861                             /* rv = */ rounded_product(dval(rv), tens[e]);
01862 #endif
01863                             goto ret;
01864                             }
01865                      }
01866 #ifndef Inaccurate_Divide
01867               else if (e >= -Ten_pmax) {
01868 #ifdef Honor_FLT_ROUNDS
01869                      /* round correctly FLT_ROUNDS = 2 or 3 */
01870                      if (sign) {
01871                             rv = -rv;
01872                             sign = 0;
01873                             }
01874 #endif
01875                      /* rv = */ rounded_quotient(dval(rv), tens[-e]);
01876                      goto ret;
01877                      }
01878 #endif
01879               }
01880        e1 += nd - k;
01881 
01882 #ifdef IEEE_Arith
01883 #ifdef SET_INEXACT
01884        inexact = 1;
01885        if (k <= DBL_DIG)
01886               oldinexact = get_inexact();
01887 #endif
01888 #ifdef Avoid_Underflow
01889        scale = 0;
01890 #endif
01891 #ifdef Honor_FLT_ROUNDS
01892        if ((rounding = Flt_Rounds) >= 2) {
01893               if (sign)
01894                      rounding = rounding == 2 ? 0 : 2;
01895               else
01896                      if (rounding != 2)
01897                             rounding = 0;
01898               }
01899 #endif
01900 #endif /*IEEE_Arith*/
01901 
01902        /* Get starting approximation = rv * 10**e1 */
01903 
01904        if (e1 > 0) {
01905               if (i = e1 & 15)
01906                      dval(rv) *= tens[i];
01907               if (e1 &= ~15) {
01908                      if (e1 > DBL_MAX_10_EXP) {
01909  ovfl:
01910 #ifndef NO_ERRNO
01911                             PR_SetError(PR_RANGE_ERROR, 0);
01912 #endif
01913                             /* Can't trust HUGE_VAL */
01914 #ifdef IEEE_Arith
01915 #ifdef Honor_FLT_ROUNDS
01916                             switch(rounding) {
01917                               case 0: /* toward 0 */
01918                               case 3: /* toward -infinity */
01919                                    word0(rv) = Big0;
01920                                    word1(rv) = Big1;
01921                                    break;
01922                               default:
01923                                    word0(rv) = Exp_mask;
01924                                    word1(rv) = 0;
01925                               }
01926 #else /*Honor_FLT_ROUNDS*/
01927                             word0(rv) = Exp_mask;
01928                             word1(rv) = 0;
01929 #endif /*Honor_FLT_ROUNDS*/
01930 #ifdef SET_INEXACT
01931                             /* set overflow bit */
01932                             dval(rv0) = 1e300;
01933                             dval(rv0) *= dval(rv0);
01934 #endif
01935 #else /*IEEE_Arith*/
01936                             word0(rv) = Big0;
01937                             word1(rv) = Big1;
01938 #endif /*IEEE_Arith*/
01939                             if (bd0)
01940                                    goto retfree;
01941                             goto ret;
01942                             }
01943                      e1 >>= 4;
01944                      for(j = 0; e1 > 1; j++, e1 >>= 1)
01945                             if (e1 & 1)
01946                                    dval(rv) *= bigtens[j];
01947               /* The last multiplication could overflow. */
01948                      word0(rv) -= P*Exp_msk1;
01949                      dval(rv) *= bigtens[j];
01950                      if ((z = word0(rv) & Exp_mask)
01951                       > Exp_msk1*(DBL_MAX_EXP+Bias-P))
01952                             goto ovfl;
01953                      if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
01954                             /* set to largest number */
01955                             /* (Can't trust DBL_MAX) */
01956                             word0(rv) = Big0;
01957                             word1(rv) = Big1;
01958                             }
01959                      else
01960                             word0(rv) += P*Exp_msk1;
01961                      }
01962               }
01963        else if (e1 < 0) {
01964               e1 = -e1;
01965               if (i = e1 & 15)
01966                      dval(rv) /= tens[i];
01967               if (e1 >>= 4) {
01968                      if (e1 >= 1 << n_bigtens)
01969                             goto undfl;
01970 #ifdef Avoid_Underflow
01971                      if (e1 & Scale_Bit)
01972                             scale = 2*P;
01973                      for(j = 0; e1 > 0; j++, e1 >>= 1)
01974                             if (e1 & 1)
01975                                    dval(rv) *= tinytens[j];
01976                      if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
01977                                           >> Exp_shift)) > 0) {
01978                             /* scaled rv is denormal; zap j low bits */
01979                             if (j >= 32) {
01980                                    word1(rv) = 0;
01981                                    if (j >= 53)
01982                                     word0(rv) = (P+2)*Exp_msk1;
01983                                    else
01984                                     word0(rv) &= 0xffffffff << j-32;
01985                                    }
01986                             else
01987                                    word1(rv) &= 0xffffffff << j;
01988                             }
01989 #else
01990                      for(j = 0; e1 > 1; j++, e1 >>= 1)
01991                             if (e1 & 1)
01992                                    dval(rv) *= tinytens[j];
01993                      /* The last multiplication could underflow. */
01994                      dval(rv0) = dval(rv);
01995                      dval(rv) *= tinytens[j];
01996                      if (!dval(rv)) {
01997                             dval(rv) = 2.*dval(rv0);
01998                             dval(rv) *= tinytens[j];
01999 #endif
02000                             if (!dval(rv)) {
02001  undfl:
02002                                    dval(rv) = 0.;
02003 #ifndef NO_ERRNO
02004                                    PR_SetError(PR_RANGE_ERROR, 0);
02005 #endif
02006                                    if (bd0)
02007                                           goto retfree;
02008                                    goto ret;
02009                                    }
02010 #ifndef Avoid_Underflow
02011                             word0(rv) = Tiny0;
02012                             word1(rv) = Tiny1;
02013                             /* The refinement below will clean
02014                              * this approximation up.
02015                              */
02016                             }
02017 #endif
02018                      }
02019               }
02020 
02021        /* Now the hard part -- adjusting rv to the correct value.*/
02022 
02023        /* Put digits into bd: true value = bd * 10^e */
02024 
02025        bd0 = s2b(s0, nd0, nd, y);
02026 
02027        for(;;) {
02028               bd = Balloc(bd0->k);
02029               Bcopy(bd, bd0);
02030               bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
02031               bs = i2b(1);
02032 
02033               if (e >= 0) {
02034                      bb2 = bb5 = 0;
02035                      bd2 = bd5 = e;
02036                      }
02037               else {
02038                      bb2 = bb5 = -e;
02039                      bd2 = bd5 = 0;
02040                      }
02041               if (bbe >= 0)
02042                      bb2 += bbe;
02043               else
02044                      bd2 -= bbe;
02045               bs2 = bb2;
02046 #ifdef Honor_FLT_ROUNDS
02047               if (rounding != 1)
02048                      bs2++;
02049 #endif
02050 #ifdef Avoid_Underflow
02051               j = bbe - scale;
02052               i = j + bbbits - 1;  /* logb(rv) */
02053               if (i < Emin) /* denormal */
02054                      j += P - Emin;
02055               else
02056                      j = P + 1 - bbbits;
02057 #else /*Avoid_Underflow*/
02058 #ifdef Sudden_Underflow
02059 #ifdef IBM
02060               j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
02061 #else
02062               j = P + 1 - bbbits;
02063 #endif
02064 #else /*Sudden_Underflow*/
02065               j = bbe;
02066               i = j + bbbits - 1;  /* logb(rv) */
02067               if (i < Emin) /* denormal */
02068                      j += P - Emin;
02069               else
02070                      j = P + 1 - bbbits;
02071 #endif /*Sudden_Underflow*/
02072 #endif /*Avoid_Underflow*/
02073               bb2 += j;
02074               bd2 += j;
02075 #ifdef Avoid_Underflow
02076               bd2 += scale;
02077 #endif
02078               i = bb2 < bd2 ? bb2 : bd2;
02079               if (i > bs2)
02080                      i = bs2;
02081               if (i > 0) {
02082                      bb2 -= i;
02083                      bd2 -= i;
02084                      bs2 -= i;
02085                      }
02086               if (bb5 > 0) {
02087                      bs = pow5mult(bs, bb5);
02088                      bb1 = mult(bs, bb);
02089                      Bfree(bb);
02090                      bb = bb1;
02091                      }
02092               if (bb2 > 0)
02093                      bb = lshift(bb, bb2);
02094               if (bd5 > 0)
02095                      bd = pow5mult(bd, bd5);
02096               if (bd2 > 0)
02097                      bd = lshift(bd, bd2);
02098               if (bs2 > 0)
02099                      bs = lshift(bs, bs2);
02100               delta = diff(bb, bd);
02101               dsign = delta->sign;
02102               delta->sign = 0;
02103               i = cmp(delta, bs);
02104 #ifdef Honor_FLT_ROUNDS
02105               if (rounding != 1) {
02106                      if (i < 0) {
02107                             /* Error is less than an ulp */
02108                             if (!delta->x[0] && delta->wds <= 1) {
02109                                    /* exact */
02110 #ifdef SET_INEXACT
02111                                    inexact = 0;
02112 #endif
02113                                    break;
02114                                    }
02115                             if (rounding) {
02116                                    if (dsign) {
02117                                           adj = 1.;
02118                                           goto apply_adj;
02119                                           }
02120                                    }
02121                             else if (!dsign) {
02122                                    adj = -1.;
02123                                    if (!word1(rv)
02124                                     && !(word0(rv) & Frac_mask)) {
02125                                           y = word0(rv) & Exp_mask;
02126 #ifdef Avoid_Underflow
02127                                           if (!scale || y > 2*P*Exp_msk1)
02128 #else
02129                                           if (y)
02130 #endif
02131                                             {
02132                                             delta = lshift(delta,Log2P);
02133                                             if (cmp(delta, bs) <= 0)
02134                                                  adj = -0.5;
02135                                             }
02136                                           }
02137  apply_adj:
02138 #ifdef Avoid_Underflow
02139                                    if (scale && (y = word0(rv) & Exp_mask)
02140                                           <= 2*P*Exp_msk1)
02141                                      word0(adj) += (2*P+1)*Exp_msk1 - y;
02142 #else
02143 #ifdef Sudden_Underflow
02144                                    if ((word0(rv) & Exp_mask) <=
02145                                                  P*Exp_msk1) {
02146                                           word0(rv) += P*Exp_msk1;
02147                                           dval(rv) += adj*ulp(dval(rv));
02148                                           word0(rv) -= P*Exp_msk1;
02149                                           }
02150                                    else
02151 #endif /*Sudden_Underflow*/
02152 #endif /*Avoid_Underflow*/
02153                                    dval(rv) += adj*ulp(dval(rv));
02154                                    }
02155                             break;
02156                             }
02157                      adj = ratio(delta, bs);
02158                      if (adj < 1.)
02159                             adj = 1.;
02160                      if (adj <= 0x7ffffffe) {
02161                             /* adj = rounding ? ceil(adj) : floor(adj); */
02162                             y = adj;
02163                             if (y != adj) {
02164                                    if (!((rounding>>1) ^ dsign))
02165                                           y++;
02166                                    adj = y;
02167                                    }
02168                             }
02169 #ifdef Avoid_Underflow
02170                      if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02171                             word0(adj) += (2*P+1)*Exp_msk1 - y;
02172 #else
02173 #ifdef Sudden_Underflow
02174                      if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02175                             word0(rv) += P*Exp_msk1;
02176                             adj *= ulp(dval(rv));
02177                             if (dsign)
02178                                    dval(rv) += adj;
02179                             else
02180                                    dval(rv) -= adj;
02181                             word0(rv) -= P*Exp_msk1;
02182                             goto cont;
02183                             }
02184 #endif /*Sudden_Underflow*/
02185 #endif /*Avoid_Underflow*/
02186                      adj *= ulp(dval(rv));
02187                      if (dsign)
02188                             dval(rv) += adj;
02189                      else
02190                             dval(rv) -= adj;
02191                      goto cont;
02192                      }
02193 #endif /*Honor_FLT_ROUNDS*/
02194 
02195               if (i < 0) {
02196                      /* Error is less than half an ulp -- check for
02197                       * special case of mantissa a power of two.
02198                       */
02199                      if (dsign || word1(rv) || word0(rv) & Bndry_mask
02200 #ifdef IEEE_Arith
02201 #ifdef Avoid_Underflow
02202                       || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
02203 #else
02204                       || (word0(rv) & Exp_mask) <= Exp_msk1
02205 #endif
02206 #endif
02207                             ) {
02208 #ifdef SET_INEXACT
02209                             if (!delta->x[0] && delta->wds <= 1)
02210                                    inexact = 0;
02211 #endif
02212                             break;
02213                             }
02214                      if (!delta->x[0] && delta->wds <= 1) {
02215                             /* exact result */
02216 #ifdef SET_INEXACT
02217                             inexact = 0;
02218 #endif
02219                             break;
02220                             }
02221                      delta = lshift(delta,Log2P);
02222                      if (cmp(delta, bs) > 0)
02223                             goto drop_down;
02224                      break;
02225                      }
02226               if (i == 0) {
02227                      /* exactly half-way between */
02228                      if (dsign) {
02229                             if ((word0(rv) & Bndry_mask1) == Bndry_mask1
02230                              &&  word1(rv) == (
02231 #ifdef Avoid_Underflow
02232                      (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
02233               ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
02234 #endif
02235                                              0xffffffff)) {
02236                                    /*boundary case -- increment exponent*/
02237                                    word0(rv) = (word0(rv) & Exp_mask)
02238                                           + Exp_msk1
02239 #ifdef IBM
02240                                           | Exp_msk1 >> 4
02241 #endif
02242                                           ;
02243                                    word1(rv) = 0;
02244 #ifdef Avoid_Underflow
02245                                    dsign = 0;
02246 #endif
02247                                    break;
02248                                    }
02249                             }
02250                      else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
02251  drop_down:
02252                             /* boundary case -- decrement exponent */
02253 #ifdef Sudden_Underflow /*{{*/
02254                             L = word0(rv) & Exp_mask;
02255 #ifdef IBM
02256                             if (L <  Exp_msk1)
02257 #else
02258 #ifdef Avoid_Underflow
02259                             if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
02260 #else
02261                             if (L <= Exp_msk1)
02262 #endif /*Avoid_Underflow*/
02263 #endif /*IBM*/
02264                                    goto undfl;
02265                             L -= Exp_msk1;
02266 #else /*Sudden_Underflow}{*/
02267 #ifdef Avoid_Underflow
02268                             if (scale) {
02269                                    L = word0(rv) & Exp_mask;
02270                                    if (L <= (2*P+1)*Exp_msk1) {
02271                                           if (L > (P+2)*Exp_msk1)
02272                                                  /* round even ==> */
02273                                                  /* accept rv */
02274                                                  break;
02275                                           /* rv = smallest denormal */
02276                                           goto undfl;
02277                                           }
02278                                    }
02279 #endif /*Avoid_Underflow*/
02280                             L = (word0(rv) & Exp_mask) - Exp_msk1;
02281 #endif /*Sudden_Underflow}}*/
02282                             word0(rv) = L | Bndry_mask1;
02283                             word1(rv) = 0xffffffff;
02284 #ifdef IBM
02285                             goto cont;
02286 #else
02287                             break;
02288 #endif
02289                             }
02290 #ifndef ROUND_BIASED
02291                      if (!(word1(rv) & LSB))
02292                             break;
02293 #endif
02294                      if (dsign)
02295                             dval(rv) += ulp(dval(rv));
02296 #ifndef ROUND_BIASED
02297                      else {
02298                             dval(rv) -= ulp(dval(rv));
02299 #ifndef Sudden_Underflow
02300                             if (!dval(rv))
02301                                    goto undfl;
02302 #endif
02303                             }
02304 #ifdef Avoid_Underflow
02305                      dsign = 1 - dsign;
02306 #endif
02307 #endif
02308                      break;
02309                      }
02310               if ((aadj = ratio(delta, bs)) <= 2.) {
02311                      if (dsign)
02312                             aadj = aadj1 = 1.;
02313                      else if (word1(rv) || word0(rv) & Bndry_mask) {
02314 #ifndef Sudden_Underflow
02315                             if (word1(rv) == Tiny1 && !word0(rv))
02316                                    goto undfl;
02317 #endif
02318                             aadj = 1.;
02319                             aadj1 = -1.;
02320                             }
02321                      else {
02322                             /* special case -- power of FLT_RADIX to be */
02323                             /* rounded down... */
02324 
02325                             if (aadj < 2./FLT_RADIX)
02326                                    aadj = 1./FLT_RADIX;
02327                             else
02328                                    aadj *= 0.5;
02329                             aadj1 = -aadj;
02330                             }
02331                      }
02332               else {
02333                      aadj *= 0.5;
02334                      aadj1 = dsign ? aadj : -aadj;
02335 #ifdef Check_FLT_ROUNDS
02336                      switch(Rounding) {
02337                             case 2: /* towards +infinity */
02338                                    aadj1 -= 0.5;
02339                                    break;
02340                             case 0: /* towards 0 */
02341                             case 3: /* towards -infinity */
02342                                    aadj1 += 0.5;
02343                             }
02344 #else
02345                      if (Flt_Rounds == 0)
02346                             aadj1 += 0.5;
02347 #endif /*Check_FLT_ROUNDS*/
02348                      }
02349               y = word0(rv) & Exp_mask;
02350 
02351               /* Check for overflow */
02352 
02353               if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
02354                      dval(rv0) = dval(rv);
02355                      word0(rv) -= P*Exp_msk1;
02356                      adj = aadj1 * ulp(dval(rv));
02357                      dval(rv) += adj;
02358                      if ((word0(rv) & Exp_mask) >=
02359                                    Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
02360                             if (word0(rv0) == Big0 && word1(rv0) == Big1)
02361                                    goto ovfl;
02362                             word0(rv) = Big0;
02363                             word1(rv) = Big1;
02364                             goto cont;
02365                             }
02366                      else
02367                             word0(rv) += P*Exp_msk1;
02368                      }
02369               else {
02370 #ifdef Avoid_Underflow
02371                      if (scale && y <= 2*P*Exp_msk1) {
02372                             if (aadj <= 0x7fffffff) {
02373                                    if ((z = aadj) <= 0)
02374                                           z = 1;
02375                                    aadj = z;
02376                                    aadj1 = dsign ? aadj : -aadj;
02377                                    }
02378                             word0(aadj1) += (2*P+1)*Exp_msk1 - y;
02379                             }
02380                      adj = aadj1 * ulp(dval(rv));
02381                      dval(rv) += adj;
02382 #else
02383 #ifdef Sudden_Underflow
02384                      if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
02385                             dval(rv0) = dval(rv);
02386                             word0(rv) += P*Exp_msk1;
02387                             adj = aadj1 * ulp(dval(rv));
02388                             dval(rv) += adj;
02389 #ifdef IBM
02390                             if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
02391 #else
02392                             if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
02393 #endif
02394                                    {
02395                                    if (word0(rv0) == Tiny0
02396                                     && word1(rv0) == Tiny1)
02397                                           goto undfl;
02398                                    word0(rv) = Tiny0;
02399                                    word1(rv) = Tiny1;
02400                                    goto cont;
02401                                    }
02402                             else
02403                                    word0(rv) -= P*Exp_msk1;
02404                             }
02405                      else {
02406                             adj = aadj1 * ulp(dval(rv));
02407                             dval(rv) += adj;
02408                             }
02409 #else /*Sudden_Underflow*/
02410                      /* Compute adj so that the IEEE rounding rules will
02411                       * correctly round rv + adj in some half-way cases.
02412                       * If rv * ulp(rv) is denormalized (i.e.,
02413                       * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
02414                       * trouble from bits lost to denormalization;
02415                       * example: 1.2e-307 .
02416                       */
02417                      if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
02418                             aadj1 = (double)(int)(aadj + 0.5);
02419                             if (!dsign)
02420                                    aadj1 = -aadj1;
02421                             }
02422                      adj = aadj1 * ulp(dval(rv));
02423                      dval(rv) += adj;
02424 #endif /*Sudden_Underflow*/
02425 #endif /*Avoid_Underflow*/
02426                      }
02427               z = word0(rv) & Exp_mask;
02428 #ifndef SET_INEXACT
02429 #ifdef Avoid_Underflow
02430               if (!scale)
02431 #endif
02432               if (y == z) {
02433                      /* Can we stop now? */
02434                      L = (Long)aadj;
02435                      aadj -= L;
02436                      /* The tolerances below are conservative. */
02437                      if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
02438                             if (aadj < .4999999 || aadj > .5000001)
02439                                    break;
02440                             }
02441                      else if (aadj < .4999999/FLT_RADIX)
02442                             break;
02443                      }
02444 #endif
02445  cont:
02446               Bfree(bb);
02447               Bfree(bd);
02448               Bfree(bs);
02449               Bfree(delta);
02450               }
02451 #ifdef SET_INEXACT
02452        if (inexact) {
02453               if (!oldinexact) {
02454                      word0(rv0) = Exp_1 + (70 << Exp_shift);
02455                      word1(rv0) = 0;
02456                      dval(rv0) += 1.;
02457                      }
02458               }
02459        else if (!oldinexact)
02460               clear_inexact();
02461 #endif
02462 #ifdef Avoid_Underflow
02463        if (scale) {
02464               word0(rv0) = Exp_1 - 2*P*Exp_msk1;
02465               word1(rv0) = 0;
02466               dval(rv) *= dval(rv0);
02467 #ifndef NO_ERRNO
02468               /* try to avoid the bug of testing an 8087 register value */
02469               if (word0(rv) == 0 && word1(rv) == 0)
02470                      PR_SetError(PR_RANGE_ERROR, 0);
02471 #endif
02472               }
02473 #endif /* Avoid_Underflow */
02474 #ifdef SET_INEXACT
02475        if (inexact && !(word0(rv) & Exp_mask)) {
02476               /* set underflow bit */
02477               dval(rv0) = 1e-300;
02478               dval(rv0) *= dval(rv0);
02479               }
02480 #endif
02481  retfree:
02482        Bfree(bb);
02483        Bfree(bd);
02484        Bfree(bs);
02485        Bfree(bd0);
02486        Bfree(delta);
02487  ret:
02488        if (se)
02489               *se = (char *)s;
02490        return sign ? -dval(rv) : dval(rv);
02491        }
02492 
02493  static int
02494 quorem
02495 #ifdef KR_headers
02496        (b, S) Bigint *b, *S;
02497 #else
02498        (Bigint *b, Bigint *S)
02499 #endif
02500 {
02501        int n;
02502        ULong *bx, *bxe, q, *sx, *sxe;
02503 #ifdef ULLong
02504        ULLong borrow, carry, y, ys;
02505 #else
02506        ULong borrow, carry, y, ys;
02507 #ifdef Pack_32
02508        ULong si, z, zs;
02509 #endif
02510 #endif
02511 
02512        n = S->wds;
02513 #ifdef DEBUG
02514        /*debug*/ if (b->wds > n)
02515        /*debug*/     Bug("oversize b in quorem");
02516 #endif
02517        if (b->wds < n)
02518               return 0;
02519        sx = S->x;
02520        sxe = sx + --n;
02521        bx = b->x;
02522        bxe = bx + n;
02523        q = *bxe / (*sxe + 1);      /* ensure q <= true quotient */
02524 #ifdef DEBUG
02525        /*debug*/ if (q > 9)
02526        /*debug*/     Bug("oversized quotient in quorem");
02527 #endif
02528        if (q) {
02529               borrow = 0;
02530               carry = 0;
02531               do {
02532 #ifdef ULLong
02533                      ys = *sx++ * (ULLong)q + carry;
02534                      carry = ys >> 32;
02535                      y = *bx - (ys & FFFFFFFF) - borrow;
02536                      borrow = y >> 32 & (ULong)1;
02537                      *bx++ = y & FFFFFFFF;
02538 #else
02539 #ifdef Pack_32
02540                      si = *sx++;
02541                      ys = (si & 0xffff) * q + carry;
02542                      zs = (si >> 16) * q + (ys >> 16);
02543                      carry = zs >> 16;
02544                      y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02545                      borrow = (y & 0x10000) >> 16;
02546                      z = (*bx >> 16) - (zs & 0xffff) - borrow;
02547                      borrow = (z & 0x10000) >> 16;
02548                      Storeinc(bx, z, y);
02549 #else
02550                      ys = *sx++ * q + carry;
02551                      carry = ys >> 16;
02552                      y = *bx - (ys & 0xffff) - borrow;
02553                      borrow = (y & 0x10000) >> 16;
02554                      *bx++ = y & 0xffff;
02555 #endif
02556 #endif
02557                      }
02558                      while(sx <= sxe);
02559               if (!*bxe) {
02560                      bx = b->x;
02561                      while(--bxe > bx && !*bxe)
02562                             --n;
02563                      b->wds = n;
02564                      }
02565               }
02566        if (cmp(b, S) >= 0) {
02567               q++;
02568               borrow = 0;
02569               carry = 0;
02570               bx = b->x;
02571               sx = S->x;
02572               do {
02573 #ifdef ULLong
02574                      ys = *sx++ + carry;
02575                      carry = ys >> 32;
02576                      y = *bx - (ys & FFFFFFFF) - borrow;
02577                      borrow = y >> 32 & (ULong)1;
02578                      *bx++ = y & FFFFFFFF;
02579 #else
02580 #ifdef Pack_32
02581                      si = *sx++;
02582                      ys = (si & 0xffff) + carry;
02583                      zs = (si >> 16) + (ys >> 16);
02584                      carry = zs >> 16;
02585                      y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
02586                      borrow = (y & 0x10000) >> 16;
02587                      z = (*bx >> 16) - (zs & 0xffff) - borrow;
02588                      borrow = (z & 0x10000) >> 16;
02589                      Storeinc(bx, z, y);
02590 #else
02591                      ys = *sx++ + carry;
02592                      carry = ys >> 16;
02593                      y = *bx - (ys & 0xffff) - borrow;
02594                      borrow = (y & 0x10000) >> 16;
02595                      *bx++ = y & 0xffff;
02596 #endif
02597 #endif
02598                      }
02599                      while(sx <= sxe);
02600               bx = b->x;
02601               bxe = bx + n;
02602               if (!*bxe) {
02603                      while(--bxe > bx && !*bxe)
02604                             --n;
02605                      b->wds = n;
02606                      }
02607               }
02608        return q;
02609        }
02610 
02611 #ifndef MULTIPLE_THREADS
02612  static char *dtoa_result;
02613 #endif
02614 
02615  static char *
02616 #ifdef KR_headers
02617 rv_alloc(i) int i;
02618 #else
02619 rv_alloc(int i)
02620 #endif
02621 {
02622        int j, k, *r;
02623 
02624        j = sizeof(ULong);
02625        for(k = 0;
02626               sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
02627               j <<= 1)
02628                      k++;
02629        r = (int*)Balloc(k);
02630        *r = k;
02631        return
02632 #ifndef MULTIPLE_THREADS
02633        dtoa_result =
02634 #endif
02635               (char *)(r+1);
02636        }
02637 
02638  static char *
02639 #ifdef KR_headers
02640 nrv_alloc(s, rve, n) char *s, **rve; int n;
02641 #else
02642 nrv_alloc(char *s, char **rve, int n)
02643 #endif
02644 {
02645        char *rv, *t;
02646 
02647        t = rv = rv_alloc(n);
02648        while(*t = *s++) t++;
02649        if (rve)
02650               *rve = t;
02651        return rv;
02652        }
02653 
02654 /* freedtoa(s) must be used to free values s returned by dtoa
02655  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
02656  * but for consistency with earlier versions of dtoa, it is optional
02657  * when MULTIPLE_THREADS is not defined.
02658  */
02659 
02660  void
02661 #ifdef KR_headers
02662 freedtoa(s) char *s;
02663 #else
02664 freedtoa(char *s)
02665 #endif
02666 {
02667        Bigint *b = (Bigint *)((int *)s - 1);
02668        b->maxwds = 1 << (b->k = *(int*)b);
02669        Bfree(b);
02670 #ifndef MULTIPLE_THREADS
02671        if (s == dtoa_result)
02672               dtoa_result = 0;
02673 #endif
02674        }
02675 
02676 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
02677  *
02678  * Inspired by "How to Print Floating-Point Numbers Accurately" by
02679  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
02680  *
02681  * Modifications:
02682  *     1. Rather than iterating, we use a simple numeric overestimate
02683  *        to determine k = floor(log10(d)).  We scale relevant
02684  *        quantities using O(log2(k)) rather than O(k) multiplications.
02685  *     2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
02686  *        try to generate digits strictly left to right.  Instead, we
02687  *        compute with fewer bits and propagate the carry if necessary
02688  *        when rounding the final digit up.  This is often faster.
02689  *     3. Under the assumption that input will be rounded nearest,
02690  *        mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
02691  *        That is, we allow equality in stopping tests when the
02692  *        round-nearest rule will give the same floating-point value
02693  *        as would satisfaction of the stopping test with strict
02694  *        inequality.
02695  *     4. We remove common factors of powers of 2 from relevant
02696  *        quantities.
02697  *     5. When converting floating-point integers less than 1e16,
02698  *        we use floating-point arithmetic rather than resorting
02699  *        to multiple-precision integers.
02700  *     6. When asked to produce fewer than 15 digits, we first try
02701  *        to get by with floating-point arithmetic; we resort to
02702  *        multiple-precision integer arithmetic only if we cannot
02703  *        guarantee that the floating-point calculation has given
02704  *        the correctly rounded result.  For k requested digits and
02705  *        "uniformly" distributed input, the probability is
02706  *        something like 10^(k-15) that we must resort to the Long
02707  *        calculation.
02708  */
02709 
02710  static char *
02711 dtoa
02712 #ifdef KR_headers
02713        (d, mode, ndigits, decpt, sign, rve)
02714        double d; int mode, ndigits, *decpt, *sign; char **rve;
02715 #else
02716        (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
02717 #endif
02718 {
02719  /*    Arguments ndigits, decpt, sign are similar to those
02720        of ecvt and fcvt; trailing zeros are suppressed from
02721        the returned string.  If not null, *rve is set to point
02722        to the end of the return value.  If d is +-Infinity or NaN,
02723        then *decpt is set to 9999.
02724 
02725        mode:
02726               0 ==> shortest string that yields d when read in
02727                      and rounded to nearest.
02728               1 ==> like 0, but with Steele & White stopping rule;
02729                      e.g. with IEEE P754 arithmetic , mode 0 gives
02730                      1e23 whereas mode 1 gives 9.999999999999999e22.
02731               2 ==> max(1,ndigits) significant digits.  This gives a
02732                      return value similar to that of ecvt, except
02733                      that trailing zeros are suppressed.
02734               3 ==> through ndigits past the decimal point.  This
02735                      gives a return value similar to that from fcvt,
02736                      except that trailing zeros are suppressed, and
02737                      ndigits can be negative.
02738               4,5 ==> similar to 2 and 3, respectively, but (in
02739                      round-nearest mode) with the tests of mode 0 to
02740                      possibly return a shorter string that rounds to d.
02741                      With IEEE arithmetic and compilation with
02742                      -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
02743                      as modes 2 and 3 when FLT_ROUNDS != 1.
02744               6-9 ==> Debugging modes similar to mode - 4:  don't try
02745                      fast floating-point estimate (if applicable).
02746 
02747               Values of mode other than 0-9 are treated as mode 0.
02748 
02749               Sufficient space is allocated to the return value
02750               to hold the suppressed trailing zeros.
02751        */
02752 
02753        int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
02754               j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
02755               spec_case, try_quick;
02756        Long L;
02757 #ifndef Sudden_Underflow
02758        int denorm;
02759        ULong x;
02760 #endif
02761        Bigint *b, *b1, *delta, *mlo, *mhi, *S;
02762        double d2, ds, eps;
02763        char *s, *s0;
02764 #ifdef Honor_FLT_ROUNDS
02765        int rounding;
02766 #endif
02767 #ifdef SET_INEXACT
02768        int inexact, oldinexact;
02769 #endif
02770 
02771 #ifndef MULTIPLE_THREADS
02772        if (dtoa_result) {
02773               freedtoa(dtoa_result);
02774               dtoa_result = 0;
02775               }
02776 #endif
02777 
02778        if (word0(d) & Sign_bit) {
02779               /* set sign for everything, including 0's and NaNs */
02780               *sign = 1;
02781               word0(d) &= ~Sign_bit;      /* clear sign bit */
02782               }
02783        else
02784               *sign = 0;
02785 
02786 #if defined(IEEE_Arith) + defined(VAX)
02787 #ifdef IEEE_Arith
02788        if ((word0(d) & Exp_mask) == Exp_mask)
02789 #else
02790        if (word0(d)  == 0x8000)
02791 #endif
02792               {
02793               /* Infinity or NaN */
02794               *decpt = 9999;
02795 #ifdef IEEE_Arith
02796               if (!word1(d) && !(word0(d) & 0xfffff))
02797                      return nrv_alloc("Infinity", rve, 8);
02798 #endif
02799               return nrv_alloc("NaN", rve, 3);
02800               }
02801 #endif
02802 #ifdef IBM
02803        dval(d) += 0; /* normalize */
02804 #endif
02805        if (!dval(d)) {
02806               *decpt = 1;
02807               return nrv_alloc("0", rve, 1);
02808               }
02809 
02810 #ifdef SET_INEXACT
02811        try_quick = oldinexact = get_inexact();
02812        inexact = 1;
02813 #endif
02814 #ifdef Honor_FLT_ROUNDS
02815        if ((rounding = Flt_Rounds) >= 2) {
02816               if (*sign)
02817                      rounding = rounding == 2 ? 0 : 2;
02818               else
02819                      if (rounding != 2)
02820                             rounding = 0;
02821               }
02822 #endif
02823 
02824        b = d2b(dval(d), &be, &bbits);
02825 #ifdef Sudden_Underflow
02826        i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
02827 #else
02828        if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
02829 #endif
02830               dval(d2) = dval(d);
02831               word0(d2) &= Frac_mask1;
02832               word0(d2) |= Exp_11;
02833 #ifdef IBM
02834               if (j = 11 - hi0bits(word0(d2) & Frac_mask))
02835                      dval(d2) /= 1 << j;
02836 #endif
02837 
02838               /* log(x)     ~=~ log(1.5) + (x-1.5)/1.5
02839                * log10(x)    =  log(x) / log(10)
02840                *            ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
02841                * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
02842                *
02843                * This suggests computing an approximation k to log10(d) by
02844                *
02845                * k = (i - Bias)*0.301029995663981
02846                *     + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
02847                *
02848                * We want k to be too large rather than too small.
02849                * The error in the first-order Taylor series approximation
02850                * is in our favor, so we just round up the constant enough
02851                * to compensate for any error in the multiplication of
02852                * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
02853                * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
02854                * adding 1e-13 to the constant term more than suffices.
02855                * Hence we adjust the constant term to 0.1760912590558.
02856                * (We could get a more accurate k by invoking log10,
02857                *  but this is probably not worthwhile.)
02858                */
02859 
02860               i -= Bias;
02861 #ifdef IBM
02862               i <<= 2;
02863               i += j;
02864 #endif
02865 #ifndef Sudden_Underflow
02866               denorm = 0;
02867               }
02868        else {
02869               /* d is denormalized */
02870 
02871               i = bbits + be + (Bias + (P-1) - 1);
02872               x = i > 32  ? word0(d) << 64 - i | word1(d) >> i - 32
02873                          : word1(d) << 32 - i;
02874               dval(d2) = x;
02875               word0(d2) -= 31*Exp_msk1; /* adjust exponent */
02876               i -= (Bias + (P-1) - 1) + 1;
02877               denorm = 1;
02878               }
02879 #endif
02880        ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
02881        k = (int)ds;
02882        if (ds < 0. && ds != k)
02883               k--;   /* want k = floor(ds) */
02884        k_check = 1;
02885        if (k >= 0 && k <= Ten_pmax) {
02886               if (dval(d) < tens[k])
02887                      k--;
02888               k_check = 0;
02889               }
02890        j = bbits - i - 1;
02891        if (j >= 0) {
02892               b2 = 0;
02893               s2 = j;
02894               }
02895        else {
02896               b2 = -j;
02897               s2 = 0;
02898               }
02899        if (k >= 0) {
02900               b5 = 0;
02901               s5 = k;
02902               s2 += k;
02903               }
02904        else {
02905               b2 -= k;
02906               b5 = -k;
02907               s5 = 0;
02908               }
02909        if (mode < 0 || mode > 9)
02910               mode = 0;
02911 
02912 #ifndef SET_INEXACT
02913 #ifdef Check_FLT_ROUNDS
02914        try_quick = Rounding == 1;
02915 #else
02916        try_quick = 1;
02917 #endif
02918 #endif /*SET_INEXACT*/
02919 
02920        if (mode > 5) {
02921               mode -= 4;
02922               try_quick = 0;
02923               }
02924        leftright = 1;
02925        switch(mode) {
02926               case 0:
02927               case 1:
02928                      ilim = ilim1 = -1;
02929                      i = 18;
02930                      ndigits = 0;
02931                      break;
02932               case 2:
02933                      leftright = 0;
02934                      /* no break */
02935               case 4:
02936                      if (ndigits <= 0)
02937                             ndigits = 1;
02938                      ilim = ilim1 = i = ndigits;
02939                      break;
02940               case 3:
02941                      leftright = 0;
02942                      /* no break */
02943               case 5:
02944                      i = ndigits + k + 1;
02945                      ilim = i;
02946                      ilim1 = i - 1;
02947                      if (i <= 0)
02948                             i = 1;
02949               }
02950        s = s0 = rv_alloc(i);
02951 
02952 #ifdef Honor_FLT_ROUNDS
02953        if (mode > 1 && rounding != 1)
02954               leftright = 0;
02955 #endif
02956 
02957        if (ilim >= 0 && ilim <= Quick_max && try_quick) {
02958 
02959               /* Try to get by with floating-point arithmetic. */
02960 
02961               i = 0;
02962               dval(d2) = dval(d);
02963               k0 = k;
02964               ilim0 = ilim;
02965               ieps = 2; /* conservative */
02966               if (k > 0) {
02967                      ds = tens[k&0xf];
02968                      j = k >> 4;
02969                      if (j & Bletch) {
02970                             /* prevent overflows */
02971                             j &= Bletch - 1;
02972                             dval(d) /= bigtens[n_bigtens-1];
02973                             ieps++;
02974                             }
02975                      for(; j; j >>= 1, i++)
02976                             if (j & 1) {
02977                                    ieps++;
02978                                    ds *= bigtens[i];
02979                                    }
02980                      dval(d) /= ds;
02981                      }
02982               else if (j1 = -k) {
02983                      dval(d) *= tens[j1 & 0xf];
02984                      for(j = j1 >> 4; j; j >>= 1, i++)
02985                             if (j & 1) {
02986                                    ieps++;
02987                                    dval(d) *= bigtens[i];
02988                                    }
02989                      }
02990               if (k_check && dval(d) < 1. && ilim > 0) {
02991                      if (ilim1 <= 0)
02992                             goto fast_failed;
02993                      ilim = ilim1;
02994                      k--;
02995                      dval(d) *= 10.;
02996                      ieps++;
02997                      }
02998               dval(eps) = ieps*dval(d) + 7.;
02999               word0(eps) -= (P-1)*Exp_msk1;
03000               if (ilim == 0) {
03001                      S = mhi = 0;
03002                      dval(d) -= 5.;
03003                      if (dval(d) > dval(eps))
03004                             goto one_digit;
03005                      if (dval(d) < -dval(eps))
03006                             goto no_digits;
03007                      goto fast_failed;
03008                      }
03009 #ifndef No_leftright
03010               if (leftright) {
03011                      /* Use Steele & White method of only
03012                       * generating digits needed.
03013                       */
03014                      dval(eps) = 0.5/tens[ilim-1] - dval(eps);
03015                      for(i = 0;;) {
03016                             L = dval(d);
03017                             dval(d) -= L;
03018                             *s++ = '0' + (int)L;
03019                             if (dval(d) < dval(eps))
03020                                    goto ret1;
03021                             if (1. - dval(d) < dval(eps))
03022                                    goto bump_up;
03023                             if (++i >= ilim)
03024                                    break;
03025                             dval(eps) *= 10.;
03026                             dval(d) *= 10.;
03027                             }
03028                      }
03029               else {
03030 #endif
03031                      /* Generate ilim digits, then fix them up. */
03032                      dval(eps) *= tens[ilim-1];
03033                      for(i = 1;; i++, dval(d) *= 10.) {
03034                             L = (Long)(dval(d));
03035                             if (!(dval(d) -= L))
03036                                    ilim = i;
03037                             *s++ = '0' + (int)L;
03038                             if (i == ilim) {
03039                                    if (dval(d) > 0.5 + dval(eps))
03040                                           goto bump_up;
03041                                    else if (dval(d) < 0.5 - dval(eps)) {
03042                                           while(*--s == '0');
03043                                           s++;
03044                                           goto ret1;
03045                                           }
03046                                    break;
03047                                    }
03048                             }
03049 #ifndef No_leftright
03050                      }
03051 #endif
03052  fast_failed:
03053               s = s0;
03054               dval(d) = dval(d2);
03055               k = k0;
03056               ilim = ilim0;
03057               }
03058 
03059        /* Do we have a "small" integer? */
03060 
03061        if (be >= 0 && k <= Int_max) {
03062               /* Yes. */
03063               ds = tens[k];
03064               if (ndigits < 0 && ilim <= 0) {
03065                      S = mhi = 0;
03066                      if (ilim < 0 || dval(d) <= 5*ds)
03067                             goto no_digits;
03068                      goto one_digit;
03069                      }
03070               for(i = 1; i <= k+1; i++, dval(d) *= 10.) {
03071                      L = (Long)(dval(d) / ds);
03072                      dval(d) -= L*ds;
03073 #ifdef Check_FLT_ROUNDS
03074                      /* If FLT_ROUNDS == 2, L will usually be high by 1 */
03075                      if (dval(d) < 0) {
03076                             L--;
03077                             dval(d) += ds;
03078                             }
03079 #endif
03080                      *s++ = '0' + (int)L;
03081                      if (!dval(d)) {
03082 #ifdef SET_INEXACT
03083                             inexact = 0;
03084 #endif
03085                             break;
03086                             }
03087                      if (i == ilim) {
03088 #ifdef Honor_FLT_ROUNDS
03089                             if (mode > 1)
03090                             switch(rounding) {
03091                               case 0: goto ret1;
03092                               case 2: goto bump_up;
03093                               }
03094 #endif
03095                             dval(d) += dval(d);
03096                             if (dval(d) > ds || dval(d) == ds && L & 1) {
03097  bump_up:
03098                                    while(*--s == '9')
03099                                           if (s == s0) {
03100                                                  k++;
03101                                                  *s = '0';
03102                                                  break;
03103                                                  }
03104                                    ++*s++;
03105                                    }
03106                             break;
03107                             }
03108                      }
03109               goto ret1;
03110               }
03111 
03112        m2 = b2;
03113        m5 = b5;
03114        mhi = mlo = 0;
03115        if (leftright) {
03116               i =
03117 #ifndef Sudden_Underflow
03118                      denorm ? be + (Bias + (P-1) - 1 + 1) :
03119 #endif
03120 #ifdef IBM
03121                      1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
03122 #else
03123                      1 + P - bbits;
03124 #endif
03125               b2 += i;
03126               s2 += i;
03127               mhi = i2b(1);
03128               }
03129        if (m2 > 0 && s2 > 0) {
03130               i = m2 < s2 ? m2 : s2;
03131               b2 -= i;
03132               m2 -= i;
03133               s2 -= i;
03134               }
03135        if (b5 > 0) {
03136               if (leftright) {
03137                      if (m5 > 0) {
03138                             mhi = pow5mult(mhi, m5);
03139                             b1 = mult(mhi, b);
03140                             Bfree(b);
03141                             b = b1;
03142                             }
03143                      if (j = b5 - m5)
03144                             b = pow5mult(b, j);
03145                      }
03146               else
03147                      b = pow5mult(b, b5);
03148               }
03149        S = i2b(1);
03150        if (s5 > 0)
03151               S = pow5mult(S, s5);
03152 
03153        /* Check for special case that d is a normalized power of 2. */
03154 
03155        spec_case = 0;
03156        if ((mode < 2 || leftright)
03157 #ifdef Honor_FLT_ROUNDS
03158                      && rounding == 1
03159 #endif
03160                             ) {
03161               if (!word1(d) && !(word0(d) & Bndry_mask)
03162 #ifndef Sudden_Underflow
03163                && word0(d) & (Exp_mask & ~Exp_msk1)
03164 #endif
03165                             ) {
03166                      /* The special case */
03167                      b2 += Log2P;
03168                      s2 += Log2P;
03169                      spec_case = 1;
03170                      }
03171               }
03172 
03173        /* Arrange for convenient computation of quotients:
03174         * shift left if necessary so divisor has 4 leading 0 bits.
03175         *
03176         * Perhaps we should just compute leading 28 bits of S once
03177         * and for all and pass them and a shift to quorem, so it
03178         * can do shifts and ors to compute the numerator for q.
03179         */
03180 #ifdef Pack_32
03181        if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
03182               i = 32 - i;
03183 #else
03184        if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
03185               i = 16 - i;
03186 #endif
03187        if (i > 4) {
03188               i -= 4;
03189               b2 += i;
03190               m2 += i;
03191               s2 += i;
03192               }
03193        else if (i < 4) {
03194               i += 28;
03195               b2 += i;
03196               m2 += i;
03197               s2 += i;
03198               }
03199        if (b2 > 0)
03200               b = lshift(b, b2);
03201        if (s2 > 0)
03202               S = lshift(S, s2);
03203        if (k_check) {
03204               if (cmp(b,S) < 0) {
03205                      k--;
03206                      b = multadd(b, 10, 0);      /* we botched the k estimate */
03207                      if (leftright)
03208                             mhi = multadd(mhi, 10, 0);
03209                      ilim = ilim1;
03210                      }
03211               }
03212        if (ilim <= 0 && (mode == 3 || mode == 5)) {
03213               if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
03214                      /* no digits, fcvt style */
03215  no_digits:
03216                      k = -1 - ndigits;
03217                      goto ret;
03218                      }
03219  one_digit:
03220               *s++ = '1';
03221               k++;
03222               goto ret;
03223               }
03224        if (leftright) {
03225               if (m2 > 0)
03226                      mhi = lshift(mhi, m2);
03227 
03228               /* Compute mlo -- check for special case
03229                * that d is a normalized power of 2.
03230                */
03231 
03232               mlo = mhi;
03233               if (spec_case) {
03234                      mhi = Balloc(mhi->k);
03235                      Bcopy(mhi, mlo);
03236                      mhi = lshift(mhi, Log2P);
03237                      }
03238 
03239               for(i = 1;;i++) {
03240                      dig = quorem(b,S) + '0';
03241                      /* Do we yet have the shortest decimal string
03242                       * that will round to d?
03243                       */
03244                      j = cmp(b, mlo);
03245                      delta = diff(S, mhi);
03246                      j1 = delta->sign ? 1 : cmp(b, delta);
03247                      Bfree(delta);
03248 #ifndef ROUND_BIASED
03249                      if (j1 == 0 && mode != 1 && !(word1(d) & 1)
03250 #ifdef Honor_FLT_ROUNDS
03251                             && rounding >= 1
03252 #endif
03253                                                            ) {
03254                             if (dig == '9')
03255                                    goto round_9_up;
03256                             if (j > 0)
03257                                    dig++;
03258 #ifdef SET_INEXACT
03259                             else if (!b->x[0] && b->wds <= 1)
03260                                    inexact = 0;
03261 #endif
03262                             *s++ = dig;
03263                             goto ret;
03264                             }
03265 #endif
03266                      if (j < 0 || j == 0 && mode != 1
03267 #ifndef ROUND_BIASED
03268                                                  && !(word1(d) & 1)
03269 #endif
03270                                    ) {
03271                             if (!b->x[0] && b->wds <= 1) {
03272 #ifdef SET_INEXACT
03273                                    inexact = 0;
03274 #endif
03275                                    goto accept_dig;
03276                                    }
03277 #ifdef Honor_FLT_ROUNDS
03278                             if (mode > 1)
03279                              switch(rounding) {
03280                               case 0: goto accept_dig;
03281                               case 2: goto keep_dig;
03282                               }
03283 #endif /*Honor_FLT_ROUNDS*/
03284                             if (j1 > 0) {
03285                                    b = lshift(b, 1);
03286                                    j1 = cmp(b, S);
03287                                    if ((j1 > 0 || j1 == 0 && dig & 1)
03288                                    && dig++ == '9')
03289                                           goto round_9_up;
03290                                    }
03291  accept_dig:
03292                             *s++ = dig;
03293                             goto ret;
03294                             }
03295                      if (j1 > 0) {
03296 #ifdef Honor_FLT_ROUNDS
03297                             if (!rounding)
03298                                    goto accept_dig;
03299 #endif
03300                             if (dig == '9') { /* possible if i == 1 */
03301  round_9_up:
03302                                    *s++ = '9';
03303                                    goto roundoff;
03304                                    }
03305                             *s++ = dig + 1;
03306                             goto ret;
03307                             }
03308 #ifdef Honor_FLT_ROUNDS
03309  keep_dig:
03310 #endif
03311                      *s++ = dig;
03312                      if (i == ilim)
03313                             break;
03314                      b = multadd(b, 10, 0);
03315                      if (mlo == mhi)
03316                             mlo = mhi = multadd(mhi, 10, 0);
03317                      else {
03318                             mlo = multadd(mlo, 10, 0);
03319                             mhi = multadd(mhi, 10, 0);
03320                             }
03321                      }
03322               }
03323        else
03324               for(i = 1;; i++) {
03325                      *s++ = dig = quorem(b,S) + '0';
03326                      if (!b->x[0] && b->wds <= 1) {
03327 #ifdef SET_INEXACT
03328                             inexact = 0;
03329 #endif
03330                             goto ret;
03331                             }
03332                      if (i >= ilim)
03333                             break;
03334                      b = multadd(b, 10, 0);
03335                      }
03336 
03337        /* Round off last digit */
03338 
03339 #ifdef Honor_FLT_ROUNDS
03340        switch(rounding) {
03341          case 0: goto trimzeros;
03342          case 2: goto roundoff;
03343          }
03344 #endif
03345        b = lshift(b, 1);
03346        j = cmp(b, S);
03347        if (j > 0 || j == 0 && dig & 1) {
03348  roundoff:
03349               while(*--s == '9')
03350                      if (s == s0) {
03351                             k++;
03352                             *s++ = '1';
03353                             goto ret;
03354                             }
03355               ++*s++;
03356               }
03357        else {
03358  trimzeros:
03359               while(*--s == '0');
03360               s++;
03361               }
03362  ret:
03363        Bfree(S);
03364        if (mhi) {
03365               if (mlo && mlo != mhi)
03366                      Bfree(mlo);
03367               Bfree(mhi);
03368               }
03369  ret1:
03370 #ifdef SET_INEXACT
03371        if (inexact) {
03372               if (!oldinexact) {
03373                      word0(d) = Exp_1 + (70 << Exp_shift);
03374                      word1(d) = 0;
03375                      dval(d) += 1.;
03376                      }
03377               }
03378        else if (!oldinexact)
03379               clear_inexact();
03380 #endif
03381        Bfree(b);
03382        *s = 0;
03383        *decpt = k + 1;
03384        if (rve)
03385               *rve = s;
03386        return s0;
03387        }
03388 #ifdef __cplusplus
03389 }
03390 #endif
03391 
03392 PR_IMPLEMENT(PRStatus)
03393 PR_dtoa(PRFloat64 d, PRIntn mode, PRIntn ndigits,
03394        PRIntn *decpt, PRIntn *sign, char **rve, char *buf, PRSize bufsize)
03395 {
03396     char *result;
03397     PRSize resultlen;
03398     PRStatus rv = PR_FAILURE;
03399 
03400     if (!_pr_initialized) _PR_ImplicitInitialization();
03401 
03402     if (mode < 0 || mode > 3) {
03403         PR_SetError(PR_INVALID_ARGUMENT_ERROR, 0);
03404         return rv;
03405     }
03406     result = dtoa(d, mode, ndigits, decpt, sign, rve);
03407     if (!result) {
03408         PR_SetError(PR_OUT_OF_MEMORY_ERROR, 0);
03409         return rv;
03410     }
03411     resultlen = strlen(result)+1;
03412     if (bufsize < resultlen) {
03413         PR_SetError(PR_BUFFER_OVERFLOW_ERROR, 0);
03414     } else {
03415         memcpy(buf, result, resultlen);
03416         if (rve) {
03417             *rve = buf + (*rve - result);
03418         }
03419         rv = PR_SUCCESS;
03420     }
03421     freedtoa(result);
03422     return rv;  
03423 }
03424 
03425 /*
03426 ** conversion routines for floating point
03427 ** prcsn - number of digits of precision to generate floating
03428 ** point value.
03429 ** This should be reparameterized so that you can send in a
03430 **   prcn for the positive and negative ranges.  For now, 
03431 **   conform to the ECMA JavaScript spec which says numbers
03432 **   less than 1e-6 are in scientific notation.
03433 ** Also, the ECMA spec says that there should always be a
03434 **   '+' or '-' after the 'e' in scientific notation
03435 */
03436 PR_IMPLEMENT(void)
03437 PR_cnvtf(char *buf,int bufsz, int prcsn,double fval)
03438 {
03439     PRIntn decpt, sign, numdigits;
03440     char *num, *nump;
03441     char *bufp = buf;
03442     char *endnum;
03443 
03444     /* If anything fails, we store an empty string in 'buf' */
03445     num = (char*)PR_MALLOC(bufsz);
03446     if (num == NULL) {
03447         buf[0] = '\0';
03448         return;
03449     }
03450     /* XXX Why use mode 1? */
03451     if (PR_dtoa(dval(fval),1,prcsn,&decpt,&sign,&endnum,num,bufsz)
03452             == PR_FAILURE) {
03453         buf[0] = '\0';
03454         goto done;
03455     }
03456     numdigits = endnum - num;
03457     nump = num;
03458 
03459     if (sign &&
03460         !(word0(fval) == Sign_bit && word1(fval) == 0) &&
03461         !((word0(fval) & Exp_mask) == Exp_mask &&
03462           (word1(fval) || (word0(fval) & 0xfffff)))) {
03463         *bufp++ = '-';
03464     }
03465 
03466     if (decpt == 9999) {
03467         while ((*bufp++ = *nump++) != 0) {} /* nothing to execute */
03468         goto done;
03469     }
03470 
03471     if (decpt > (prcsn+1) || decpt < -(prcsn-1) || decpt < -5) {
03472         *bufp++ = *nump++;
03473         if (numdigits != 1) {
03474             *bufp++ = '.';
03475         }
03476 
03477         while (*nump != '\0') {
03478             *bufp++ = *nump++;
03479         }
03480         *bufp++ = 'e';
03481         PR_snprintf(bufp, bufsz - (bufp - buf), "%+d", decpt-1);
03482     } else if (decpt >= 0) {
03483         if (decpt == 0) {
03484             *bufp++ = '0';
03485         } else {
03486             while (decpt--) {
03487                 if (*nump != '\0') {
03488                     *bufp++ = *nump++;
03489                 } else {
03490                     *bufp++ = '0';
03491                 }
03492             }
03493         }
03494         if (*nump != '\0') {
03495             *bufp++ = '.';
03496             while (*nump != '\0') {
03497                 *bufp++ = *nump++;
03498             }
03499         }
03500         *bufp++ = '\0';
03501     } else if (decpt < 0) {
03502         *bufp++ = '0';
03503         *bufp++ = '.';
03504         while (decpt++) {
03505             *bufp++ = '0';
03506         }
03507 
03508         while (*nump != '\0') {
03509             *bufp++ = *nump++;
03510         }
03511         *bufp++ = '\0';
03512     }
03513 done:
03514     PR_DELETE(num);
03515 }