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