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