My Project
OSdtoa.cpp
Go to the documentation of this file.
1/******************************************************************
2 *
3 * This code is almost entirely copied from David Gay's package ASL,
4 * with the following minor exceptions: At the beginning of the file,
5 * there are some definitions, and after line 47 we changed the names
6 * of four routines in order to avoid conflicts with ASL:
7 *
8 * We changed the names of dtoa to os_dtoa
9 * strtod to os_strtod and
10 * freedtoa to os_freedtoa
11 * gethex to os_gethex
12 *
13 ******************************************************************/
14
21#include "OSConfig.h"
22#include "OSdtoa.h"
23#include "OSParameters.h"
24
25#ifdef WORDS_BIGENDIAN
26#define IEEE_MC68k
27#else
28#define IEEE_8087
29#endif
30
31#define INFNAN_CHECK
32
33#define NO_LONG_LONG
34#define Just_16
35
36#if SIZEOF_LONG == 2*SIZEOF_INT
37#define Long int
38#define Intcast (int)(long)
39#endif
40
47/****************************************************************
48 *
49 * The author of this software is David M. Gay.
50 *
51 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
52 *
53 * Permission to use, copy, modify, and distribute this software for any
54 * purpose without fee is hereby granted, provided that this entire notice
55 * is included in all copies of any software which is or includes a copy
56 * or modification of this software and in all copies of the supporting
57 * documentation for such software.
58 *
59 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
60 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
61 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
62 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
63 *
64 ***************************************************************/
65
66/* Please send bug reports to David M. Gay (dmg at acm dot org,
67 * with " at " changed to "@" and " dot " changed to "."). */
68
69/* On a machine with IEEE extended-precision registers, it is
70 * necessary to specify double-precision (53-bit) rounding precision
71 * before invoking strtod or dtoa. If the machine uses (the equivalent
72 * of) Intel 80x87 arithmetic, the call
73 * _control87(PC_53, MCW_PC);
74 * does this with many compilers. Whether this or another call is
75 * appropriate depends on the compiler; for this to work, it may be
76 * necessary to #include "float.h" or another system-dependent header
77 * file.
78 */
79
80/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
81 * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.)
82 *
83 * This strtod returns a nearest machine number to the input decimal
84 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
85 * broken by the IEEE round-even rule. Otherwise ties are broken by
86 * biased rounding (add half and chop).
87 *
88 * Inspired loosely by William D. Clinger's paper "How to Read Floating
89 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
90 *
91 * Modifications:
92 *
93 * 1. We only require IEEE, IBM, or VAX double-precision
94 * arithmetic (not IEEE double-extended).
95 * 2. We get by with floating-point arithmetic in a case that
96 * Clinger missed -- when we're computing d * 10^n
97 * for a small integer d and the integer n is not too
98 * much larger than 22 (the maximum integer k for which
99 * we can represent 10^k exactly), we may be able to
100 * compute (d*10^k) * 10^(e-k) with just one roundoff.
101 * 3. Rather than a bit-at-a-time adjustment of the binary
102 * result in the hard case, we use floating-point
103 * arithmetic to determine the adjustment to within
104 * one bit; only in really hard cases do we need to
105 * compute a second residual.
106 * 4. Because of 3., we don't need a large table of powers of 10
107 * for ten-to-e (just some small tables, e.g. of 10^k
108 * for 0 <= k <= 22).
109 */
110
111/*
112 * #define IEEE_8087 for IEEE-arithmetic machines where the least
113 * significant byte has the lowest address.
114 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
115 * significant byte has the lowest address.
116 * #define Long int on machines with 32-bit ints and 64-bit longs.
117 * #define IBM for IBM mainframe-style floating-point arithmetic.
118 * #define VAX for VAX-style floating-point arithmetic (D_floating).
119 * #define No_leftright to omit left-right logic in fast floating-point
120 * computation of dtoa. This will cause dtoa modes 4 and 5 to be
121 * treated the same as modes 2 and 3 for some inputs.
122 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
123 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
124 * is also #defined, fegetround() will be queried for the rounding mode.
125 * Note that both FLT_ROUNDS and fegetround() are specified by the C99
126 * standard (and are specified to be consistent, with fesetround()
127 * affecting the value of FLT_ROUNDS), but that some (Linux) systems
128 * do not work correctly in this regard, so using fegetround() is more
129 * portable than using FLT_ROUNDS directly.
130 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
131 * and Honor_FLT_ROUNDS is not #defined.
132 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
133 * that use extended-precision instructions to compute rounded
134 * products and quotients) with IBM.
135 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
136 * that rounds toward +Infinity.
137 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
138 * rounding when the underlying floating-point arithmetic uses
139 * unbiased rounding. This prevent using ordinary floating-point
140 * arithmetic when the result could be computed with one rounding error.
141 * #define Inaccurate_Divide for IEEE-format with correctly rounded
142 * products but inaccurate quotients, e.g., for Intel i860.
143 * #define NO_LONG_LONG on machines that do not have a "long long"
144 * integer type (of >= 64 bits). On such machines, you can
145 * #define Just_16 to store 16 bits per 32-bit Long when doing
146 * high-precision integer arithmetic. Whether this speeds things
147 * up or slows things down depends on the machine and the number
148 * being converted. If long long is available and the name is
149 * something other than "long long", #define Llong to be the name,
150 * and if "unsigned Llong" does not work as an unsigned version of
151 * Llong, #define #ULLong to be the corresponding unsigned type.
152 * #define KR_headers for old-style C function headers.
153 * #define Bad_float_h if your system lacks a float.h or if it does not
154 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
155 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
156 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
157 * if memory is available and otherwise does something you deem
158 * appropriate. If MALLOC is undefined, malloc will be invoked
159 * directly -- and assumed always to succeed. Similarly, if you
160 * want something other than the system's free() to be called to
161 * recycle memory acquired from MALLOC, #define FREE to be the
162 * name of the alternate routine. (FREE or free is only called in
163 * pathological cases, e.g., in a dtoa call after a dtoa return in
164 * mode 3 with thousands of digits requested.)
165 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
166 * memory allocations from a private pool of memory when possible.
167 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
168 * unless #defined to be a different length. This default length
169 * suffices to get rid of MALLOC calls except for unusual cases,
170 * such as decimal-to-binary conversion of a very long string of
171 * digits. The longest string dtoa can return is about 751 bytes
172 * long. For conversions by strtod of strings of 800 digits and
173 * all dtoa conversions in single-threaded executions with 8-byte
174 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
175 * pointers, PRIVATE_MEM >= 7112 appears adequate.
176 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
177 * #defined automatically on IEEE systems. On such systems,
178 * when INFNAN_CHECK is #defined, strtod checks
179 * for Infinity and NaN (case insensitively). On some systems
180 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
181 * appropriately -- to the most significant word of a quiet NaN.
182 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
183 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
184 * strtod also accepts (case insensitively) strings of the form
185 * NaN(x), where x is a string of hexadecimal digits and spaces;
186 * if there is only one string of hexadecimal digits, it is taken
187 * for the 52 fraction bits of the resulting NaN; if there are two
188 * or more strings of hex digits, the first is for the high 20 bits,
189 * the second and subsequent for the low 32 bits, with intervening
190 * white space ignored; but if this results in none of the 52
191 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
192 * and NAN_WORD1 are used instead.
193 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
194 * multiple threads. In this case, you must provide (or suitably
195 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
196 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
197 * in pow5mult, ensures lazy evaluation of only one copy of high
198 * powers of 5; omitting this lock would introduce a small
199 * probability of wasting memory, but would otherwise be harmless.)
200 * You must also invoke freedtoa(s) to free the value s returned by
201 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
202 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
203 * avoids underflows on inputs whose result does not underflow.
204 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
205 * floating-point numbers and flushes underflows to zero rather
206 * than implementing gradual underflow, then you must also #define
207 * Sudden_Underflow.
208 * #define USE_LOCALE to use the current locale's decimal_point value.
209 * #define SET_INEXACT if IEEE arithmetic is being used and extra
210 * computation should be done to set the inexact flag when the
211 * result is inexact and avoid setting inexact when the result
212 * is exact. In this case, dtoa.c must be compiled in
213 * an environment, perhaps provided by #include "dtoa.c" in a
214 * suitable wrapper, that defines two functions,
215 * int get_inexact(void);
216 * void clear_inexact(void);
217 * such that get_inexact() returns a nonzero value if the
218 * inexact bit is already set, and clear_inexact() sets the
219 * inexact bit to 0. When SET_INEXACT is #defined, strtod
220 * also does extra computations to set the underflow and overflow
221 * flags when appropriate (i.e., when the result is tiny and
222 * inexact or when it is a numeric value rounded to +-infinity).
223 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
224 * the result overflows to +-Infinity or underflows to 0.
225 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
226 * values by strtod.
227 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
228 * to disable logic for "fast" testing of very long input strings
229 * to strtod. This testing proceeds by initially truncating the
230 * input string, then if necessary comparing the whole string with
231 * a decimal expansion to decide close cases. This logic is only
232 * used for input more than STRTOD_DIGLIM digits long (default 40).
233 */
234
235#ifndef Long
236#define Long long
237#endif
238#ifndef ULong
239typedef unsigned Long ULong;
240#endif
241
242
243
244#ifdef DEBUG
245#include "stdio.h"
246#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
247#endif
248
249#include "stdlib.h"
250#include "string.h"
251
252#ifdef USE_LOCALE
253#include "locale.h"
254#endif
255
256#ifdef Honor_FLT_ROUNDS
257#ifndef Trust_FLT_ROUNDS
258#include <fenv.h>
259#endif
260#endif
261
262#ifdef MALLOC
263#ifdef KR_headers
264extern char *MALLOC();
265#else
266extern void *MALLOC(size_t);
267#endif
268#else
269#define MALLOC malloc
270#endif
271
272#ifndef Omit_Private_Memory
273#ifndef PRIVATE_MEM
274#define PRIVATE_MEM 2304
275#endif
276#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
278#endif
279
280#undef IEEE_Arith
281#undef Avoid_Underflow
282#ifdef IEEE_MC68k
283#define IEEE_Arith
284#endif
285#ifdef IEEE_8087
286#define IEEE_Arith
287#endif
288
289#ifdef IEEE_Arith
290#ifndef NO_INFNAN_CHECK
291#undef INFNAN_CHECK
292#define INFNAN_CHECK
293#endif
294#else
295#undef INFNAN_CHECK
296#define NO_STRTOD_BIGCOMP
297#endif
298
299#include "errno.h"
300
301#ifdef Bad_float_h
302
303#ifdef IEEE_Arith
304#define DBL_DIG 15
305#define DBL_MAX_10_EXP 308
306#define DBL_MAX_EXP 1024
307#define FLT_RADIX 2
308#endif /*IEEE_Arith*/
309
310#ifdef IBM
311#define DBL_DIG 16
312#define DBL_MAX_10_EXP 75
313#define DBL_MAX_EXP 63
314#define FLT_RADIX 16
315#define DBL_MAX 7.2370055773322621e+75
316#endif
317
318#ifdef VAX
319#define DBL_DIG 16
320#define DBL_MAX_10_EXP 38
321#define DBL_MAX_EXP 127
322#define FLT_RADIX 2
323#define DBL_MAX 1.7014118346046923e+38
324#endif
325
326#ifndef LONG_MAX
327#define LONG_MAX 2147483647
328#endif
329
330#else /* ifndef Bad_float_h */
331#include "float.h"
332#endif /* Bad_float_h */
333
334#ifndef __MATH_H__
335#include "math.h"
336#endif
337
338#ifdef __cplusplus
339extern "C" {
340#endif
341
342#ifndef CONST
343#ifdef KR_headers
344#define CONST /* blank */
345#else
346#define CONST const
347#endif
348#endif
349
350//#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
351//Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
352//#endif
353
354typedef union { double d; ULong L[2]; } U;
355
356#ifdef IEEE_8087
357#define word0(x) (x)->L[1]
358#define word1(x) (x)->L[0]
359#else
360#define word0(x) (x)->L[0]
361#define word1(x) (x)->L[1]
362#endif
363#define dval(x) (x)->d
364
365#ifndef STRTOD_DIGLIM
366#define STRTOD_DIGLIM 40
367#endif
368
369#ifdef DIGLIM_DEBUG
370extern int strtod_diglim;
371#else
372#define strtod_diglim STRTOD_DIGLIM
373#endif
374
375/* The following definition of Storeinc is appropriate for MIPS processors.
376 * An alternative that might be better on some machines is
377 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
378 */
379#if defined(IEEE_8087) + defined(VAX)
380#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
381((unsigned short *)a)[0] = (unsigned short)c, a++)
382#else
383#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
384((unsigned short *)a)[1] = (unsigned short)c, a++)
385#endif
386
387/* #define P DBL_MANT_DIG */
388/* Ten_pmax = floor(P*log(2)/log(5)) */
389/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
390/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
391/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
392
393#ifdef IEEE_Arith
394#define Exp_shift 20
395#define Exp_shift1 20
396#define Exp_msk1 0x100000
397#define Exp_msk11 0x100000
398#define Exp_mask 0x7ff00000
399#define P 53
400#define Nbits 53
401#define Bias 1023
402#define Emax 1023
403#define Emin (-1022)
404#define Exp_1 0x3ff00000
405#define Exp_11 0x3ff00000
406#define Ebits 11
407#define Frac_mask 0xfffff
408#define Frac_mask1 0xfffff
409#define Ten_pmax 22
410#define Bletch 0x10
411#define Bndry_mask 0xfffff
412#define Bndry_mask1 0xfffff
413#define LSB 1
414#define Sign_bit 0x80000000
415#define Log2P 1
416#define Tiny0 0
417#define Tiny1 1
418#define Quick_max 14
419#define Int_max 14
420#ifndef NO_IEEE_Scale
421#define Avoid_Underflow
422#ifdef Flush_Denorm /* debugging option */
423#undef Sudden_Underflow
424#endif
425#endif
426
427#ifndef Flt_Rounds
428#ifdef FLT_ROUNDS
429#define Flt_Rounds FLT_ROUNDS
430#else
431#define Flt_Rounds 1
432#endif
433#endif /*Flt_Rounds*/
434
435#ifdef Honor_FLT_ROUNDS
436#undef Check_FLT_ROUNDS
437#define Check_FLT_ROUNDS
438#else
439#define Rounding Flt_Rounds
440#endif
441
442#else /* ifndef IEEE_Arith */
443#undef Check_FLT_ROUNDS
444#undef Honor_FLT_ROUNDS
445#undef SET_INEXACT
446#undef Sudden_Underflow
447#define Sudden_Underflow
448#ifdef IBM
449#undef Flt_Rounds
450#define Flt_Rounds 0
451#define Exp_shift 24
452#define Exp_shift1 24
453#define Exp_msk1 0x1000000
454#define Exp_msk11 0x1000000
455#define Exp_mask 0x7f000000
456#define P 14
457#define Nbits 56
458#define Bias 65
459#define Emax 248
460#define Emin (-260)
461#define Exp_1 0x41000000
462#define Exp_11 0x41000000
463#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
464#define Frac_mask 0xffffff
465#define Frac_mask1 0xffffff
466#define Bletch 4
467#define Ten_pmax 22
468#define Bndry_mask 0xefffff
469#define Bndry_mask1 0xffffff
470#define LSB 1
471#define Sign_bit 0x80000000
472#define Log2P 4
473#define Tiny0 0x100000
474#define Tiny1 0
475#define Quick_max 14
476#define Int_max 15
477#else /* VAX */
478#undef Flt_Rounds
479#define Flt_Rounds 1
480#define Exp_shift 23
481#define Exp_shift1 7
482#define Exp_msk1 0x80
483#define Exp_msk11 0x800000
484#define Exp_mask 0x7f80
485#define P 56
486#define Nbits 56
487#define Bias 129
488#define Emax 126
489#define Emin (-129)
490#define Exp_1 0x40800000
491#define Exp_11 0x4080
492#define Ebits 8
493#define Frac_mask 0x7fffff
494#define Frac_mask1 0xffff007f
495#define Ten_pmax 24
496#define Bletch 2
497#define Bndry_mask 0xffff007f
498#define Bndry_mask1 0xffff007f
499#define LSB 0x10000
500#define Sign_bit 0x8000
501#define Log2P 1
502#define Tiny0 0x80
503#define Tiny1 0
504#define Quick_max 15
505#define Int_max 15
506#endif /* IBM, VAX */
507#endif /* IEEE_Arith */
508
509#ifndef IEEE_Arith
510#define ROUND_BIASED
511#else
512#ifdef ROUND_BIASED_without_Round_Up
513#undef ROUND_BIASED
514#define ROUND_BIASED
515#endif
516#endif
517
518#ifdef RND_PRODQUOT
519#define rounded_product(a,b) a = rnd_prod(a, b)
520#define rounded_quotient(a,b) a = rnd_quot(a, b)
521#ifdef KR_headers
522extern double rnd_prod(), rnd_quot();
523#else
524extern double rnd_prod(double, double), rnd_quot(double, double);
525#endif
526#else
527#define rounded_product(a,b) a *= b
528#define rounded_quotient(a,b) a /= b
529#endif
530
531#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
532#define Big1 0xffffffff
533
534#ifndef Pack_32
535#define Pack_32
536#endif
537
538typedef struct BCinfo BCinfo;
541
542#ifdef KR_headers
543#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
544#else
545#define FFFFFFFF 0xffffffffUL
546#endif
547
548#ifdef NO_LONG_LONG
549#undef ULLong
550#ifdef Just_16
551#undef Pack_32
552/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
553 * This makes some inner loops simpler and sometimes saves work
554 * during multiplications, but it often seems to make things slightly
555 * slower. Hence the default is now to store 32 bits per Long.
556 */
557#endif
558#else /* long long available */
559#ifndef Llong
560#define Llong long long
561#endif
562#ifndef ULLong
563#define ULLong unsigned Llong
564#endif
565#endif /* NO_LONG_LONG */
566
567#ifndef MULTIPLE_THREADS
568#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
569#define FREE_DTOA_LOCK(n) /*nothing*/
570#endif
571
572#define Kmax 7
573
574#ifdef __cplusplus
575extern "C" double os_strtod(const char *s00, char **se);
576extern "C" char *os_dtoa(double d, int mode, int ndigits,
577 int *decpt, int *sign, char **rve);
578#endif
579
580 struct
581Bigint {
582 struct Bigint *next;
583 int k, maxwds, sign, wds;
585 };
586
587 typedef struct Bigint Bigint;
588
589 static Bigint *freelist[Kmax+1];
590
591 static Bigint *
593#ifdef KR_headers
594 (k) int k;
595#else
596 (int k)
597#endif
598{
599 int x;
600 Bigint *rv;
601#ifndef Omit_Private_Memory
602 unsigned int len;
603#endif
604
606 /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
607 /* but this case seems very unlikely. */
608 if (k <= Kmax && (rv = freelist[k]))
609 freelist[k] = rv->next;
610 else {
611 x = 1 << k;
612#ifdef Omit_Private_Memory
613 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
614#else
615 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
616 /sizeof(double);
617 if (k <= Kmax && pmem_next - private_mem + len <= (long)PRIVATE_mem) {
618 rv = (Bigint*)pmem_next;
619 pmem_next += len;
620 }
621 else
622 rv = (Bigint*)MALLOC(len*sizeof(double));
623#endif
624 rv->k = k;
625 rv->maxwds = x;
626 }
628 rv->sign = rv->wds = 0;
629 return rv;
630 }
631
632 static void
634#ifdef KR_headers
635 (v) Bigint *v;
636#else
637 (Bigint *v)
638#endif
639{
640 if (v) {
641 if (v->k > Kmax)
642#ifdef FREE
643 FREE((void*)v);
644#else
645 free((void*)v);
646#endif
647 else {
649 v->next = freelist[v->k];
650 freelist[v->k] = v;
652 }
653 }
654 }
655
656#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
657y->wds*sizeof(Long) + 2*sizeof(int))
658
659 static Bigint *
661#ifdef KR_headers
662 (b, m, a) Bigint *b; int m, a;
663#else
664 (Bigint *b, int m, int a) /* multiply by m and add a */
665#endif
666{
667 int i, wds;
668#ifdef ULLong
669 ULong *x;
670 ULLong carry, y;
671#else
672 ULong carry, *x, y;
673#ifdef Pack_32
674 ULong xi, z;
675#endif
676#endif
677 Bigint *b1;
678
679 wds = b->wds;
680 x = b->x;
681 i = 0;
682 carry = a;
683 do {
684#ifdef ULLong
685 y = *x * (ULLong)m + carry;
686 carry = y >> 32;
687 *x++ = y & FFFFFFFF;
688#else
689#ifdef Pack_32
690 xi = *x;
691 y = (xi & 0xffff) * m + carry;
692 z = (xi >> 16) * m + (y >> 16);
693 carry = z >> 16;
694 *x++ = (z << 16) + (y & 0xffff);
695#else
696 y = *x * m + carry;
697 carry = y >> 16;
698 *x++ = y & 0xffff;
699#endif
700#endif
701 }
702 while(++i < wds);
703 if (carry) {
704 if (wds >= b->maxwds) {
705 b1 = Balloc(b->k+1);
706 Bcopy(b1, b);
707 Bfree(b);
708 b = b1;
709 }
710 b->x[wds++] = carry;
711 b->wds = wds;
712 }
713 return b;
714 }
715
716 static Bigint *
718#ifdef KR_headers
719 (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
720#else
721 (const char *s, int nd0, int nd, ULong y9, int dplen)
722#endif
723{
724 Bigint *b;
725 int i, k;
726 Long x, y;
727
728 x = (nd + 8) / 9;
729 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
730#ifdef Pack_32
731 b = Balloc(k);
732 b->x[0] = y9;
733 b->wds = 1;
734#else
735 b = Balloc(k+1);
736 b->x[0] = y9 & 0xffff;
737 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
738#endif
739
740 i = 9;
741 if (9 < nd0) {
742 s += 9;
743 do b = multadd(b, 10, *s++ - '0');
744 while(++i < nd0);
745 s += dplen;
746 }
747 else
748 s += dplen + 9;
749 for(; i < nd; i++)
750 b = multadd(b, 10, *s++ - '0');
751 return b;
752 }
753
754 static int
756#ifdef KR_headers
757 (x) ULong x;
758#else
759 (ULong x)
760#endif
761{
762 int k = 0;
763
764 if (!(x & 0xffff0000)) {
765 k = 16;
766 x <<= 16;
767 }
768 if (!(x & 0xff000000)) {
769 k += 8;
770 x <<= 8;
771 }
772 if (!(x & 0xf0000000)) {
773 k += 4;
774 x <<= 4;
775 }
776 if (!(x & 0xc0000000)) {
777 k += 2;
778 x <<= 2;
779 }
780 if (!(x & 0x80000000)) {
781 k++;
782 if (!(x & 0x40000000))
783 return 32;
784 }
785 return k;
786 }
787
788 static int
790#ifdef KR_headers
791 (y) ULong *y;
792#else
793 (ULong *y)
794#endif
795{
796 int k;
797 ULong x = *y;
798
799 if (x & 7) {
800 if (x & 1)
801 return 0;
802 if (x & 2) {
803 *y = x >> 1;
804 return 1;
805 }
806 *y = x >> 2;
807 return 2;
808 }
809 k = 0;
810 if (!(x & 0xffff)) {
811 k = 16;
812 x >>= 16;
813 }
814 if (!(x & 0xff)) {
815 k += 8;
816 x >>= 8;
817 }
818 if (!(x & 0xf)) {
819 k += 4;
820 x >>= 4;
821 }
822 if (!(x & 0x3)) {
823 k += 2;
824 x >>= 2;
825 }
826 if (!(x & 1)) {
827 k++;
828 x >>= 1;
829 if (!x)
830 return 32;
831 }
832 *y = x;
833 return k;
834 }
835
836 static Bigint *
838#ifdef KR_headers
839 (i) int i;
840#else
841 (int i)
842#endif
843{
844 Bigint *b;
845
846 b = Balloc(1);
847 b->x[0] = i;
848 b->wds = 1;
849 return b;
850 }
851
852 static Bigint *
854#ifdef KR_headers
855 (a, b) Bigint *a, *b;
856#else
857 (Bigint *a, Bigint *b)
858#endif
859{
860 Bigint *c;
861 int k, wa, wb, wc;
862 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
863 ULong y;
864#ifdef ULLong
865 ULLong carry, z;
866#else
867 ULong carry, z;
868#ifdef Pack_32
869 ULong z2;
870#endif
871#endif
872
873 if (a->wds < b->wds) {
874 c = a;
875 a = b;
876 b = c;
877 }
878 k = a->k;
879 wa = a->wds;
880 wb = b->wds;
881 wc = wa + wb;
882 if (wc > a->maxwds)
883 k++;
884 c = Balloc(k);
885 for(x = c->x, xa = x + wc; x < xa; x++)
886 *x = 0;
887 xa = a->x;
888 xae = xa + wa;
889 xb = b->x;
890 xbe = xb + wb;
891 xc0 = c->x;
892#ifdef ULLong
893 for(; xb < xbe; xc0++) {
894 if ((y = *xb++)) {
895 x = xa;
896 xc = xc0;
897 carry = 0;
898 do {
899 z = *x++ * (ULLong)y + *xc + carry;
900 carry = z >> 32;
901 *xc++ = z & FFFFFFFF;
902 }
903 while(x < xae);
904 *xc = carry;
905 }
906 }
907#else
908#ifdef Pack_32
909 for(; xb < xbe; xb++, xc0++) {
910 if (y = *xb & 0xffff) {
911 x = xa;
912 xc = xc0;
913 carry = 0;
914 do {
915 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
916 carry = z >> 16;
917 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
918 carry = z2 >> 16;
919 Storeinc(xc, z2, z);
920 }
921 while(x < xae);
922 *xc = carry;
923 }
924 if (y = *xb >> 16) {
925 x = xa;
926 xc = xc0;
927 carry = 0;
928 z2 = *xc;
929 do {
930 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
931 carry = z >> 16;
932 Storeinc(xc, z, z2);
933 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
934 carry = z2 >> 16;
935 }
936 while(x < xae);
937 *xc = z2;
938 }
939 }
940#else
941 for(; xb < xbe; xc0++) {
942 if (y = *xb++) {
943 x = xa;
944 xc = xc0;
945 carry = 0;
946 do {
947 z = *x++ * y + *xc + carry;
948 carry = z >> 16;
949 *xc++ = z & 0xffff;
950 }
951 while(x < xae);
952 *xc = carry;
953 }
954 }
955#endif
956#endif
957 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
958 c->wds = wc;
959 return c;
960 }
961
962 static Bigint *p5s;
963
964 static Bigint *
966#ifdef KR_headers
967 (b, k) Bigint *b; int k;
968#else
969 (Bigint *b, int k)
970#endif
971{
972 Bigint *b1, *p5, *p51;
973 int i;
974 static int p05[3] = { 5, 25, 125 };
975
976 if ((i = k & 3))
977 b = multadd(b, p05[i-1], 0);
978
979 if (!(k >>= 2))
980 return b;
981 if (!(p5 = p5s)) {
982 /* first time */
983#ifdef MULTIPLE_THREADS
985 if (!(p5 = p5s)) {
986 p5 = p5s = i2b(625);
987 p5->next = 0;
988 }
990#else
991 p5 = p5s = i2b(625);
992 p5->next = 0;
993#endif
994 }
995 for(;;) {
996 if (k & 1) {
997 b1 = mult(b, p5);
998 Bfree(b);
999 b = b1;
1000 }
1001 if (!(k >>= 1))
1002 break;
1003 if (!(p51 = p5->next)) {
1004#ifdef MULTIPLE_THREADS
1006 if (!(p51 = p5->next)) {
1007 p51 = p5->next = mult(p5,p5);
1008 p51->next = 0;
1009 }
1010 FREE_DTOA_LOCK(1);
1011#else
1012 p51 = p5->next = mult(p5,p5);
1013 p51->next = 0;
1014#endif
1015 }
1016 p5 = p51;
1017 }
1018 return b;
1019 }
1020
1021 static Bigint *
1023#ifdef KR_headers
1024 (b, k) Bigint *b; int k;
1025#else
1026 (Bigint *b, int k)
1027#endif
1028{
1029 int i, k1, n, n1;
1030 Bigint *b1;
1031 ULong *x, *x1, *xe, z;
1032
1033#ifdef Pack_32
1034 n = k >> 5;
1035#else
1036 n = k >> 4;
1037#endif
1038 k1 = b->k;
1039 n1 = n + b->wds + 1;
1040 for(i = b->maxwds; n1 > i; i <<= 1)
1041 k1++;
1042 b1 = Balloc(k1);
1043 x1 = b1->x;
1044 for(i = 0; i < n; i++)
1045 *x1++ = 0;
1046 x = b->x;
1047 xe = x + b->wds;
1048#ifdef Pack_32
1049 if (k &= 0x1f) {
1050 k1 = 32 - k;
1051 z = 0;
1052 do {
1053 *x1++ = *x << k | z;
1054 z = *x++ >> k1;
1055 }
1056 while(x < xe);
1057 if ((*x1 = z))
1058 ++n1;
1059 }
1060#else
1061 if (k &= 0xf) {
1062 k1 = 16 - k;
1063 z = 0;
1064 do {
1065 *x1++ = *x << k & 0xffff | z;
1066 z = *x++ >> k1;
1067 }
1068 while(x < xe);
1069 if (*x1 = z)
1070 ++n1;
1071 }
1072#endif
1073 else do
1074 *x1++ = *x++;
1075 while(x < xe);
1076 b1->wds = n1 - 1;
1077 Bfree(b);
1078 return b1;
1079 }
1080
1081 static int
1083#ifdef KR_headers
1084 (a, b) Bigint *a, *b;
1085#else
1086 (Bigint *a, Bigint *b)
1087#endif
1088{
1089 ULong *xa, *xa0, *xb, *xb0;
1090 int i, j;
1091
1092 i = a->wds;
1093 j = b->wds;
1094#ifdef DEBUG
1095 if (i > 1 && !a->x[i-1])
1096 Bug("cmp called with a->x[a->wds-1] == 0");
1097 if (j > 1 && !b->x[j-1])
1098 Bug("cmp called with b->x[b->wds-1] == 0");
1099#endif
1100 if (i -= j)
1101 return i;
1102 xa0 = a->x;
1103 xa = xa0 + j;
1104 xb0 = b->x;
1105 xb = xb0 + j;
1106 for(;;) {
1107 if (*--xa != *--xb)
1108 return *xa < *xb ? -1 : 1;
1109 if (xa <= xa0)
1110 break;
1111 }
1112 return 0;
1113 }
1114
1115 static Bigint *
1117#ifdef KR_headers
1118 (a, b) Bigint *a, *b;
1119#else
1120 (Bigint *a, Bigint *b)
1121#endif
1122{
1123 Bigint *c;
1124 int i, wa, wb;
1125 ULong *xa, *xae, *xb, *xbe, *xc;
1126#ifdef ULLong
1127 ULLong borrow, y;
1128#else
1129 ULong borrow, y;
1130#ifdef Pack_32
1131 ULong z;
1132#endif
1133#endif
1134
1135 i = cmp(a,b);
1136 if (!i) {
1137 c = Balloc(0);
1138 c->wds = 1;
1139 c->x[0] = 0;
1140 return c;
1141 }
1142 if (i < 0) {
1143 c = a;
1144 a = b;
1145 b = c;
1146 i = 1;
1147 }
1148 else
1149 i = 0;
1150 c = Balloc(a->k);
1151 c->sign = i;
1152 wa = a->wds;
1153 xa = a->x;
1154 xae = xa + wa;
1155 wb = b->wds;
1156 xb = b->x;
1157 xbe = xb + wb;
1158 xc = c->x;
1159 borrow = 0;
1160#ifdef ULLong
1161 do {
1162 y = (ULLong)*xa++ - *xb++ - borrow;
1163 borrow = y >> 32 & (ULong)1;
1164 *xc++ = y & FFFFFFFF;
1165 }
1166 while(xb < xbe);
1167 while(xa < xae) {
1168 y = *xa++ - borrow;
1169 borrow = y >> 32 & (ULong)1;
1170 *xc++ = y & FFFFFFFF;
1171 }
1172#else
1173#ifdef Pack_32
1174 do {
1175 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1176 borrow = (y & 0x10000) >> 16;
1177 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1178 borrow = (z & 0x10000) >> 16;
1179 Storeinc(xc, z, y);
1180 }
1181 while(xb < xbe);
1182 while(xa < xae) {
1183 y = (*xa & 0xffff) - borrow;
1184 borrow = (y & 0x10000) >> 16;
1185 z = (*xa++ >> 16) - borrow;
1186 borrow = (z & 0x10000) >> 16;
1187 Storeinc(xc, z, y);
1188 }
1189#else
1190 do {
1191 y = *xa++ - *xb++ - borrow;
1192 borrow = (y & 0x10000) >> 16;
1193 *xc++ = y & 0xffff;
1194 }
1195 while(xb < xbe);
1196 while(xa < xae) {
1197 y = *xa++ - borrow;
1198 borrow = (y & 0x10000) >> 16;
1199 *xc++ = y & 0xffff;
1200 }
1201#endif
1202#endif
1203 while(!*--xc)
1204 wa--;
1205 c->wds = wa;
1206 return c;
1207 }
1208
1209 static double
1211#ifdef KR_headers
1212 (x) U *x;
1213#else
1214 (U *x)
1215#endif
1216{
1217 Long L;
1218 U u;
1219
1220 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1221#ifndef Avoid_Underflow
1222#ifndef Sudden_Underflow
1223 if (L > 0) {
1224#endif
1225#endif
1226#ifdef IBM
1227 L |= Exp_msk1 >> 4;
1228#endif
1229 word0(&u) = L;
1230 word1(&u) = 0;
1231#ifndef Avoid_Underflow
1232#ifndef Sudden_Underflow
1233 }
1234 else {
1235 L = -L >> Exp_shift;
1236 if (L < Exp_shift) {
1237 word0(&u) = 0x80000 >> L;
1238 word1(&u) = 0;
1239 }
1240 else {
1241 word0(&u) = 0;
1242 L -= Exp_shift;
1243 word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1244 }
1245 }
1246#endif
1247#endif
1248 return dval(&u);
1249 }
1250
1251 static double
1253#ifdef KR_headers
1254 (a, e) Bigint *a; int *e;
1255#else
1256 (Bigint *a, int *e)
1257#endif
1258{
1259 ULong *xa, *xa0, w, y, z;
1260 int k;
1261 U d;
1262#ifdef VAX
1263 ULong d0, d1;
1264#else
1265#define d0 word0(&d)
1266#define d1 word1(&d)
1267#endif
1268
1269 xa0 = a->x;
1270 xa = xa0 + a->wds;
1271 y = *--xa;
1272#ifdef DEBUG
1273 if (!y) Bug("zero y in b2d");
1274#endif
1275 k = hi0bits(y);
1276 *e = 32 - k;
1277#ifdef Pack_32
1278 if (k < Ebits) {
1279 d0 = Exp_1 | y >> (Ebits - k);
1280 w = xa > xa0 ? *--xa : 0;
1281 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1282 goto ret_d;
1283 }
1284 z = xa > xa0 ? *--xa : 0;
1285 if (k -= Ebits) {
1286 d0 = Exp_1 | y << k | z >> (32 - k);
1287 y = xa > xa0 ? *--xa : 0;
1288 d1 = z << k | y >> (32 - k);
1289 }
1290 else {
1291 d0 = Exp_1 | y;
1292 d1 = z;
1293 }
1294#else
1295 if (k < Ebits + 16) {
1296 z = xa > xa0 ? *--xa : 0;
1297 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1298 w = xa > xa0 ? *--xa : 0;
1299 y = xa > xa0 ? *--xa : 0;
1300 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1301 goto ret_d;
1302 }
1303 z = xa > xa0 ? *--xa : 0;
1304 w = xa > xa0 ? *--xa : 0;
1305 k -= Ebits + 16;
1306 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1307 y = xa > xa0 ? *--xa : 0;
1308 d1 = w << k + 16 | y << k;
1309#endif
1310 ret_d:
1311#ifdef VAX
1312 word0(&d) = d0 >> 16 | d0 << 16;
1313 word1(&d) = d1 >> 16 | d1 << 16;
1314#else
1315#undef d0
1316#undef d1
1317#endif
1318 return dval(&d);
1319 }
1320
1321 static Bigint *
1323#ifdef KR_headers
1324 (d, e, bits) U *d; int *e, *bits;
1325#else
1326 (U *d, int *e, int *bits)
1327#endif
1328{
1329 Bigint *b;
1330 int de, k;
1331 ULong *x, y, z;
1332#ifndef Sudden_Underflow
1333 int i;
1334#endif
1335#ifdef VAX
1336 ULong d0, d1;
1337 d0 = word0(d) >> 16 | word0(d) << 16;
1338 d1 = word1(d) >> 16 | word1(d) << 16;
1339#else
1340#define d0 word0(d)
1341#define d1 word1(d)
1342#endif
1343
1344#ifdef Pack_32
1345 b = Balloc(1);
1346#else
1347 b = Balloc(2);
1348#endif
1349 x = b->x;
1350
1351 z = d0 & Frac_mask;
1352 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1353#ifdef Sudden_Underflow
1354 de = (int)(d0 >> Exp_shift);
1355#ifndef IBM
1356 z |= Exp_msk11;
1357#endif
1358#else
1359 if ((de = (int)(d0 >> Exp_shift)))
1360 z |= Exp_msk1;
1361#endif
1362#ifdef Pack_32
1363 if ((y = d1)) {
1364 if ((k = lo0bits(&y))) {
1365 x[0] = y | z << (32 - k);
1366 z >>= k;
1367 }
1368 else
1369 x[0] = y;
1370#ifndef Sudden_Underflow
1371 i =
1372#endif
1373 b->wds = (x[1] = z) ? 2 : 1;
1374 }
1375 else {
1376 k = lo0bits(&z);
1377 x[0] = z;
1378#ifndef Sudden_Underflow
1379 i =
1380#endif
1381 b->wds = 1;
1382 k += 32;
1383 }
1384#else
1385 if (y = d1) {
1386 if (k = lo0bits(&y))
1387 if (k >= 16) {
1388 x[0] = y | z << 32 - k & 0xffff;
1389 x[1] = z >> k - 16 & 0xffff;
1390 x[2] = z >> k;
1391 i = 2;
1392 }
1393 else {
1394 x[0] = y & 0xffff;
1395 x[1] = y >> 16 | z << 16 - k & 0xffff;
1396 x[2] = z >> k & 0xffff;
1397 x[3] = z >> k+16;
1398 i = 3;
1399 }
1400 else {
1401 x[0] = y & 0xffff;
1402 x[1] = y >> 16;
1403 x[2] = z & 0xffff;
1404 x[3] = z >> 16;
1405 i = 3;
1406 }
1407 }
1408 else {
1409#ifdef DEBUG
1410 if (!z)
1411 Bug("Zero passed to d2b");
1412#endif
1413 k = lo0bits(&z);
1414 if (k >= 16) {
1415 x[0] = z;
1416 i = 0;
1417 }
1418 else {
1419 x[0] = z & 0xffff;
1420 x[1] = z >> 16;
1421 i = 1;
1422 }
1423 k += 32;
1424 }
1425 while(!x[i])
1426 --i;
1427 b->wds = i + 1;
1428#endif
1429#ifndef Sudden_Underflow
1430 if (de) {
1431#endif
1432#ifdef IBM
1433 *e = (de - Bias - (P-1) << 2) + k;
1434 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1435#else
1436 *e = de - Bias - (P-1) + k;
1437 *bits = P - k;
1438#endif
1439#ifndef Sudden_Underflow
1440 }
1441 else {
1442 *e = de - Bias - (P-1) + 1 + k;
1443#ifdef Pack_32
1444 *bits = 32*i - hi0bits(x[i-1]);
1445#else
1446 *bits = (i+2)*16 - hi0bits(x[i]);
1447#endif
1448 }
1449#endif
1450 return b;
1451 }
1452#undef d0
1453#undef d1
1454
1455 static double
1457#ifdef KR_headers
1458 (a, b) Bigint *a, *b;
1459#else
1460 (Bigint *a, Bigint *b)
1461#endif
1462{
1463 U da, db;
1464 int k, ka, kb;
1465
1466 dval(&da) = b2d(a, &ka);
1467 dval(&db) = b2d(b, &kb);
1468#ifdef Pack_32
1469 k = ka - kb + 32*(a->wds - b->wds);
1470#else
1471 k = ka - kb + 16*(a->wds - b->wds);
1472#endif
1473#ifdef IBM
1474 if (k > 0) {
1475 word0(&da) += (k >> 2)*Exp_msk1;
1476 if (k &= 3)
1477 dval(&da) *= 1 << k;
1478 }
1479 else {
1480 k = -k;
1481 word0(&db) += (k >> 2)*Exp_msk1;
1482 if (k &= 3)
1483 dval(&db) *= 1 << k;
1484 }
1485#else
1486 if (k > 0)
1487 word0(&da) += k*Exp_msk1;
1488 else {
1489 k = -k;
1490 word0(&db) += k*Exp_msk1;
1491 }
1492#endif
1493 return dval(&da) / dval(&db);
1494 }
1495
1496 static CONST double
1497tens[] = {
1498 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1499 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1500 1e20, 1e21, 1e22
1501#ifdef VAX
1502 , 1e23, 1e24
1503#endif
1504 };
1505
1506 static CONST double
1507#ifdef IEEE_Arith
1508bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1509static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1510#ifdef Avoid_Underflow
1511 9007199254740992.*9007199254740992.e-256
1512 /* = 2^106 * 1e-256 */
1513#else
1514 1e-256
1515#endif
1516 };
1517/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1518/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1519#define Scale_Bit 0x10
1520#define n_bigtens 5
1521#else
1522#ifdef IBM
1523bigtens[] = { 1e16, 1e32, 1e64 };
1524static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1525#define n_bigtens 3
1526#else
1527bigtens[] = { 1e16, 1e32 };
1528static CONST double tinytens[] = { 1e-16, 1e-32 };
1529#define n_bigtens 2
1530#endif
1531#endif
1532
1533#undef Need_Hexdig
1534#ifdef INFNAN_CHECK
1535#ifndef No_Hex_NaN
1536#define Need_Hexdig
1537#endif
1538#endif
1539
1540#ifndef Need_Hexdig
1541#ifndef NO_HEX_FP
1542#define Need_Hexdig
1543#endif
1544#endif
1545
1546#ifdef Need_Hexdig /*{*/
1547#if 0
1548static unsigned char hexdig[256];
1549
1550 static void
1551htinit(unsigned char *h, unsigned char *s, int inc)
1552{
1553 int i, j;
1554 for(i = 0; (j = s[i]) !=0; i++)
1555 h[j] = i + inc;
1556 }
1557
1558 static void
1559hexdig_init(void) /* Use of hexdig_init omitted 20121220 to avoid a */
1560 /* race condition when multiple threads are used. */
1561{
1562#define USC (unsigned char *)
1563 htinit(hexdig, USC "0123456789", 0x10);
1564 htinit(hexdig, USC "abcdef", 0x10 + 10);
1565 htinit(hexdig, USC "ABCDEF", 0x10 + 10);
1566 }
1567#else
1568static unsigned char hexdig[256] = {
1569 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1570 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1571 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1572 16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
1573 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1574 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1575 0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1576 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1577 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1578 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1579 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1580 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1581 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1582 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1583 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1584 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
1585 };
1586#endif
1587#endif /* } Need_Hexdig */
1588
1589#ifdef INFNAN_CHECK
1590
1591#ifndef NAN_WORD0
1592#define NAN_WORD0 0x7ff80000
1593#endif
1594
1595#ifndef NAN_WORD1
1596#define NAN_WORD1 0
1597#endif
1598
1599 static int
1601#ifdef KR_headers
1602 (sp, t) char **sp, *t;
1603#else
1604 (const char **sp, const char *t)
1605#endif
1606{
1607 int c, d;
1608 CONST char *s = *sp;
1609
1610 while((d = *t++)) {
1611 if ((c = *++s) >= 'A' && c <= 'Z')
1612 c += 'a' - 'A';
1613 if (c != d)
1614 return 0;
1615 }
1616 *sp = s + 1;
1617 return 1;
1618 }
1619
1620#ifndef No_Hex_NaN
1621 static void
1623#ifdef KR_headers
1624 (rvp, sp) U *rvp; CONST char **sp;
1625#else
1626 (U *rvp, const char **sp)
1627#endif
1628{
1629 ULong c, x[2];
1630 CONST char *s;
1631 int c1, havedig, udx0, xshift;
1632
1633 /**** if (!hexdig['0']) hexdig_init(); ****/
1634 x[0] = x[1] = 0;
1635 havedig = xshift = 0;
1636 udx0 = 1;
1637 s = *sp;
1638 /* allow optional initial 0x or 0X */
1639 while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1640 ++s;
1641 if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1642 s += 2;
1643 while((c = *(CONST unsigned char*)++s)) {
1644 if ((c1 = hexdig[c]))
1645 c = c1 & 0xf;
1646 else if (c <= ' ') {
1647 if (udx0 && havedig) {
1648 udx0 = 0;
1649 xshift = 1;
1650 }
1651 continue;
1652 }
1653#ifdef GDTOA_NON_PEDANTIC_NANCHECK
1654 else if (/*(*/ c == ')' && havedig) {
1655 *sp = s + 1;
1656 break;
1657 }
1658 else
1659 return; /* invalid form: don't change *sp */
1660#else
1661 else {
1662 do {
1663 if (/*(*/ c == ')') {
1664 *sp = s + 1;
1665 break;
1666 }
1667 } while((c = *++s));
1668 break;
1669 }
1670#endif
1671 havedig = 1;
1672 if (xshift) {
1673 xshift = 0;
1674 x[0] = x[1];
1675 x[1] = 0;
1676 }
1677 if (udx0)
1678 x[0] = (x[0] << 4) | (x[1] >> 28);
1679 x[1] = (x[1] << 4) | c;
1680 }
1681 if ((x[0] &= 0xfffff) || x[1]) {
1682 word0(rvp) = Exp_mask | x[0];
1683 word1(rvp) = x[1];
1684 }
1685 }
1686#endif /*No_Hex_NaN*/
1687#endif /* INFNAN_CHECK */
1688
1689#ifdef Pack_32
1690#define ULbits 32
1691#define kshift 5
1692#define kmask 31
1693#else
1694#define ULbits 16
1695#define kshift 4
1696#define kmask 15
1697#endif
1698
1699#if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
1700 static Bigint *
1701#ifdef KR_headers
1702increment(b) Bigint *b;
1703#else
1705#endif
1706{
1707 ULong *x, *xe;
1708 Bigint *b1;
1709
1710 x = b->x;
1711 xe = x + b->wds;
1712 do {
1713 if (*x < (ULong)0xffffffffL) {
1714 ++*x;
1715 return b;
1716 }
1717 *x++ = 0;
1718 } while(x < xe);
1719 {
1720 if (b->wds >= b->maxwds) {
1721 b1 = Balloc(b->k+1);
1722 Bcopy(b1,b);
1723 Bfree(b);
1724 b = b1;
1725 }
1726 b->x[b->wds++] = 1;
1727 }
1728 return b;
1729 }
1730
1731#endif /*}*/
1732
1733#ifndef NO_HEX_FP /*{*/
1734
1735 static void
1736#ifdef KR_headers
1737rshift(b, k) Bigint *b; int k;
1738#else
1740#endif
1741{
1742 ULong *x, *x1, *xe, y;
1743 int n;
1744
1745 x = x1 = b->x;
1746 n = k >> kshift;
1747 if (n < b->wds) {
1748 xe = x + b->wds;
1749 x += n;
1750 if (k &= kmask) {
1751 n = 32 - k;
1752 y = *x++ >> k;
1753 while(x < xe) {
1754 *x1++ = (y | (*x << n)) & 0xffffffff;
1755 y = *x++ >> k;
1756 }
1757 if ((*x1 = y) !=0)
1758 x1++;
1759 }
1760 else
1761 while(x < xe)
1762 *x1++ = *x++;
1763 }
1764 if ((b->wds = x1 - b->x) == 0)
1765 b->x[0] = 0;
1766 }
1767
1768 static ULong
1769#ifdef KR_headers
1770any_on(b, k) Bigint *b; int k;
1771#else
1773#endif
1774{
1775 int n, nwds;
1776 ULong *x, *x0, x1, x2;
1777
1778 x = b->x;
1779 nwds = b->wds;
1780 n = k >> kshift;
1781 if (n > nwds)
1782 n = nwds;
1783 else if (n < nwds && (k &= kmask)) {
1784 x1 = x2 = x[n];
1785 x1 >>= k;
1786 x1 <<= k;
1787 if (x1 != x2)
1788 return 1;
1789 }
1790 x0 = x;
1791 x += n;
1792 while(x > x0)
1793 if (*--x)
1794 return 1;
1795 return 0;
1796 }
1797
1798enum { /* rounding values: same as FLT_ROUNDS */
1802 Round_down = 3
1804
1805 void
1806#ifdef KR_headers
1807os_gethex(sp, rvp, rounding, sign)
1808 CONST char **sp; U *rvp; int rounding, sign;
1809#else
1810os_gethex( CONST char **sp, U *rvp, int rounding, int sign)
1811#endif
1812{
1813 Bigint *b;
1814 CONST unsigned char *decpt, *s0, *s, *s1;
1815 Long e, e1;
1816 ULong L, lostbits, *x;
1817 int big, denorm, esign, havedig, k, n, nbits, up, zret;
1818#ifdef IBM
1819 int j;
1820#endif
1821 enum {
1822#ifdef IEEE_Arith /*{{*/
1823 emax = 0x7fe - Bias - P + 1,
1824 emin = Emin - P + 1
1825#else /*}{*/
1826 emin = Emin - P,
1827#ifdef VAX
1828 emax = 0x7ff - Bias - P + 1
1829#endif
1830#ifdef IBM
1831 emax = 0x7f - Bias - P
1832#endif
1833#endif /*}}*/
1834 };
1835#ifdef USE_LOCALE
1836 int i;
1837#ifdef NO_LOCALE_CACHE
1838 const unsigned char *decimalpoint = (unsigned char*)
1839 localeconv()->decimal_point;
1840#else
1841 const unsigned char *decimalpoint;
1842 static unsigned char *decimalpoint_cache;
1843 if (!(s0 = decimalpoint_cache)) {
1844 s0 = (unsigned char*)localeconv()->decimal_point;
1845 if ((decimalpoint_cache = (unsigned char*)
1846 MALLOC(strlen((CONST char*)s0) + 1))) {
1847 strcpy((char*)decimalpoint_cache, (CONST char*)s0);
1848 s0 = decimalpoint_cache;
1849 }
1850 }
1851 decimalpoint = s0;
1852#endif
1853#endif
1854
1855 /**** if (!hexdig['0']) hexdig_init(); ****/
1856 havedig = 0;
1857 s0 = *(CONST unsigned char **)sp + 2;
1858 while(s0[havedig] == '0')
1859 havedig++;
1860 s0 += havedig;
1861 s = s0;
1862 decpt = 0;
1863 zret = 0;
1864 e = 0;
1865 if (hexdig[*s])
1866 havedig++;
1867 else {
1868 zret = 1;
1869#ifdef USE_LOCALE
1870 for(i = 0; decimalpoint[i]; ++i) {
1871 if (s[i] != decimalpoint[i])
1872 goto pcheck;
1873 }
1874 decpt = s += i;
1875#else
1876 if (*s != '.')
1877 goto pcheck;
1878 decpt = ++s;
1879#endif
1880 if (!hexdig[*s])
1881 goto pcheck;
1882 while(*s == '0')
1883 s++;
1884 if (hexdig[*s])
1885 zret = 0;
1886 havedig = 1;
1887 s0 = s;
1888 }
1889 while(hexdig[*s])
1890 s++;
1891#ifdef USE_LOCALE
1892 if (*s == *decimalpoint && !decpt) {
1893 for(i = 1; decimalpoint[i]; ++i) {
1894 if (s[i] != decimalpoint[i])
1895 goto pcheck;
1896 }
1897 decpt = s += i;
1898#else
1899 if (*s == '.' && !decpt) {
1900 decpt = ++s;
1901#endif
1902 while(hexdig[*s])
1903 s++;
1904 }/*}*/
1905 if (decpt)
1906 e = -(((Long)(s-decpt)) << 2);
1907 pcheck:
1908 s1 = s;
1909 big = esign = 0;
1910 switch(*s) {
1911 case 'p':
1912 case 'P':
1913 switch(*++s) {
1914 case '-':
1915 esign = 1;
1916 /* no break */
1917 case '+':
1918 s++;
1919 }
1920 if ((n = hexdig[*s]) == 0 || n > 0x19) {
1921 s = s1;
1922 break;
1923 }
1924 e1 = n - 0x10;
1925 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
1926 if (e1 & 0xf8000000)
1927 big = 1;
1928 e1 = 10*e1 + n - 0x10;
1929 }
1930 if (esign)
1931 e1 = -e1;
1932 e += e1;
1933 }
1934 *sp = (char*)s;
1935 if (!havedig)
1936 *sp = (char*)s0 - 1;
1937 if (zret)
1938 goto retz1;
1939 if (big) {
1940 if (esign) {
1941#ifdef IEEE_Arith
1942 switch(rounding) {
1943 case Round_up:
1944 if (sign)
1945 break;
1946 goto ret_tiny;
1947 case Round_down:
1948 if (!sign)
1949 break;
1950 goto ret_tiny;
1951 }
1952#endif
1953 goto retz;
1954#ifdef IEEE_Arith
1955 ret_tinyf:
1956 Bfree(b);
1957 ret_tiny:
1958#ifndef NO_ERRNO
1959 errno = ERANGE;
1960#endif
1961 word0(rvp) = 0;
1962 word1(rvp) = 1;
1963 return;
1964#endif /* IEEE_Arith */
1965 }
1966 switch(rounding) {
1967 case Round_near:
1968 goto ovfl1;
1969 case Round_up:
1970 if (!sign)
1971 goto ovfl1;
1972 goto ret_big;
1973 case Round_down:
1974 if (sign)
1975 goto ovfl1;
1976 goto ret_big;
1977 }
1978 ret_big:
1979 word0(rvp) = Big0;
1980 word1(rvp) = Big1;
1981 return;
1982 }
1983 n = s1 - s0 - 1;
1984 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
1985 k++;
1986 b = Balloc(k);
1987 x = b->x;
1988 n = 0;
1989 L = 0;
1990#ifdef USE_LOCALE
1991 for(i = 0; decimalpoint[i+1]; ++i);
1992#endif
1993 while(s1 > s0) {
1994#ifdef USE_LOCALE
1995 if (*--s1 == decimalpoint[i]) {
1996 s1 -= i;
1997 continue;
1998 }
1999#else
2000 if (*--s1 == '.')
2001 continue;
2002#endif
2003 if (n == ULbits) {
2004 *x++ = L;
2005 L = 0;
2006 n = 0;
2007 }
2008 L |= (hexdig[*s1] & 0x0f) << n;
2009 n += 4;
2010 }
2011 *x++ = L;
2012 b->wds = n = x - b->x;
2013 n = ULbits*n - hi0bits(L);
2014 nbits = Nbits;
2015 lostbits = 0;
2016 x = b->x;
2017 if (n > nbits) {
2018 n -= nbits;
2019 if (any_on(b,n)) {
2020 lostbits = 1;
2021 k = n - 1;
2022 if (x[k>>kshift] & 1 << (k & kmask)) {
2023 lostbits = 2;
2024 if (k > 0 && any_on(b,k))
2025 lostbits = 3;
2026 }
2027 }
2028 rshift(b, n);
2029 e += n;
2030 }
2031 else if (n < nbits) {
2032 n = nbits - n;
2033 b = lshift(b, n);
2034 e -= n;
2035 x = b->x;
2036 }
2037 if (e > Emax) {
2038 ovfl:
2039 Bfree(b);
2040 ovfl1:
2041#ifndef NO_ERRNO
2042 errno = ERANGE;
2043#endif
2044 word0(rvp) = Exp_mask;
2045 word1(rvp) = 0;
2046 return;
2047 }
2048 denorm = 0;
2049 if (e < emin) {
2050 denorm = 1;
2051 n = emin - e;
2052 if (n >= nbits) {
2053#ifdef IEEE_Arith /*{*/
2054 switch (rounding) {
2055 case Round_near:
2056 if (n == nbits && (n < 2 || any_on(b,n-1)))
2057 goto ret_tinyf;
2058 break;
2059 case Round_up:
2060 if (!sign)
2061 goto ret_tinyf;
2062 break;
2063 case Round_down:
2064 if (sign)
2065 goto ret_tinyf;
2066 }
2067#endif /* } IEEE_Arith */
2068 Bfree(b);
2069 retz:
2070#ifndef NO_ERRNO
2071 errno = ERANGE;
2072#endif
2073 retz1:
2074 rvp->d = 0.;
2075 return;
2076 }
2077 k = n - 1;
2078 if (lostbits)
2079 lostbits = 1;
2080 else if (k > 0)
2081 lostbits = any_on(b,k);
2082 if (x[k>>kshift] & 1 << (k & kmask))
2083 lostbits |= 2;
2084 nbits -= n;
2085 rshift(b,n);
2086 e = emin;
2087 }
2088 if (lostbits) {
2089 up = 0;
2090 switch(rounding) {
2091 case Round_zero:
2092 break;
2093 case Round_near:
2094 if (lostbits & 2
2095 && (lostbits & 1) | (x[0] & 1))
2096 up = 1;
2097 break;
2098 case Round_up:
2099 up = 1 - sign;
2100 break;
2101 case Round_down:
2102 up = sign;
2103 }
2104 if (up) {
2105 k = b->wds;
2106 b = increment(b);
2107 x = b->x;
2108 if (denorm) {
2109#if 0
2110 if (nbits == Nbits - 1
2111 && x[nbits >> kshift] & 1 << (nbits & kmask))
2112 denorm = 0; /* not currently used */
2113#endif
2114 }
2115 else if (b->wds > k
2116 || ((n = nbits & kmask) !=0
2117 && hi0bits(x[k-1]) < 32-n)) {
2118 rshift(b,1);
2119 if (++e > Emax)
2120 goto ovfl;
2121 }
2122 }
2123 }
2124#ifdef IEEE_Arith
2125 if (denorm)
2126 word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
2127 else
2128 word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2129 word1(rvp) = b->x[0];
2130#endif
2131#ifdef IBM
2132 if ((j = e & 3)) {
2133 k = b->x[0] & ((1 << j) - 1);
2134 rshift(b,j);
2135 if (k) {
2136 switch(rounding) {
2137 case Round_up:
2138 if (!sign)
2139 increment(b);
2140 break;
2141 case Round_down:
2142 if (sign)
2143 increment(b);
2144 break;
2145 case Round_near:
2146 j = 1 << (j-1);
2147 if (k & j && ((k & (j-1)) | lostbits))
2148 increment(b);
2149 }
2150 }
2151 }
2152 e >>= 2;
2153 word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
2154 word1(rvp) = b->x[0];
2155#endif
2156#ifdef VAX
2157 /* The next two lines ignore swap of low- and high-order 2 bytes. */
2158 /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2159 /* word1(rvp) = b->x[0]; */
2160 word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
2161 word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
2162#endif
2163 Bfree(b);
2164 }
2165#endif
2167 static int
2168#ifdef KR_headers
2169dshift(b, p2) Bigint *b; int p2;
2170#else
2171dshift(Bigint *b, int p2)
2172#endif
2173{
2174 int rv = hi0bits(b->x[b->wds-1]) - 4;
2175 if (p2 > 0)
2176 rv -= p2;
2177 return rv & kmask;
2178 }
2179
2180 static int
2182#ifdef KR_headers
2183 (b, S) Bigint *b, *S;
2184#else
2185 (Bigint *b, Bigint *S)
2186#endif
2187{
2188 int n;
2189 ULong *bx, *bxe, q, *sx, *sxe;
2190#ifdef ULLong
2191 ULLong borrow, carry, y, ys;
2192#else
2193 ULong borrow, carry, y, ys;
2194#ifdef Pack_32
2195 ULong si, z, zs;
2196#endif
2197#endif
2198
2199 n = S->wds;
2200#ifdef DEBUG
2201 /*debug*/ if (b->wds > n)
2202 /*debug*/ Bug("oversize b in quorem");
2203#endif
2204 if (b->wds < n)
2205 return 0;
2206 sx = S->x;
2207 sxe = sx + --n;
2208 bx = b->x;
2209 bxe = bx + n;
2210 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2211#ifdef DEBUG
2212#ifdef NO_STRTOD_BIGCOMP
2213 /*debug*/ if (q > 9)
2214#else
2215 /* An oversized q is possible when quorem is called from bigcomp and */
2216 /* the input is near, e.g., twice the smallest denormalized number. */
2217 /*debug*/ if (q > 15)
2218#endif
2219 /*debug*/ Bug("oversized quotient in quorem");
2220#endif
2221 if (q) {
2222 borrow = 0;
2223 carry = 0;
2224 do {
2225#ifdef ULLong
2226 ys = *sx++ * (ULLong)q + carry;
2227 carry = ys >> 32;
2228 y = *bx - (ys & FFFFFFFF) - borrow;
2229 borrow = y >> 32 & (ULong)1;
2230 *bx++ = y & FFFFFFFF;
2231#else
2232#ifdef Pack_32
2233 si = *sx++;
2234 ys = (si & 0xffff) * q + carry;
2235 zs = (si >> 16) * q + (ys >> 16);
2236 carry = zs >> 16;
2237 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2238 borrow = (y & 0x10000) >> 16;
2239 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2240 borrow = (z & 0x10000) >> 16;
2241 Storeinc(bx, z, y);
2242#else
2243 ys = *sx++ * q + carry;
2244 carry = ys >> 16;
2245 y = *bx - (ys & 0xffff) - borrow;
2246 borrow = (y & 0x10000) >> 16;
2247 *bx++ = y & 0xffff;
2248#endif
2249#endif
2250 }
2251 while(sx <= sxe);
2252 if (!*bxe) {
2253 bx = b->x;
2254 while(--bxe > bx && !*bxe)
2255 --n;
2256 b->wds = n;
2257 }
2258 }
2259 if (cmp(b, S) >= 0) {
2260 q++;
2261 borrow = 0;
2262 carry = 0;
2263 bx = b->x;
2264 sx = S->x;
2265 do {
2266#ifdef ULLong
2267 ys = *sx++ + carry;
2268 carry = ys >> 32;
2269 y = *bx - (ys & FFFFFFFF) - borrow;
2270 borrow = y >> 32 & (ULong)1;
2271 *bx++ = y & FFFFFFFF;
2272#else
2273#ifdef Pack_32
2274 si = *sx++;
2275 ys = (si & 0xffff) + carry;
2276 zs = (si >> 16) + (ys >> 16);
2277 carry = zs >> 16;
2278 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2279 borrow = (y & 0x10000) >> 16;
2280 z = (*bx >> 16) - (zs & 0xffff) - borrow;
2281 borrow = (z & 0x10000) >> 16;
2282 Storeinc(bx, z, y);
2283#else
2284 ys = *sx++ + carry;
2285 carry = ys >> 16;
2286 y = *bx - (ys & 0xffff) - borrow;
2287 borrow = (y & 0x10000) >> 16;
2288 *bx++ = y & 0xffff;
2289#endif
2290#endif
2291 }
2292 while(sx <= sxe);
2293 bx = b->x;
2294 bxe = bx + n;
2295 if (!*bxe) {
2296 while(--bxe > bx && !*bxe)
2297 --n;
2298 b->wds = n;
2299 }
2300 }
2301 return q;
2302 }
2303
2304#if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
2305 static double
2307#ifdef KR_headers
2308 (x, bc) U *x; BCinfo *bc;
2309#else
2310 (U *x, BCinfo *bc)
2311#endif
2312{
2313 U u;
2314 double rv;
2315 int i;
2316
2317 rv = ulp(x);
2318 if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
2319 return rv; /* Is there an example where i <= 0 ? */
2320 word0(&u) = Exp_1 + (i << Exp_shift);
2321 word1(&u) = 0;
2322 return rv * u.d;
2323 }
2324#endif /*}*/
2325
2326#ifndef NO_STRTOD_BIGCOMP
2327 static void
2329#ifdef KR_headers
2330 (rv, s0, bc)
2331 U *rv; CONST char *s0; BCinfo *bc;
2332#else
2333 (U *rv, const char *s0, BCinfo *bc)
2334#endif
2335{
2336 Bigint *b, *d;
2337 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2338
2339 dsign = bc->dsign;
2340 nd = bc->nd;
2341 nd0 = bc->nd0;
2342 p5 = nd + bc->e0 - 1;
2343 speccase = 0;
2344#ifndef Sudden_Underflow
2345 if (rv->d == 0.) { /* special case: value near underflow-to-zero */
2346 /* threshold was rounded to zero */
2347 b = i2b(1);
2348 p2 = Emin - P + 1;
2349 bbits = 1;
2350#ifdef Avoid_Underflow
2351 word0(rv) = (P+2) << Exp_shift;
2352#else
2353 word1(rv) = 1;
2354#endif
2355 i = 0;
2356#ifdef Honor_FLT_ROUNDS
2357 if (bc->rounding == 1)
2358#endif
2359 {
2360 speccase = 1;
2361 --p2;
2362 dsign = 0;
2363 goto have_i;
2364 }
2365 }
2366 else
2367#endif
2368 b = d2b(rv, &p2, &bbits);
2369#ifdef Avoid_Underflow
2370 p2 -= bc->scale;
2371#endif
2372 /* floor(log2(rv)) == bbits - 1 + p2 */
2373 /* Check for denormal case. */
2374 i = P - bbits;
2375 if (i > (j = P - Emin - 1 + p2)) {
2376#ifdef Sudden_Underflow
2377 Bfree(b);
2378 b = i2b(1);
2379 p2 = Emin;
2380 i = P - 1;
2381#ifdef Avoid_Underflow
2382 word0(rv) = (1 + bc->scale) << Exp_shift;
2383#else
2384 word0(rv) = Exp_msk1;
2385#endif
2386 word1(rv) = 0;
2387#else
2388 i = j;
2389#endif
2390 }
2391#ifdef Honor_FLT_ROUNDS
2392 if (bc->rounding != 1) {
2393 if (i > 0)
2394 b = lshift(b, i);
2395 if (dsign)
2396 b = increment(b);
2397 }
2398 else
2399#endif
2400 {
2401 b = lshift(b, ++i);
2402 b->x[0] |= 1;
2403 }
2404#ifndef Sudden_Underflow
2405 have_i:
2406#endif
2407 p2 -= p5 + i;
2408 d = i2b(1);
2409 /* Arrange for convenient computation of quotients:
2410 * shift left if necessary so divisor has 4 leading 0 bits.
2411 */
2412 if (p5 > 0)
2413 d = pow5mult(d, p5);
2414 else if (p5 < 0)
2415 b = pow5mult(b, -p5);
2416 if (p2 > 0) {
2417 b2 = p2;
2418 d2 = 0;
2419 }
2420 else {
2421 b2 = 0;
2422 d2 = -p2;
2423 }
2424 i = dshift(d, d2);
2425 if ((b2 += i) > 0)
2426 b = lshift(b, b2);
2427 if ((d2 += i) > 0)
2428 d = lshift(d, d2);
2429
2430 /* Now b/d = exactly half-way between the two floating-point values */
2431 /* on either side of the input string. Compute first digit of b/d. */
2432
2433 if (!(dig = quorem(b,d))) {
2434 b = multadd(b, 10, 0); /* very unlikely */
2435 dig = quorem(b,d);
2436 }
2437
2438 /* Compare b/d with s0 */
2439
2440 for(i = 0; i < nd0; ) {
2441 if ((dd = s0[i++] - '0' - dig))
2442 goto ret;
2443 if (!b->x[0] && b->wds == 1) {
2444 if (i < nd)
2445 dd = 1;
2446 goto ret;
2447 }
2448 b = multadd(b, 10, 0);
2449 dig = quorem(b,d);
2450 }
2451 for(j = bc->dp1; i++ < nd;) {
2452 if ((dd = s0[j++] - '0' - dig))
2453 goto ret;
2454 if (!b->x[0] && b->wds == 1) {
2455 if (i < nd)
2456 dd = 1;
2457 goto ret;
2458 }
2459 b = multadd(b, 10, 0);
2460 dig = quorem(b,d);
2461 }
2462 if (dig > 0 || b->x[0] || b->wds > 1)
2463 dd = -1;
2464 ret:
2465 Bfree(b);
2466 Bfree(d);
2467#ifdef Honor_FLT_ROUNDS
2468 if (bc->rounding != 1) {
2469 if (dd < 0) {
2470 if (bc->rounding == 0) {
2471 if (!dsign)
2472 goto retlow1;
2473 }
2474 else if (dsign)
2475 goto rethi1;
2476 }
2477 else if (dd > 0) {
2478 if (bc->rounding == 0) {
2479 if (dsign)
2480 goto rethi1;
2481 goto ret1;
2482 }
2483 if (!dsign)
2484 goto rethi1;
2485 dval(rv) += 2.*sulp(rv,bc);
2486 }
2487 else {
2488 bc->inexact = 0;
2489 if (dsign)
2490 goto rethi1;
2491 }
2492 }
2493 else
2494#endif
2495 if (speccase) {
2496 if (dd <= 0)
2497 rv->d = 0.;
2498 }
2499 else if (dd < 0) {
2500 if (!dsign) /* does not happen for round-near */
2501retlow1:
2502 dval(rv) -= sulp(rv,bc);
2503 }
2504 else if (dd > 0) {
2505 if (dsign) {
2506 rethi1:
2507 dval(rv) += sulp(rv,bc);
2508 }
2509 }
2510 else {
2511 /* Exact half-way case: apply round-even rule. */
2512 if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
2513 i = 1 - j;
2514 if (i <= 31) {
2515 if (word1(rv) & (0x1 << i))
2516 goto odd;
2517 }
2518 else if (word0(rv) & (0x1 << (i-32)))
2519 goto odd;
2520 }
2521 else if (word1(rv) & 1) {
2522 odd:
2523 if (dsign)
2524 goto rethi1;
2525 goto retlow1;
2526 }
2527 }
2528
2529#ifdef Honor_FLT_ROUNDS
2530 ret1:
2531#endif
2532 return;
2533 }
2534#endif /* NO_STRTOD_BIGCOMP */
2535
2536 double
2538#ifdef KR_headers
2539 (s00, se) CONST char *s00; char **se;
2540#else
2541 (const char *s00, char **se)
2542#endif
2543{
2544 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2545 int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
2546 CONST char *s, *s0, *s1;
2547 double aadj, aadj1;
2548 Long L;
2549 U aadj2, adj, rv, rv0;
2550 ULong y, z;
2551 BCinfo bc;
2552 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2553#ifdef Avoid_Underflow
2554 ULong Lsb, Lsb1;
2555#endif
2556#ifdef SET_INEXACT
2557 int oldinexact;
2558#endif
2559#ifndef NO_STRTOD_BIGCOMP
2560 int req_bigcomp = 0;
2561#endif
2562#ifdef Honor_FLT_ROUNDS /*{*/
2563#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2564 bc.rounding = Flt_Rounds;
2565#else /*}{*/
2566 bc.rounding = 1;
2567 switch(fegetround()) {
2568 case FE_TOWARDZERO: bc.rounding = 0; break;
2569 case FE_UPWARD: bc.rounding = 2; break;
2570 case FE_DOWNWARD: bc.rounding = 3;
2571 }
2572#endif /*}}*/
2573#endif /*}*/
2574#ifdef USE_LOCALE
2575 CONST char *s2;
2576#endif
2577
2578 sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
2579 dval(&rv) = 0.;
2580 for(s = s00;;s++) switch(*s) {
2581 case '-':
2582 sign = 1;
2583 /* no break */
2584 case '+':
2585 if (*++s)
2586 goto break2;
2587 /* no break */
2588 case 0:
2589 goto ret0;
2590 case '\t':
2591 case '\n':
2592 case '\v':
2593 case '\f':
2594 case '\r':
2595 case ' ':
2596 continue;
2597 default:
2598 goto break2;
2599 }
2600 break2:
2601 if (*s == '0') {
2602#ifndef NO_HEX_FP /*{*/
2603 switch(s[1]) {
2604 case 'x':
2605 case 'X':
2606#ifdef Honor_FLT_ROUNDS
2607 os_gethex(&s, &rv, bc.rounding, sign);
2608#else
2609 os_gethex(&s, &rv, 1, sign);
2610#endif
2611 goto ret;
2612 }
2613#endif /*}*/
2614 nz0 = 1;
2615 while(*++s == '0') ;
2616 if (!*s)
2617 goto ret;
2618 }
2619 s0 = s;
2620 y = z = 0;
2621 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2622 if (nd < 9)
2623 y = 10*y + c - '0';
2624 else if (nd < 16)
2625 z = 10*z + c - '0';
2626 nd0 = nd;
2627 bc.dp0 = bc.dp1 = s - s0;
2628 for(s1 = s; s1 > s0 && *--s1 == '0'; )
2629 ++nz1;
2630#ifdef USE_LOCALE
2631 s1 = localeconv()->decimal_point;
2632 if (c == *s1) {
2633 c = '.';
2634 if (*++s1) {
2635 s2 = s;
2636 for(;;) {
2637 if (*++s2 != *s1) {
2638 c = 0;
2639 break;
2640 }
2641 if (!*++s1) {
2642 s = s2;
2643 break;
2644 }
2645 }
2646 }
2647 }
2648#endif
2649 if (c == '.') {
2650 c = *++s;
2651 bc.dp1 = s - s0;
2652 bc.dplen = bc.dp1 - bc.dp0;
2653 if (!nd) {
2654 for(; c == '0'; c = *++s)
2655 nz++;
2656 if (c > '0' && c <= '9') {
2657 bc.dp0 = s0 - s;
2658 bc.dp1 = bc.dp0 + bc.dplen;
2659 s0 = s;
2660 nf += nz;
2661 nz = 0;
2662 goto have_dig;
2663 }
2664 goto dig_done;
2665 }
2666 for(; c >= '0' && c <= '9'; c = *++s) {
2667 have_dig:
2668 nz++;
2669 if (c -= '0') {
2670 nf += nz;
2671 for(i = 1; i < nz; i++)
2672 if (nd++ < 9)
2673 y *= 10;
2674 else if (nd <= DBL_DIG + 1)
2675 z *= 10;
2676 if (nd++ < 9)
2677 y = 10*y + c;
2678 else if (nd <= DBL_DIG + 1)
2679 z = 10*z + c;
2680 nz = nz1 = 0;
2681 }
2682 }
2683 }
2684 dig_done:
2685 e = 0;
2686 if (c == 'e' || c == 'E') {
2687 if (!nd && !nz && !nz0) {
2688 goto ret0;
2689 }
2690 s00 = s;
2691 esign = 0;
2692 switch(c = *++s) {
2693 case '-':
2694 esign = 1;
2695 case '+':
2696 c = *++s;
2697 }
2698 if (c >= '0' && c <= '9') {
2699 while(c == '0')
2700 c = *++s;
2701 if (c > '0' && c <= '9') {
2702 L = c - '0';
2703 s1 = s;
2704 while((c = *++s) >= '0' && c <= '9')
2705 L = 10*L + c - '0';
2706 if (s - s1 > 8 || L > 19999)
2707 /* Avoid confusion from exponents
2708 * so large that e might overflow.
2709 */
2710 e = 19999; /* safe for 16 bit ints */
2711 else
2712 e = (int)L;
2713 if (esign)
2714 e = -e;
2715 }
2716 else
2717 e = 0;
2718 }
2719 else
2720 s = s00;
2721 }
2722 if (!nd) {
2723 if (!nz && !nz0) {
2724#ifdef INFNAN_CHECK
2725 /* Check for Nan and Infinity */
2726 if (!bc.dplen)
2727 switch(c) {
2728 case 'i':
2729 case 'I':
2730 if (match(&s,"nf")) {
2731 --s;
2732 if (!match(&s,"inity"))
2733 ++s;
2734 word0(&rv) = 0x7ff00000;
2735 word1(&rv) = 0;
2736 goto ret;
2737 }
2738 break;
2739 case 'n':
2740 case 'N':
2741 if (match(&s, "an")) {
2742 word0(&rv) = NAN_WORD0;
2743 word1(&rv) = NAN_WORD1;
2744#ifndef No_Hex_NaN
2745 if (*s == '(') /*)*/
2746 hexnan(&rv, &s);
2747#endif
2748 goto ret;
2749 }
2750 }
2751#endif /* INFNAN_CHECK */
2752 ret0:
2753 s = s00;
2754 sign = 0;
2755 }
2756 goto ret;
2757 }
2758 bc.e0 = e1 = e -= nf;
2759
2760 /* Now we have nd0 digits, starting at s0, followed by a
2761 * decimal point, followed by nd-nd0 digits. The number we're
2762 * after is the integer represented by those digits times
2763 * 10**e */
2764
2765 if (!nd0)
2766 nd0 = nd;
2767 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2768 dval(&rv) = y;
2769 if (k > 9) {
2770#ifdef SET_INEXACT
2771 if (k > DBL_DIG)
2772 oldinexact = get_inexact();
2773#endif
2774 dval(&rv) = tens[k - 9] * dval(&rv) + z;
2775 }
2776 bd0 = 0;
2777 if (nd <= DBL_DIG
2778#ifndef RND_PRODQUOT
2779#ifndef Honor_FLT_ROUNDS
2780 && Flt_Rounds == 1
2781#endif
2782#endif
2783 ) {
2784 if (!e)
2785 goto ret;
2786#ifndef ROUND_BIASED_without_Round_Up
2787 if (e > 0) {
2788 if (e <= Ten_pmax) {
2789#ifdef VAX
2790 goto vax_ovfl_check;
2791#else
2792#ifdef Honor_FLT_ROUNDS
2793 /* round correctly FLT_ROUNDS = 2 or 3 */
2794 if (sign) {
2795 rv.d = -rv.d;
2796 sign = 0;
2797 }
2798#endif
2799 /* rv = */ rounded_product(dval(&rv), tens[e]);
2800 goto ret;
2801#endif
2802 }
2803 i = DBL_DIG - nd;
2804 if (e <= Ten_pmax + i) {
2805 /* A fancier test would sometimes let us do
2806 * this for larger i values.
2807 */
2808#ifdef Honor_FLT_ROUNDS
2809 /* round correctly FLT_ROUNDS = 2 or 3 */
2810 if (sign) {
2811 rv.d = -rv.d;
2812 sign = 0;
2813 }
2814#endif
2815 e -= i;
2816 dval(&rv) *= tens[i];
2817#ifdef VAX
2818 /* VAX exponent range is so narrow we must
2819 * worry about overflow here...
2820 */
2821 vax_ovfl_check:
2822 word0(&rv) -= P*Exp_msk1;
2823 /* rv = */ rounded_product(dval(&rv), tens[e]);
2824 if ((word0(&rv) & Exp_mask)
2825 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2826 goto ovfl;
2827 word0(&rv) += P*Exp_msk1;
2828#else
2829 /* rv = */ rounded_product(dval(&rv), tens[e]);
2830#endif
2831 goto ret;
2832 }
2833 }
2834#ifndef Inaccurate_Divide
2835 else if (e >= -Ten_pmax) {
2836#ifdef Honor_FLT_ROUNDS
2837 /* round correctly FLT_ROUNDS = 2 or 3 */
2838 if (sign) {
2839 rv.d = -rv.d;
2840 sign = 0;
2841 }
2842#endif
2843 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2844 goto ret;
2845 }
2846#endif
2847#endif /* ROUND_BIASED_without_Round_Up */
2848 }
2849 e1 += nd - k;
2850
2851#ifdef IEEE_Arith
2852#ifdef SET_INEXACT
2853 bc.inexact = 1;
2854 if (k <= DBL_DIG)
2855 oldinexact = get_inexact();
2856#endif
2857#ifdef Avoid_Underflow
2858 bc.scale = 0;
2859#endif
2860#ifdef Honor_FLT_ROUNDS
2861 if (bc.rounding >= 2) {
2862 if (sign)
2863 bc.rounding = bc.rounding == 2 ? 0 : 2;
2864 else
2865 if (bc.rounding != 2)
2866 bc.rounding = 0;
2867 }
2868#endif
2869#endif /*IEEE_Arith*/
2870
2871 /* Get starting approximation = rv * 10**e1 */
2872
2873 if (e1 > 0) {
2874 if ((i = e1 & 15))
2875 dval(&rv) *= tens[i];
2876 if (e1 &= ~15) {
2877 if (e1 > DBL_MAX_10_EXP) {
2878 ovfl:
2879 /* Can't trust HUGE_VAL */
2880#ifdef IEEE_Arith
2881#ifdef Honor_FLT_ROUNDS
2882 switch(bc.rounding) {
2883 case 0: /* toward 0 */
2884 case 3: /* toward -infinity */
2885 word0(&rv) = Big0;
2886 word1(&rv) = Big1;
2887 break;
2888 default:
2889 word0(&rv) = Exp_mask;
2890 word1(&rv) = 0;
2891 }
2892#else /*Honor_FLT_ROUNDS*/
2893 word0(&rv) = Exp_mask;
2894 word1(&rv) = 0;
2895#endif /*Honor_FLT_ROUNDS*/
2896#ifdef SET_INEXACT
2897 /* set overflow bit */
2898 dval(&rv0) = 1e300;
2899 dval(&rv0) *= dval(&rv0);
2900#endif
2901#else /*IEEE_Arith*/
2902 word0(&rv) = Big0;
2903 word1(&rv) = Big1;
2904#endif /*IEEE_Arith*/
2905 range_err:
2906 if (bd0) {
2907 Bfree(bb);
2908 Bfree(bd);
2909 Bfree(bs);
2910 Bfree(bd0);
2911 Bfree(delta);
2912 }
2913#ifndef NO_ERRNO
2914 errno = ERANGE;
2915#endif
2916 goto ret;
2917 }
2918 e1 >>= 4;
2919 for(j = 0; e1 > 1; j++, e1 >>= 1)
2920 if (e1 & 1)
2921 dval(&rv) *= bigtens[j];
2922 /* The last multiplication could overflow. */
2923 word0(&rv) -= P*Exp_msk1;
2924 dval(&rv) *= bigtens[j];
2925 if ((z = word0(&rv) & Exp_mask)
2926 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2927 goto ovfl;
2928 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2929 /* set to largest number */
2930 /* (Can't trust DBL_MAX) */
2931 word0(&rv) = Big0;
2932 word1(&rv) = Big1;
2933 }
2934 else
2935 word0(&rv) += P*Exp_msk1;
2936 }
2937 }
2938 else if (e1 < 0) {
2939 e1 = -e1;
2940 if ((i = e1 & 15))
2941 dval(&rv) /= tens[i];
2942 if (e1 >>= 4) {
2943 if (e1 >= 1 << n_bigtens)
2944 goto undfl;
2945#ifdef Avoid_Underflow
2946 if (e1 & Scale_Bit)
2947 bc.scale = 2*P;
2948 for(j = 0; e1 > 0; j++, e1 >>= 1)
2949 if (e1 & 1)
2950 dval(&rv) *= tinytens[j];
2951 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2952 >> Exp_shift)) > 0) {
2953 /* scaled rv is denormal; clear j low bits */
2954 if (j >= 32) {
2955 if (j > 54)
2956 goto undfl;
2957 word1(&rv) = 0;
2958 if (j >= 53)
2959 word0(&rv) = (P+2)*Exp_msk1;
2960 else
2961 word0(&rv) &= 0xffffffff << (j-32);
2962 }
2963 else
2964 word1(&rv) &= 0xffffffff << j;
2965 }
2966#else
2967 for(j = 0; e1 > 1; j++, e1 >>= 1)
2968 if (e1 & 1)
2969 dval(&rv) *= tinytens[j];
2970 /* The last multiplication could underflow. */
2971 dval(&rv0) = dval(&rv);
2972 dval(&rv) *= tinytens[j];
2973 if (!dval(&rv)) {
2974 dval(&rv) = 2.*dval(&rv0);
2975 dval(&rv) *= tinytens[j];
2976#endif
2977 if (!dval(&rv)) {
2978 undfl:
2979 dval(&rv) = 0.;
2980 goto range_err;
2981 }
2982#ifndef Avoid_Underflow
2983 word0(&rv) = Tiny0;
2984 word1(&rv) = Tiny1;
2985 /* The refinement below will clean
2986 * this approximation up.
2987 */
2988 }
2989#endif
2990 }
2991 }
2992
2993 /* Now the hard part -- adjusting rv to the correct value.*/
2994
2995 /* Put digits into bd: true value = bd * 10^e */
2996
2997 bc.nd = nd - nz1;
2998#ifndef NO_STRTOD_BIGCOMP
2999 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
3000 /* to silence an erroneous warning about bc.nd0 */
3001 /* possibly not being initialized. */
3002 if (nd > strtod_diglim) {
3003 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
3004 /* minimum number of decimal digits to distinguish double values */
3005 /* in IEEE arithmetic. */
3006 i = j = 18;
3007 if (i > nd0)
3008 j += bc.dplen;
3009 for(;;) {
3010 if (--j < bc.dp1 && j >= bc.dp0)
3011 j = bc.dp0 - 1;
3012 if (s0[j] != '0')
3013 break;
3014 --i;
3015 }
3016 e += nd - i;
3017 nd = i;
3018 if (nd0 > nd)
3019 nd0 = nd;
3020 if (nd < 9) { /* must recompute y */
3021 y = 0;
3022 for(i = 0; i < nd0; ++i)
3023 y = 10*y + s0[i] - '0';
3024 for(j = bc.dp1; i < nd; ++i)
3025 y = 10*y + s0[j++] - '0';
3026 }
3027 }
3028#endif
3029 bd0 = s2b(s0, nd0, nd, y, bc.dplen);
3030
3031 for(;;) {
3032 bd = Balloc(bd0->k);
3033 Bcopy(bd, bd0);
3034 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
3035 bs = i2b(1);
3036
3037 if (e >= 0) {
3038 bb2 = bb5 = 0;
3039 bd2 = bd5 = e;
3040 }
3041 else {
3042 bb2 = bb5 = -e;
3043 bd2 = bd5 = 0;
3044 }
3045 if (bbe >= 0)
3046 bb2 += bbe;
3047 else
3048 bd2 -= bbe;
3049 bs2 = bb2;
3050#ifdef Honor_FLT_ROUNDS
3051 if (bc.rounding != 1)
3052 bs2++;
3053#endif
3054#ifdef Avoid_Underflow
3055 Lsb = LSB;
3056 Lsb1 = 0;
3057 j = bbe - bc.scale;
3058 i = j + bbbits - 1; /* logb(rv) */
3059 j = P + 1 - bbbits;
3060 if (i < Emin) { /* denormal */
3061 i = Emin - i;
3062 j -= i;
3063 if (i < 32)
3064 Lsb <<= i;
3065 else if (i < 52)
3066 Lsb1 = Lsb << (i-32);
3067 else
3068 Lsb1 = Exp_mask;
3069 }
3070#else /*Avoid_Underflow*/
3071#ifdef Sudden_Underflow
3072#ifdef IBM
3073 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
3074#else
3075 j = P + 1 - bbbits;
3076#endif
3077#else /*Sudden_Underflow*/
3078 j = bbe;
3079 i = j + bbbits - 1; /* logb(rv) */
3080 if (i < Emin) /* denormal */
3081 j += P - Emin;
3082 else
3083 j = P + 1 - bbbits;
3084#endif /*Sudden_Underflow*/
3085#endif /*Avoid_Underflow*/
3086 bb2 += j;
3087 bd2 += j;
3088#ifdef Avoid_Underflow
3089 bd2 += bc.scale;
3090#endif
3091 i = bb2 < bd2 ? bb2 : bd2;
3092 if (i > bs2)
3093 i = bs2;
3094 if (i > 0) {
3095 bb2 -= i;
3096 bd2 -= i;
3097 bs2 -= i;
3098 }
3099 if (bb5 > 0) {
3100 bs = pow5mult(bs, bb5);
3101 bb1 = mult(bs, bb);
3102 Bfree(bb);
3103 bb = bb1;
3104 }
3105 if (bb2 > 0)
3106 bb = lshift(bb, bb2);
3107 if (bd5 > 0)
3108 bd = pow5mult(bd, bd5);
3109 if (bd2 > 0)
3110 bd = lshift(bd, bd2);
3111 if (bs2 > 0)
3112 bs = lshift(bs, bs2);
3113 delta = diff(bb, bd);
3114 bc.dsign = delta->sign;
3115 delta->sign = 0;
3116 i = cmp(delta, bs);
3117#ifndef NO_STRTOD_BIGCOMP /*{*/
3118 if (bc.nd > nd && i <= 0) {
3119 if (bc.dsign) {
3120 /* Must use bigcomp(). */
3121 req_bigcomp = 1;
3122 break;
3123 }
3124#ifdef Honor_FLT_ROUNDS
3125 if (bc.rounding != 1) {
3126 if (i < 0) {
3127 req_bigcomp = 1;
3128 break;
3129 }
3130 }
3131 else
3132#endif
3133 i = -1; /* Discarded digits make delta smaller. */
3134 }
3135#endif /*}*/
3136#ifdef Honor_FLT_ROUNDS /*{*/
3137 if (bc.rounding != 1) {
3138 if (i < 0) {
3139 /* Error is less than an ulp */
3140 if (!delta->x[0] && delta->wds <= 1) {
3141 /* exact */
3142#ifdef SET_INEXACT
3143 bc.inexact = 0;
3144#endif
3145 break;
3146 }
3147 if (bc.rounding) {
3148 if (bc.dsign) {
3149 adj.d = 1.;
3150 goto apply_adj;
3151 }
3152 }
3153 else if (!bc.dsign) {
3154 adj.d = -1.;
3155 if (!word1(&rv)
3156 && !(word0(&rv) & Frac_mask)) {
3157 y = word0(&rv) & Exp_mask;
3158#ifdef Avoid_Underflow
3159 if (!bc.scale || y > 2*P*Exp_msk1)
3160#else
3161 if (y)
3162#endif
3163 {
3164 delta = lshift(delta,Log2P);
3165 if (cmp(delta, bs) <= 0)
3166 adj.d = -0.5;
3167 }
3168 }
3169 apply_adj:
3170#ifdef Avoid_Underflow /*{*/
3171 if (bc.scale && (y = word0(&rv) & Exp_mask)
3172 <= 2*P*Exp_msk1)
3173 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3174#else
3175#ifdef Sudden_Underflow
3176 if ((word0(&rv) & Exp_mask) <=
3177 P*Exp_msk1) {
3178 word0(&rv) += P*Exp_msk1;
3179 dval(&rv) += adj.d*ulp(dval(&rv));
3180 word0(&rv) -= P*Exp_msk1;
3181 }
3182 else
3183#endif /*Sudden_Underflow*/
3184#endif /*Avoid_Underflow}*/
3185 dval(&rv) += adj.d*ulp(&rv);
3186 }
3187 break;
3188 }
3189 adj.d = ratio(delta, bs);
3190 if (adj.d < 1.)
3191 adj.d = 1.;
3192 if (adj.d <= 0x7ffffffe) {
3193 /* adj = rounding ? ceil(adj) : floor(adj); */
3194 y = adj.d;
3195 if (y != adj.d) {
3196 if (!((bc.rounding>>1) ^ bc.dsign))
3197 y++;
3198 adj.d = y;
3199 }
3200 }
3201#ifdef Avoid_Underflow /*{*/
3202 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3203 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3204#else
3205#ifdef Sudden_Underflow
3206 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3207 word0(&rv) += P*Exp_msk1;
3208 adj.d *= ulp(dval(&rv));
3209 if (bc.dsign)
3210 dval(&rv) += adj.d;
3211 else
3212 dval(&rv) -= adj.d;
3213 word0(&rv) -= P*Exp_msk1;
3214 goto cont;
3215 }
3216#endif /*Sudden_Underflow*/
3217#endif /*Avoid_Underflow}*/
3218 adj.d *= ulp(&rv);
3219 if (bc.dsign) {
3220 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3221 goto ovfl;
3222 dval(&rv) += adj.d;
3223 }
3224 else
3225 dval(&rv) -= adj.d;
3226 goto cont;
3227 }
3228#endif /*}Honor_FLT_ROUNDS*/
3229
3230 if (i < 0) {
3231 /* Error is less than half an ulp -- check for
3232 * special case of mantissa a power of two.
3233 */
3234 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3235#ifdef IEEE_Arith /*{*/
3236#ifdef Avoid_Underflow
3237 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3238#else
3239 || (word0(&rv) & Exp_mask) <= Exp_msk1
3240#endif
3241#endif /*}*/
3242 ) {
3243#ifdef SET_INEXACT
3244 if (!delta->x[0] && delta->wds <= 1)
3245 bc.inexact = 0;
3246#endif
3247 break;
3248 }
3249 if (!delta->x[0] && delta->wds <= 1) {
3250 /* exact result */
3251#ifdef SET_INEXACT
3252 bc.inexact = 0;
3253#endif
3254 break;
3255 }
3256 delta = lshift(delta,Log2P);
3257 if (cmp(delta, bs) > 0)
3258 goto drop_down;
3259 break;
3260 }
3261 if (i == 0) {
3262 /* exactly half-way between */
3263 if (bc.dsign) {
3264 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3265 && word1(&rv) == (
3266#ifdef Avoid_Underflow
3267 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3268 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3269#endif
3270 0xffffffff)) {
3271 /*boundary case -- increment exponent*/
3272 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3273 goto ovfl;
3274 word0(&rv) = (word0(&rv) & Exp_mask)
3275 + Exp_msk1
3276#ifdef IBM
3277 | Exp_msk1 >> 4
3278#endif
3279 ;
3280 word1(&rv) = 0;
3281#ifdef Avoid_Underflow
3282 bc.dsign = 0;
3283#endif
3284 break;
3285 }
3286 }
3287 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3288 drop_down:
3289 /* boundary case -- decrement exponent */
3290#ifdef Sudden_Underflow /*{{*/
3291 L = word0(&rv) & Exp_mask;
3292#ifdef IBM
3293 if (L < Exp_msk1)
3294#else
3295#ifdef Avoid_Underflow
3296 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3297#else
3298 if (L <= Exp_msk1)
3299#endif /*Avoid_Underflow*/
3300#endif /*IBM*/
3301 {
3302 if (bc.nd >nd) {
3303 bc.uflchk = 1;
3304 break;
3305 }
3306 goto undfl;
3307 }
3308 L -= Exp_msk1;
3309#else /*Sudden_Underflow}{*/
3310#ifdef Avoid_Underflow
3311 if (bc.scale) {
3312 L = word0(&rv) & Exp_mask;
3313 if (L <= (2*P+1)*Exp_msk1) {
3314 if (L > (P+2)*Exp_msk1)
3315 /* round even ==> */
3316 /* accept rv */
3317 break;
3318 /* rv = smallest denormal */
3319 if (bc.nd >nd) {
3320 bc.uflchk = 1;
3321 break;
3322 }
3323 goto undfl;
3324 }
3325 }
3326#endif /*Avoid_Underflow*/
3327 L = (word0(&rv) & Exp_mask) - Exp_msk1;
3328#endif /*Sudden_Underflow}}*/
3329 word0(&rv) = L | Bndry_mask1;
3330 word1(&rv) = 0xffffffff;
3331#ifdef IBM
3332 goto cont;
3333#else
3334#ifndef NO_STRTOD_BIGCOMP
3335 if (bc.nd > nd)
3336 goto cont;
3337#endif
3338 break;
3339#endif
3340 }
3341#ifndef ROUND_BIASED
3342#ifdef Avoid_Underflow
3343 if (Lsb1) {
3344 if (!(word0(&rv) & Lsb1))
3345 break;
3346 }
3347 else if (!(word1(&rv) & Lsb))
3348 break;
3349#else
3350 if (!(word1(&rv) & LSB))
3351 break;
3352#endif
3353#endif
3354 if (bc.dsign)
3355#ifdef Avoid_Underflow
3356 dval(&rv) += sulp(&rv, &bc);
3357#else
3358 dval(&rv) += ulp(&rv);
3359#endif
3360#ifndef ROUND_BIASED
3361 else {
3362#ifdef Avoid_Underflow
3363 dval(&rv) -= sulp(&rv, &bc);
3364#else
3365 dval(&rv) -= ulp(&rv);
3366#endif
3367#ifndef Sudden_Underflow
3368 if (!dval(&rv)) {
3369 if (bc.nd >nd) {
3370 bc.uflchk = 1;
3371 break;
3372 }
3373 goto undfl;
3374 }
3375#endif
3376 }
3377#ifdef Avoid_Underflow
3378 bc.dsign = 1 - bc.dsign;
3379#endif
3380#endif
3381 break;
3382 }
3383 if ((aadj = ratio(delta, bs)) <= 2.) {
3384 if (bc.dsign)
3385 aadj = aadj1 = 1.;
3386 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3387#ifndef Sudden_Underflow
3388 if (word1(&rv) == Tiny1 && !word0(&rv)) {
3389 if (bc.nd >nd) {
3390 bc.uflchk = 1;
3391 break;
3392 }
3393 goto undfl;
3394 }
3395#endif
3396 aadj = 1.;
3397 aadj1 = -1.;
3398 }
3399 else {
3400 /* special case -- power of FLT_RADIX to be */
3401 /* rounded down... */
3402
3403 if (aadj < 2./FLT_RADIX)
3404 aadj = 1./FLT_RADIX;
3405 else
3406 aadj *= 0.5;
3407 aadj1 = -aadj;
3408 }
3409 }
3410 else {
3411 aadj *= 0.5;
3412 aadj1 = bc.dsign ? aadj : -aadj;
3413#ifdef Check_FLT_ROUNDS
3414 switch(bc.rounding) {
3415 case 2: /* towards +infinity */
3416 aadj1 -= 0.5;
3417 break;
3418 case 0: /* towards 0 */
3419 case 3: /* towards -infinity */
3420 aadj1 += 0.5;
3421 }
3422#else
3423 if (Flt_Rounds == 0)
3424 aadj1 += 0.5;
3425#endif /*Check_FLT_ROUNDS*/
3426 }
3427 y = word0(&rv) & Exp_mask;
3428
3429 /* Check for overflow */
3430
3431 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3432 dval(&rv0) = dval(&rv);
3433 word0(&rv) -= P*Exp_msk1;
3434 adj.d = aadj1 * ulp(&rv);
3435 dval(&rv) += adj.d;
3436 if ((word0(&rv) & Exp_mask) >=
3437 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3438 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3439 goto ovfl;
3440 word0(&rv) = Big0;
3441 word1(&rv) = Big1;
3442 goto cont;
3443 }
3444 else
3445 word0(&rv) += P*Exp_msk1;
3446 }
3447 else {
3448#ifdef Avoid_Underflow
3449 if (bc.scale && y <= 2*P*Exp_msk1) {
3450 if (aadj <= 0x7fffffff) {
3451 if ((z = aadj) <= 0)
3452 z = 1;
3453 aadj = z;
3454 aadj1 = bc.dsign ? aadj : -aadj;
3455 }
3456 dval(&aadj2) = aadj1;
3457 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3458 aadj1 = dval(&aadj2);
3459 adj.d = aadj1 * ulp(&rv);
3460 dval(&rv) += adj.d;
3461 if (rv.d == 0.)
3462#ifdef NO_STRTOD_BIGCOMP
3463 goto undfl;
3464#else
3465 {
3466 if (bc.nd > nd)
3467 bc.dsign = 1;
3468 break;
3469 }
3470#endif
3471 }
3472 else {
3473 adj.d = aadj1 * ulp(&rv);
3474 dval(&rv) += adj.d;
3475 }
3476#else
3477#ifdef Sudden_Underflow
3478 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3479 dval(&rv0) = dval(&rv);
3480 word0(&rv) += P*Exp_msk1;
3481 adj.d = aadj1 * ulp(&rv);
3482 dval(&rv) += adj.d;
3483#ifdef IBM
3484 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
3485#else
3486 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3487#endif
3488 {
3489 if (word0(&rv0) == Tiny0
3490 && word1(&rv0) == Tiny1) {
3491 if (bc.nd >nd) {
3492 bc.uflchk = 1;
3493 break;
3494 }
3495 goto undfl;
3496 }
3497 word0(&rv) = Tiny0;
3498 word1(&rv) = Tiny1;
3499 goto cont;
3500 }
3501 else
3502 word0(&rv) -= P*Exp_msk1;
3503 }
3504 else {
3505 adj.d = aadj1 * ulp(&rv);
3506 dval(&rv) += adj.d;
3507 }
3508#else /*Sudden_Underflow*/
3509 /* Compute adj so that the IEEE rounding rules will
3510 * correctly round rv + adj in some half-way cases.
3511 * If rv * ulp(rv) is denormalized (i.e.,
3512 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3513 * trouble from bits lost to denormalization;
3514 * example: 1.2e-307 .
3515 */
3516 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3517 aadj1 = (double)(int)(aadj + 0.5);
3518 if (!bc.dsign)
3519 aadj1 = -aadj1;
3520 }
3521 adj.d = aadj1 * ulp(&rv);
3522 dval(&rv) += adj.d;
3523#endif /*Sudden_Underflow*/
3524#endif /*Avoid_Underflow*/
3525 }
3526 z = word0(&rv) & Exp_mask;
3527#ifndef SET_INEXACT
3528 if (bc.nd == nd) {
3529#ifdef Avoid_Underflow
3530 if (!bc.scale)
3531#endif
3532 if (y == z) {
3533 /* Can we stop now? */
3534 L = (Long)aadj;
3535 aadj -= L;
3536 /* The tolerances below are conservative. */
3537 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3538 if (aadj < .4999999 || aadj > .5000001)
3539 break;
3540 }
3541 else if (aadj < .4999999/FLT_RADIX)
3542 break;
3543 }
3544 }
3545#endif
3546 cont:
3547 Bfree(bb);
3548 Bfree(bd);
3549 Bfree(bs);
3550 Bfree(delta);
3551 }
3552 Bfree(bb);
3553 Bfree(bd);
3554 Bfree(bs);
3555 Bfree(bd0);
3556 Bfree(delta);
3557#ifndef NO_STRTOD_BIGCOMP
3558 if (req_bigcomp) {
3559 bd0 = 0;
3560 bc.e0 += nz1;
3561 bigcomp(&rv, s0, &bc);
3562 y = word0(&rv) & Exp_mask;
3563 if (y == Exp_mask)
3564 goto ovfl;
3565 if (y == 0 && rv.d == 0.)
3566 goto undfl;
3567 }
3568#endif
3569#ifdef SET_INEXACT
3570 if (bc.inexact) {
3571 if (!oldinexact) {
3572 word0(&rv0) = Exp_1 + (70 << Exp_shift);
3573 word1(&rv0) = 0;
3574 dval(&rv0) += 1.;
3575 }
3576 }
3577 else if (!oldinexact)
3578 clear_inexact();
3579#endif
3580#ifdef Avoid_Underflow
3581 if (bc.scale) {
3582 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3583 word1(&rv0) = 0;
3584 dval(&rv) *= dval(&rv0);
3585#ifndef NO_ERRNO
3586 /* try to avoid the bug of testing an 8087 register value */
3587#ifdef IEEE_Arith
3588 if (!(word0(&rv) & Exp_mask))
3589#else
3590 if (word0(&rv) == 0 && word1(&rv) == 0)
3591#endif
3592 errno = ERANGE;
3593#endif
3594 }
3595#endif /* Avoid_Underflow */
3596#ifdef SET_INEXACT
3597 if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3598 /* set underflow bit */
3599 dval(&rv0) = 1e-300;
3600 dval(&rv0) *= dval(&rv0);
3601 }
3602#endif
3603 ret:
3604 if (se)
3605 *se = (char *)s;
3606 return sign ? -dval(&rv) : dval(&rv);
3607 }
3608
3609#ifndef MULTIPLE_THREADS
3610 static char *dtoa_result;
3611#endif
3612
3613 static char *
3614#ifdef KR_headers
3615rv_alloc(i) int i;
3616#else
3618#endif
3619{
3620 int j, k, *r;
3621
3622 j = sizeof(ULong);
3623 for(k = 0;
3624 (int)(sizeof(Bigint) - sizeof(ULong) - sizeof(int)) + j <= i;
3625 j <<= 1)
3626 k++;
3627 r = (int*)Balloc(k);
3628 *r = k;
3629 return
3630#ifndef MULTIPLE_THREADS
3631 dtoa_result =
3632#endif
3633 (char *)(r+1);
3634 }
3635
3636 static char *
3637#ifdef KR_headers
3638nrv_alloc(s, rve, n) char *s, **rve; int n;
3639#else
3640nrv_alloc(const char *s, char **rve, int n)
3641#endif
3642{
3643 char *rv, *t;
3644
3645 t = rv = rv_alloc(n);
3646 while((*t = *s++)) t++;
3647 if (rve)
3648 *rve = t;
3649 return rv;
3650 }
3651
3652/* freedtoa(s) must be used to free values s returned by dtoa
3653 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3654 * but for consistency with earlier versions of dtoa, it is optional
3655 * when MULTIPLE_THREADS is not defined.
3656 */
3657
3658 void
3659#ifdef KR_headers
3660os_freedtoa(s) char *s;
3661#else
3663#endif
3664{
3665 Bigint *b = (Bigint *)((int *)s - 1);
3666 b->maxwds = 1 << (b->k = *(int*)b);
3667 Bfree(b);
3668#ifndef MULTIPLE_THREADS
3669 if (s == dtoa_result)
3670 dtoa_result = 0;
3671#endif
3672 }
3673
3674/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3675 *
3676 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3677 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3678 *
3679 * Modifications:
3680 * 1. Rather than iterating, we use a simple numeric overestimate
3681 * to determine k = floor(log10(d)). We scale relevant
3682 * quantities using O(log2(k)) rather than O(k) multiplications.
3683 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3684 * try to generate digits strictly left to right. Instead, we
3685 * compute with fewer bits and propagate the carry if necessary
3686 * when rounding the final digit up. This is often faster.
3687 * 3. Under the assumption that input will be rounded nearest,
3688 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3689 * That is, we allow equality in stopping tests when the
3690 * round-nearest rule will give the same floating-point value
3691 * as would satisfaction of the stopping test with strict
3692 * inequality.
3693 * 4. We remove common factors of powers of 2 from relevant
3694 * quantities.
3695 * 5. When converting floating-point integers less than 1e16,
3696 * we use floating-point arithmetic rather than resorting
3697 * to multiple-precision integers.
3698 * 6. When asked to produce fewer than 15 digits, we first try
3699 * to get by with floating-point arithmetic; we resort to
3700 * multiple-precision integer arithmetic only if we cannot
3701 * guarantee that the floating-point calculation has given
3702 * the correctly rounded result. For k requested digits and
3703 * "uniformly" distributed input, the probability is
3704 * something like 10^(k-15) that we must resort to the Long
3705 * calculation.
3706 */
3707
3708 char *
3710#ifdef KR_headers
3711 (dd, mode, ndigits, decpt, sign, rve)
3712 double dd; int mode, ndigits, *decpt, *sign; char **rve;
3713#else
3714 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
3715#endif
3716{
3717 /* Arguments ndigits, decpt, sign are similar to those
3718 of ecvt and fcvt; trailing zeros are suppressed from
3719 the returned string. If not null, *rve is set to point
3720 to the end of the return value. If d is +-Infinity or NaN,
3721 then *decpt is set to 9999.
3722
3723 mode:
3724 0 ==> shortest string that yields d when read in
3725 and rounded to nearest.
3726 1 ==> like 0, but with Steele & White stopping rule;
3727 e.g. with IEEE P754 arithmetic , mode 0 gives
3728 1e23 whereas mode 1 gives 9.999999999999999e22.
3729 2 ==> max(1,ndigits) significant digits. This gives a
3730 return value similar to that of ecvt, except
3731 that trailing zeros are suppressed.
3732 3 ==> through ndigits past the decimal point. This
3733 gives a return value similar to that from fcvt,
3734 except that trailing zeros are suppressed, and
3735 ndigits can be negative.
3736 4,5 ==> similar to 2 and 3, respectively, but (in
3737 round-nearest mode) with the tests of mode 0 to
3738 possibly return a shorter string that rounds to d.
3739 With IEEE arithmetic and compilation with
3740 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3741 as modes 2 and 3 when FLT_ROUNDS != 1.
3742 6-9 ==> Debugging modes similar to mode - 4: don't try
3743 fast floating-point estimate (if applicable).
3744
3745 Values of mode other than 0-9 are treated as mode 0.
3746
3747 Sufficient space is allocated to the return value
3748 to hold the suppressed trailing zeros.
3749 */
3750
3751 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3752 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3753 spec_case, try_quick;
3754 Long L;
3755#ifndef Sudden_Underflow
3756 int denorm;
3757 ULong x;
3758#endif
3759 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3760 U d2, eps, u;
3761 double ds;
3762 char *s, *s0;
3763#ifndef No_leftright
3764#ifdef IEEE_Arith
3765 U eps1;
3766#endif
3767#endif
3768#ifdef SET_INEXACT
3769 int inexact, oldinexact;
3770#endif
3771#ifdef Honor_FLT_ROUNDS /*{*/
3772 int Rounding;
3773#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3775#else /*}{*/
3776 Rounding = 1;
3777 switch(fegetround()) {
3778 case FE_TOWARDZERO: Rounding = 0; break;
3779 case FE_UPWARD: Rounding = 2; break;
3780 case FE_DOWNWARD: Rounding = 3;
3781 }
3782#endif /*}}*/
3783#endif /*}*/
3784
3785#ifndef MULTIPLE_THREADS
3786 if (dtoa_result) {
3788 dtoa_result = 0;
3789 }
3790#endif
3791
3792 u.d = dd;
3793 if (word0(&u) & Sign_bit) {
3794 /* set sign for everything, including 0's and NaNs */
3795 *sign = 1;
3796 word0(&u) &= ~Sign_bit; /* clear sign bit */
3797 }
3798 else
3799 *sign = 0;
3800
3801#if defined(IEEE_Arith) + defined(VAX)
3802#ifdef IEEE_Arith
3803 if ((word0(&u) & Exp_mask) == Exp_mask)
3804#else
3805 if (word0(&u) == 0x8000)
3806#endif
3807 {
3808 /* Infinity or NaN */
3809 *decpt = 9999;
3810#ifdef IEEE_Arith
3811 if (!word1(&u) && !(word0(&u) & 0xfffff))
3812 return nrv_alloc("Infinity", rve, 8);
3813#endif
3814 return nrv_alloc("NaN", rve, 3);
3815 }
3816#endif
3817#ifdef IBM
3818 dval(&u) += 0; /* normalize */
3819#endif
3820 if (!dval(&u)) {
3821 *decpt = 1;
3822 return nrv_alloc("0", rve, 1);
3823 }
3824
3825#ifdef SET_INEXACT
3826 try_quick = oldinexact = get_inexact();
3827 inexact = 1;
3828#endif
3829#ifdef Honor_FLT_ROUNDS
3830 if (Rounding >= 2) {
3831 if (*sign)
3832 Rounding = Rounding == 2 ? 0 : 2;
3833 else
3834 if (Rounding != 2)
3835 Rounding = 0;
3836 }
3837#endif
3838
3839 b = d2b(&u, &be, &bbits);
3840#ifdef Sudden_Underflow
3841 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3842#else
3843 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3844#endif
3845 dval(&d2) = dval(&u);
3846 word0(&d2) &= Frac_mask1;
3847 word0(&d2) |= Exp_11;
3848#ifdef IBM
3849 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3850 dval(&d2) /= 1 << j;
3851#endif
3852
3853 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3854 * log10(x) = log(x) / log(10)
3855 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3856 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3857 *
3858 * This suggests computing an approximation k to log10(d) by
3859 *
3860 * k = (i - Bias)*0.301029995663981
3861 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3862 *
3863 * We want k to be too large rather than too small.
3864 * The error in the first-order Taylor series approximation
3865 * is in our favor, so we just round up the constant enough
3866 * to compensate for any error in the multiplication of
3867 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3868 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3869 * adding 1e-13 to the constant term more than suffices.
3870 * Hence we adjust the constant term to 0.1760912590558.
3871 * (We could get a more accurate k by invoking log10,
3872 * but this is probably not worthwhile.)
3873 */
3874
3875 i -= Bias;
3876#ifdef IBM
3877 i <<= 2;
3878 i += j;
3879#endif
3880#ifndef Sudden_Underflow
3881 denorm = 0;
3882 }
3883 else {
3884 /* d is denormalized */
3885
3886 i = bbits + be + (Bias + (P-1) - 1);
3887 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3888 : word1(&u) << (32 - i);
3889 dval(&d2) = x;
3890 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3891 i -= (Bias + (P-1) - 1) + 1;
3892 denorm = 1;
3893 }
3894#endif
3895 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3896 k = (int)ds;
3897 if (ds < 0. && ds != k)
3898 k--; /* want k = floor(ds) */
3899 k_check = 1;
3900 if (k >= 0 && k <= Ten_pmax) {
3901 if (dval(&u) < tens[k])
3902 k--;
3903 k_check = 0;
3904 }
3905 j = bbits - i - 1;
3906 if (j >= 0) {
3907 b2 = 0;
3908 s2 = j;
3909 }
3910 else {
3911 b2 = -j;
3912 s2 = 0;
3913 }
3914 if (k >= 0) {
3915 b5 = 0;
3916 s5 = k;
3917 s2 += k;
3918 }
3919 else {
3920 b2 -= k;
3921 b5 = -k;
3922 s5 = 0;
3923 }
3924 if (mode < 0 || mode > 9)
3925 mode = 0;
3926
3927#ifndef SET_INEXACT
3928#ifdef Check_FLT_ROUNDS
3929 try_quick = Rounding == 1;
3930#else
3931 try_quick = 1;
3932#endif
3933#endif /*SET_INEXACT*/
3934
3935 if (mode > 5) {
3936 mode -= 4;
3937 try_quick = 0;
3938 }
3939 leftright = 1;
3940 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
3941 /* silence erroneous "gcc -Wall" warning. */
3942 switch(mode) {
3943 case 0:
3944 case 1:
3945 i = 18;
3946 ndigits = 0;
3947 break;
3948 case 2:
3949 leftright = 0;
3950 /* no break */
3951 case 4:
3952 if (ndigits <= 0)
3953 ndigits = 1;
3954 ilim = ilim1 = i = ndigits;
3955 break;
3956 case 3:
3957 leftright = 0;
3958 /* no break */
3959 case 5:
3960 i = ndigits + k + 1;
3961 ilim = i;
3962 ilim1 = i - 1;
3963 if (i <= 0)
3964 i = 1;
3965 }
3966 s = s0 = rv_alloc(i);
3967
3968#ifdef Honor_FLT_ROUNDS
3969 if (mode > 1 && Rounding != 1)
3970 leftright = 0;
3971#endif
3972
3973 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3974
3975 /* Try to get by with floating-point arithmetic. */
3976
3977 i = 0;
3978 dval(&d2) = dval(&u);
3979 k0 = k;
3980 ilim0 = ilim;
3981 ieps = 2; /* conservative */
3982 j1 = 0;
3983 if (k > 0) {
3984 ds = tens[k&0xf];
3985 j = k >> 4;
3986 if (j & Bletch) {
3987 /* prevent overflows */
3988 j &= Bletch - 1;
3989 dval(&u) /= bigtens[n_bigtens-1];
3990 ieps++;
3991 }
3992 for(; j; j >>= 1, i++)
3993 if (j & 1) {
3994 ieps++;
3995 ds *= bigtens[i];
3996 }
3997 dval(&u) /= ds;
3998 }
3999 else if ((j1 = -k)) {
4000 dval(&u) *= tens[j1 & 0xf];
4001 for(j = j1 >> 4; j; j >>= 1, i++)
4002 if (j & 1) {
4003 ieps++;
4004 dval(&u) *= bigtens[i];
4005 }
4006 }
4007 if (k_check && dval(&u) < 1. && ilim > 0) {
4008 if (ilim1 <= 0)
4009 goto fast_failed;
4010 ilim = ilim1;
4011 k--;
4012 dval(&u) *= 10.;
4013 ieps++;
4014 }
4015 dval(&eps) = ieps*dval(&u) + 7.;
4016 word0(&eps) -= (P-1)*Exp_msk1;
4017 if (ilim == 0) {
4018 S = mhi = 0;
4019 dval(&u) -= 5.;
4020 if (dval(&u) > dval(&eps))
4021 goto one_digit;
4022 if (dval(&u) < -dval(&eps))
4023 goto no_digits;
4024 goto fast_failed;
4025 }
4026#ifndef No_leftright
4027 if (leftright) {
4028 /* Use Steele & White method of only
4029 * generating digits needed.
4030 */
4031 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
4032#ifdef IEEE_Arith
4033 if (k0 < 0 && j1 >= 307) {
4034 eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
4035 word0(&eps1) -= Exp_msk1 * (Bias+P-1);
4036 dval(&eps1) *= tens[j1 & 0xf];
4037 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
4038 if (j & 1)
4039 dval(&eps1) *= bigtens[i];
4040 if (eps.d < eps1.d)
4041 eps.d = eps1.d;
4042 }
4043#endif
4044 for(i = 0;;) {
4045 L = dval(&u);
4046 dval(&u) -= L;
4047 *s++ = '0' + (int)L;
4048 if (1. - dval(&u) < dval(&eps))
4049 goto bump_up;
4050 if (dval(&u) < dval(&eps))
4051 goto ret1;
4052 if (++i >= ilim)
4053 break;
4054 dval(&eps) *= 10.;
4055 dval(&u) *= 10.;
4056 }
4057 }
4058 else {
4059#endif
4060 /* Generate ilim digits, then fix them up. */
4061 dval(&eps) *= tens[ilim-1];
4062 for(i = 1;; i++, dval(&u) *= 10.) {
4063 L = (Long)(dval(&u));
4064 if (!(dval(&u) -= L))
4065 ilim = i;
4066 *s++ = '0' + (int)L;
4067 if (i == ilim) {
4068 if (dval(&u) > 0.5 + dval(&eps))
4069 goto bump_up;
4070 else if (dval(&u) < 0.5 - dval(&eps)) {
4071 while(*--s == '0');
4072 s++;
4073 goto ret1;
4074 }
4075 break;
4076 }
4077 }
4078#ifndef No_leftright
4079 }
4080#endif
4081 fast_failed:
4082 s = s0;
4083 dval(&u) = dval(&d2);
4084 k = k0;
4085 ilim = ilim0;
4086 }
4087
4088 /* Do we have a "small" integer? */
4089
4090 if (be >= 0 && k <= Int_max) {
4091 /* Yes. */
4092 ds = tens[k];
4093 if (ndigits < 0 && ilim <= 0) {
4094 S = mhi = 0;
4095 if (ilim < 0 || dval(&u) <= 5*ds)
4096 goto no_digits;
4097 goto one_digit;
4098 }
4099 for(i = 1;; i++, dval(&u) *= 10.) {
4100 L = (Long)(dval(&u) / ds);
4101 dval(&u) -= L*ds;
4102#ifdef Check_FLT_ROUNDS
4103 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
4104 if (dval(&u) < 0) {
4105 L--;
4106 dval(&u) += ds;
4107 }
4108#endif
4109 *s++ = '0' + (int)L;
4110 if (!dval(&u)) {
4111#ifdef SET_INEXACT
4112 inexact = 0;
4113#endif
4114 break;
4115 }
4116 if (i == ilim) {
4117#ifdef Honor_FLT_ROUNDS
4118 if (mode > 1)
4119 switch(Rounding) {
4120 case 0: goto ret1;
4121 case 2: goto bump_up;
4122 }
4123#endif
4124 dval(&u) += dval(&u);
4125#ifdef ROUND_BIASED
4126 if (dval(&u) >= ds)
4127#else
4128 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
4129#endif
4130 {
4131 bump_up:
4132 while(*--s == '9')
4133 if (s == s0) {
4134 k++;
4135 *s = '0';
4136 break;
4137 }
4138 ++*s++;
4139 }
4140 break;
4141 }
4142 }
4143 goto ret1;
4144 }
4145
4146 m2 = b2;
4147 m5 = b5;
4148 mhi = mlo = 0;
4149 if (leftright) {
4150 i =
4151#ifndef Sudden_Underflow
4152 denorm ? be + (Bias + (P-1) - 1 + 1) :
4153#endif
4154#ifdef IBM
4155 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
4156#else
4157 1 + P - bbits;
4158#endif
4159 b2 += i;
4160 s2 += i;
4161 mhi = i2b(1);
4162 }
4163 if (m2 > 0 && s2 > 0) {
4164 i = m2 < s2 ? m2 : s2;
4165 b2 -= i;
4166 m2 -= i;
4167 s2 -= i;
4168 }
4169 if (b5 > 0) {
4170 if (leftright) {
4171 if (m5 > 0) {
4172 mhi = pow5mult(mhi, m5);
4173 b1 = mult(mhi, b);
4174 Bfree(b);
4175 b = b1;
4176 }
4177 if ((j = b5 - m5))
4178 b = pow5mult(b, j);
4179 }
4180 else
4181 b = pow5mult(b, b5);
4182 }
4183 S = i2b(1);
4184 if (s5 > 0)
4185 S = pow5mult(S, s5);
4186
4187 /* Check for special case that d is a normalized power of 2. */
4188
4189 spec_case = 0;
4190 if ((mode < 2 || leftright)
4191#ifdef Honor_FLT_ROUNDS
4192 && Rounding == 1
4193#endif
4194 ) {
4195 if (!word1(&u) && !(word0(&u) & Bndry_mask)
4196#ifndef Sudden_Underflow
4197 && word0(&u) & (Exp_mask & ~Exp_msk1)
4198#endif
4199 ) {
4200 /* The special case */
4201 b2 += Log2P;
4202 s2 += Log2P;
4203 spec_case = 1;
4204 }
4205 }
4206
4207 /* Arrange for convenient computation of quotients:
4208 * shift left if necessary so divisor has 4 leading 0 bits.
4209 *
4210 * Perhaps we should just compute leading 28 bits of S once
4211 * and for all and pass them and a shift to quorem, so it
4212 * can do shifts and ors to compute the numerator for q.
4213 */
4214 i = dshift(S, s2);
4215 b2 += i;
4216 m2 += i;
4217 s2 += i;
4218 if (b2 > 0)
4219 b = lshift(b, b2);
4220 if (s2 > 0)
4221 S = lshift(S, s2);
4222 if (k_check) {
4223 if (cmp(b,S) < 0) {
4224 k--;
4225 b = multadd(b, 10, 0); /* we botched the k estimate */
4226 if (leftright)
4227 mhi = multadd(mhi, 10, 0);
4228 ilim = ilim1;
4229 }
4230 }
4231 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4232 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4233 /* no digits, fcvt style */
4234 no_digits:
4235 k = -1 - ndigits;
4236 goto ret;
4237 }
4238 one_digit:
4239 *s++ = '1';
4240 k++;
4241 goto ret;
4242 }
4243 if (leftright) {
4244 if (m2 > 0)
4245 mhi = lshift(mhi, m2);
4246
4247 /* Compute mlo -- check for special case
4248 * that d is a normalized power of 2.
4249 */
4250
4251 mlo = mhi;
4252 if (spec_case) {
4253 mhi = Balloc(mhi->k);
4254 Bcopy(mhi, mlo);
4255 mhi = lshift(mhi, Log2P);
4256 }
4257
4258 for(i = 1;;i++) {
4259 dig = quorem(b,S) + '0';
4260 /* Do we yet have the shortest decimal string
4261 * that will round to d?
4262 */
4263 j = cmp(b, mlo);
4264 delta = diff(S, mhi);
4265 j1 = delta->sign ? 1 : cmp(b, delta);
4266 Bfree(delta);
4267#ifndef ROUND_BIASED
4268 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4269#ifdef Honor_FLT_ROUNDS
4270 && Rounding >= 1
4271#endif
4272 ) {
4273 if (dig == '9')
4274 goto round_9_up;
4275 if (j > 0)
4276 dig++;
4277#ifdef SET_INEXACT
4278 else if (!b->x[0] && b->wds <= 1)
4279 inexact = 0;
4280#endif
4281 *s++ = dig;
4282 goto ret;
4283 }
4284#endif
4285 if (j < 0 || (j == 0 && mode != 1
4286#ifndef ROUND_BIASED
4287 && !(word1(&u) & 1)
4288#endif
4289 )) {
4290 if (!b->x[0] && b->wds <= 1) {
4291#ifdef SET_INEXACT
4292 inexact = 0;
4293#endif
4294 goto accept_dig;
4295 }
4296#ifdef Honor_FLT_ROUNDS
4297 if (mode > 1)
4298 switch(Rounding) {
4299 case 0: goto accept_dig;
4300 case 2: goto keep_dig;
4301 }
4302#endif /*Honor_FLT_ROUNDS*/
4303 if (j1 > 0) {
4304 b = lshift(b, 1);
4305 j1 = cmp(b, S);
4306#ifdef ROUND_BIASED
4307 if (j1 >= 0 /*)*/
4308#else
4309 if ((j1 > 0 || (j1 == 0 && dig & 1))
4310#endif
4311 && dig++ == '9')
4312 goto round_9_up;
4313 }
4314 accept_dig:
4315 *s++ = dig;
4316 goto ret;
4317 }
4318 if (j1 > 0) {
4319#ifdef Honor_FLT_ROUNDS
4320 if (!Rounding)
4321 goto accept_dig;
4322#endif
4323 if (dig == '9') { /* possible if i == 1 */
4324 round_9_up:
4325 *s++ = '9';
4326 goto roundoff;
4327 }
4328 *s++ = dig + 1;
4329 goto ret;
4330 }
4331#ifdef Honor_FLT_ROUNDS
4332 keep_dig:
4333#endif
4334 *s++ = dig;
4335 if (i == ilim)
4336 break;
4337 b = multadd(b, 10, 0);
4338 if (mlo == mhi)
4339 mlo = mhi = multadd(mhi, 10, 0);
4340 else {
4341 mlo = multadd(mlo, 10, 0);
4342 mhi = multadd(mhi, 10, 0);
4343 }
4344 }
4345 }
4346 else
4347 for(i = 1;; i++) {
4348 *s++ = dig = quorem(b,S) + '0';
4349 if (!b->x[0] && b->wds <= 1) {
4350#ifdef SET_INEXACT
4351 inexact = 0;
4352#endif
4353 goto ret;
4354 }
4355 if (i >= ilim)
4356 break;
4357 b = multadd(b, 10, 0);
4358 }
4359
4360 /* Round off last digit */
4361
4362#ifdef Honor_FLT_ROUNDS
4363 switch(Rounding) {
4364 case 0: goto trimzeros;
4365 case 2: goto roundoff;
4366 }
4367#endif
4368 b = lshift(b, 1);
4369 j = cmp(b, S);
4370#ifdef ROUND_BIASED
4371 if (j >= 0)
4372#else
4373 if (j > 0 || (j == 0 && dig & 1))
4374#endif
4375 {
4376 roundoff:
4377 while(*--s == '9')
4378 if (s == s0) {
4379 k++;
4380 *s++ = '1';
4381 goto ret;
4382 }
4383 ++*s++;
4384 }
4385 else {
4386#ifdef Honor_FLT_ROUNDS
4387 trimzeros:
4388#endif
4389 while(*--s == '0');
4390 s++;
4391 }
4392 ret:
4393 Bfree(S);
4394 if (mhi) {
4395 if (mlo && mlo != mhi)
4396 Bfree(mlo);
4397 Bfree(mhi);
4398 }
4399 ret1:
4400#ifdef SET_INEXACT
4401 if (inexact) {
4402 if (!oldinexact) {
4403 word0(&u) = Exp_1 + (70 << Exp_shift);
4404 word1(&u) = 0;
4405 dval(&u) += 1.;
4406 }
4407 }
4408 else if (!oldinexact)
4409 clear_inexact();
4410#endif
4411 Bfree(b);
4412 *s = 0;
4413 *decpt = k + 1;
4414 if (rve)
4415 *rve = s;
4416 return s0;
4417 }
4418#ifdef __cplusplus
4419}
4420#endif
void free(void *)
#define Emax
Definition OSdtoa.cpp:402
#define Exp_shift
Definition OSdtoa.cpp:394
#define rounded_product(a, b)
Definition OSdtoa.cpp:527
#define Exp_msk11
Definition OSdtoa.cpp:397
static void bigcomp(U *rv, const char *s0, BCinfo *bc)
Definition OSdtoa.cpp:2333
#define CONST
Definition OSdtoa.cpp:346
#define Tiny0
Definition OSdtoa.cpp:416
static double b2d(Bigint *a, int *e)
Definition OSdtoa.cpp:1256
#define Big1
Definition OSdtoa.cpp:532
#define Exp_msk1
Definition OSdtoa.cpp:396
#define FREE_DTOA_LOCK(n)
Definition OSdtoa.cpp:569
#define d1
static Bigint * pow5mult(Bigint *b, int k)
Definition OSdtoa.cpp:969
#define P
Definition OSdtoa.cpp:399
static int lo0bits(ULong *y)
Definition OSdtoa.cpp:793
#define NAN_WORD1
Definition OSdtoa.cpp:1596
#define IEEE_Arith
Definition OSdtoa.cpp:286
#define Nbits
Definition OSdtoa.cpp:400
#define Bndry_mask
Definition OSdtoa.cpp:411
#define Ebits
Definition OSdtoa.cpp:406
static void rshift(Bigint *b, int k)
Definition OSdtoa.cpp:1739
#define rounded_quotient(a, b)
Definition OSdtoa.cpp:528
static Bigint * p5s
Definition OSdtoa.cpp:962
static char * rv_alloc(int i)
Definition OSdtoa.cpp:3617
#define word1(x)
Definition OSdtoa.cpp:358
static Bigint * diff(Bigint *a, Bigint *b)
Definition OSdtoa.cpp:1120
#define MALLOC
Definition OSdtoa.cpp:269
#define Frac_mask
Definition OSdtoa.cpp:407
static Bigint * i2b(int i)
Definition OSdtoa.cpp:841
#define Tiny1
Definition OSdtoa.cpp:417
#define Frac_mask1
Definition OSdtoa.cpp:408
#define Long
Definition OSdtoa.cpp:37
double os_strtod(const char *s00, char **se)
Definition OSdtoa.cpp:2541
#define Quick_max
Definition OSdtoa.cpp:418
#define Bias
Definition OSdtoa.cpp:401
#define PRIVATE_mem
Definition OSdtoa.cpp:276
#define n_bigtens
Definition OSdtoa.cpp:1520
static int quorem(Bigint *b, Bigint *S)
Definition OSdtoa.cpp:2185
#define LSB
Definition OSdtoa.cpp:413
void os_freedtoa(char *s)
Definition OSdtoa.cpp:3662
#define dval(x)
Definition OSdtoa.cpp:363
#define Exp_mask
Definition OSdtoa.cpp:398
#define Kmax
Definition OSdtoa.cpp:572
#define Ten_pmax
Definition OSdtoa.cpp:409
static CONST double tinytens[]
Definition OSdtoa.cpp:1509
static int match(const char **sp, const char *t)
Definition OSdtoa.cpp:1604
#define Sign_bit
Definition OSdtoa.cpp:414
@ Round_zero
Definition OSdtoa.cpp:1799
@ Round_up
Definition OSdtoa.cpp:1801
@ Round_near
Definition OSdtoa.cpp:1800
@ Round_down
Definition OSdtoa.cpp:1802
#define kmask
Definition OSdtoa.cpp:1696
static Bigint * d2b(U *d, int *e, int *bits)
Definition OSdtoa.cpp:1326
#define Bletch
Definition OSdtoa.cpp:410
static Bigint * multadd(Bigint *b, int m, int a)
Definition OSdtoa.cpp:664
static void Bfree(Bigint *v)
Definition OSdtoa.cpp:637
static Bigint * mult(Bigint *a, Bigint *b)
Definition OSdtoa.cpp:857
void os_gethex(CONST char **sp, U *rvp, int rounding, int sign)
Definition OSdtoa.cpp:1810
static Bigint * increment(Bigint *b)
Definition OSdtoa.cpp:1704
static char * nrv_alloc(const char *s, char **rve, int n)
Definition OSdtoa.cpp:3640
#define strtod_diglim
Definition OSdtoa.cpp:372
#define Flt_Rounds
Definition OSdtoa.cpp:431
#define Storeinc(a, b, c)
Definition OSdtoa.cpp:380
#define word0(x)
Definition OSdtoa.cpp:357
#define ACQUIRE_DTOA_LOCK(n)
Definition OSdtoa.cpp:568
#define Exp_1
Definition OSdtoa.cpp:404
static int dshift(Bigint *b, int p2)
Definition OSdtoa.cpp:2171
#define NAN_WORD0
Definition OSdtoa.cpp:1592
static void hexnan(U *rvp, const char **sp)
Definition OSdtoa.cpp:1626
char * os_dtoa(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
Definition OSdtoa.cpp:3714
#define Exp_shift1
Definition OSdtoa.cpp:395
static Bigint * freelist[Kmax+1]
Definition OSdtoa.cpp:589
#define Rounding
Definition OSdtoa.cpp:439
#define kshift
Definition OSdtoa.cpp:1695
#define Emin
Definition OSdtoa.cpp:403
static char * dtoa_result
Definition OSdtoa.cpp:3610
#define Bcopy(x, y)
Definition OSdtoa.cpp:656
static int hi0bits(ULong x)
Definition OSdtoa.cpp:759
static Bigint * Balloc(int k)
Definition OSdtoa.cpp:596
static ULong any_on(Bigint *b, int k)
Definition OSdtoa.cpp:1772
#define Bndry_mask1
Definition OSdtoa.cpp:412
#define FFFFFFFF
Definition OSdtoa.cpp:545
static double sulp(U *x, BCinfo *bc)
Definition OSdtoa.cpp:2310
#define Scale_Bit
Definition OSdtoa.cpp:1519
#define d0
static Bigint * lshift(Bigint *b, int k)
Definition OSdtoa.cpp:1026
static double ulp(U *x)
Definition OSdtoa.cpp:1214
static int cmp(Bigint *a, Bigint *b)
Definition OSdtoa.cpp:1086
#define Log2P
Definition OSdtoa.cpp:415
static double ratio(Bigint *a, Bigint *b)
Definition OSdtoa.cpp:1460
static CONST double tens[]
Definition OSdtoa.cpp:1497
#define Big0
Definition OSdtoa.cpp:531
static double private_mem[PRIVATE_mem]
Definition OSdtoa.cpp:277
static Bigint * s2b(const char *s, int nd0, int nd, ULong y9, int dplen)
Definition OSdtoa.cpp:721
#define Avoid_Underflow
Definition OSdtoa.cpp:421
#define Exp_11
Definition OSdtoa.cpp:405
static CONST double bigtens[]
Definition OSdtoa.cpp:1508
#define Int_max
Definition OSdtoa.cpp:419
struct Bigint Bigint
Definition OSdtoa.cpp:587
static double * pmem_next
Definition OSdtoa.cpp:277
#define ULbits
Definition OSdtoa.cpp:1694
unsigned Long ULong
end of OS code, below is David Gay, except as noted above
Definition OSdtoa.cpp:239
static unsigned char hexdig[256]
Definition OSdtoa.cpp:1568
int e0
Definition OSdtoa.cpp:540
int dp0
Definition OSdtoa.cpp:540
int scale
Definition OSdtoa.cpp:540
int uflchk
Definition OSdtoa.cpp:540
int dplen
Definition OSdtoa.cpp:540
int nd0
Definition OSdtoa.cpp:540
int dp1
Definition OSdtoa.cpp:540
int nd
Definition OSdtoa.cpp:540
int inexact
Definition OSdtoa.cpp:540
int dsign
Definition OSdtoa.cpp:540
int rounding
Definition OSdtoa.cpp:540
int k
Definition OSdtoa.cpp:583
struct Bigint * next
Definition OSdtoa.cpp:582
int sign
Definition OSdtoa.cpp:583
int maxwds
Definition OSdtoa.cpp:583
int wds
Definition OSdtoa.cpp:583
ULong x[1]
Definition OSdtoa.cpp:584
Definition OSdtoa.cpp:354
double d
Definition OSdtoa.cpp:354