1 /*-
2 * Copyright (c) 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.redist.c%
6 */
7
8 #if defined(LIBC_SCCS) && !defined(lint)
9 static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 06/04/93";
10 #endif /* LIBC_SCCS and not lint */
11
12 /****************************************************************
13 *
14 * The author of this software is David M. Gay.
15 *
16 * Copyright (c) 1991 by AT&T.
17 *
18 * Permission to use, copy, modify, and distribute this software for any
19 * purpose without fee is hereby granted, provided that this entire notice
20 * is included in all copies of any software which is or includes a copy
21 * or modification of this software and in all copies of the supporting
22 * documentation for such software.
23 *
24 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
25 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
26 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
27 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
28 *
29 ***************************************************************/
30
31 /* Please send bug reports to
32 David M. Gay
33 AT&T Bell Laboratories, Room 2C-463
34 600 Mountain Avenue
35 Murray Hill, NJ 07974-2070
36 U.S.A.
37 dmg@research.att.com or research!dmg
38 */
39
40 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
41 *
42 * This strtod returns a nearest machine number to the input decimal
43 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
44 * broken by the IEEE round-even rule. Otherwise ties are broken by
45 * biased rounding (add half and chop).
46 *
47 * Inspired loosely by William D. Clinger's paper "How to Read Floating
48 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
49 *
50 * Modifications:
51 *
52 * 1. We only require IEEE, IBM, or VAX double-precision
53 * arithmetic (not IEEE double-extended).
54 * 2. We get by with floating-point arithmetic in a case that
55 * Clinger missed -- when we're computing d * 10^n
56 * for a small integer d and the integer n is not too
57 * much larger than 22 (the maximum integer k for which
58 * we can represent 10^k exactly), we may be able to
59 * compute (d*10^k) * 10^(e-k) with just one roundoff.
60 * 3. Rather than a bit-at-a-time adjustment of the binary
61 * result in the hard case, we use floating-point
62 * arithmetic to determine the adjustment to within
63 * one bit; only in really hard cases do we need to
64 * compute a second residual.
65 * 4. Because of 3., we don't need a large table of powers of 10
66 * for ten-to-e (just some small tables, e.g. of 10^k
67 * for 0 <= k <= 22).
68 */
69
70 /*
71 * #define IEEE_8087 for IEEE-arithmetic machines where the least
72 * significant byte has the lowest address.
73 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
74 * significant byte has the lowest address.
75 * #define Sudden_Underflow for IEEE-format machines without gradual
76 * underflow (i.e., that flush to zero on underflow).
77 * #define IBM for IBM mainframe-style floating-point arithmetic.
78 * #define VAX for VAX-style floating-point arithmetic.
79 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
80 * #define No_leftright to omit left-right logic in fast floating-point
81 * computation of dtoa.
82 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
83 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
84 * that use extended-precision instructions to compute rounded
85 * products and quotients) with IBM.
86 * #define ROUND_BIASED for IEEE-format with biased rounding.
87 * #define Inaccurate_Divide for IEEE-format with correctly rounded
88 * products but inaccurate quotients, e.g., for Intel i860.
89 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
90 * integer arithmetic. Whether this speeds things up or slows things
91 * down depends on the machine and the number being converted.
92 * #define KR_headers for old-style C function headers.
93 * #define Bad_float_h if your system lacks a float.h or if it does not
94 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
95 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
96 */
97
98 #if defined(i386) || defined(mips) && defined(MIPSEL)
99 #define IEEE_8087
100 #else
101 #define IEEE_MC68k
102 #endif
103
104 #ifdef DEBUG
105 #include "stdio.h"
106 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
107 #endif
108
109 #ifdef __cplusplus
110 #include "malloc.h"
111 #include "memory.h"
112 #else
113 #ifndef KR_headers
114 #include "stdlib.h"
115 #include "string.h"
116 #else
117 #include "malloc.h"
118 #include "memory.h"
119 #endif
120 #endif
121
122 #include "errno.h"
123 #ifdef Bad_float_h
124 #undef __STDC__
125 #ifdef IEEE_MC68k
126 #define IEEE_ARITHMETIC
127 #endif
128 #ifdef IEEE_8087
129 #define IEEE_ARITHMETIC
130 #endif
131 #ifdef IEEE_ARITHMETIC
132 #define DBL_DIG 15
133 #define DBL_MAX_10_EXP 308
134 #define DBL_MAX_EXP 1024
135 #define FLT_RADIX 2
136 #define FLT_ROUNDS 1
137 #define DBL_MAX 1.7976931348623157e+308
138 #endif
139
140 #ifdef IBM
141 #define DBL_DIG 16
142 #define DBL_MAX_10_EXP 75
143 #define DBL_MAX_EXP 63
144 #define FLT_RADIX 16
145 #define FLT_ROUNDS 0
146 #define DBL_MAX 7.2370055773322621e+75
147 #endif
148
149 #ifdef VAX
150 #define DBL_DIG 16
151 #define DBL_MAX_10_EXP 38
152 #define DBL_MAX_EXP 127
153 #define FLT_RADIX 2
154 #define FLT_ROUNDS 1
155 #define DBL_MAX 1.7014118346046923e+38
156 #endif
157
158 #ifndef LONG_MAX
159 #define LONG_MAX 2147483647
160 #endif
161 #else
162 #include "float.h"
163 #endif
164 #ifndef __MATH_H__
165 #include "math.h"
166 #endif
167
168 #ifdef __cplusplus
169 extern "C" {
170 #endif
171
172 #ifndef CONST
173 #ifdef KR_headers
174 #define CONST /* blank */
175 #else
176 #define CONST const
177 #endif
178 #endif
179
180 #ifdef Unsigned_Shifts
181 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
182 #else
183 #define Sign_Extend(a,b) /*no-op*/
184 #endif
185
186 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
187 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
188 #endif
189
190 #ifdef IEEE_8087
191 #define word0(x) ((unsigned long *)&x)[1]
192 #define word1(x) ((unsigned long *)&x)[0]
193 #else
194 #define word0(x) ((unsigned long *)&x)[0]
195 #define word1(x) ((unsigned long *)&x)[1]
196 #endif
197
198 /* The following definition of Storeinc is appropriate for MIPS processors.
199 * An alternative that might be better on some machines is
200 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
201 */
202 #if defined(IEEE_8087) + defined(VAX)
203 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
204 ((unsigned short *)a)[0] = (unsigned short)c, a++)
205 #else
206 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
207 ((unsigned short *)a)[1] = (unsigned short)c, a++)
208 #endif
209
210 /* #define P DBL_MANT_DIG */
211 /* Ten_pmax = floor(P*log(2)/log(5)) */
212 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
213 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
214 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
215
216 #if defined(IEEE_8087) + defined(IEEE_MC68k)
217 #define Exp_shift 20
218 #define Exp_shift1 20
219 #define Exp_msk1 0x100000
220 #define Exp_msk11 0x100000
221 #define Exp_mask 0x7ff00000
222 #define P 53
223 #define Bias 1023
224 #define IEEE_Arith
225 #define Emin (-1022)
226 #define Exp_1 0x3ff00000
227 #define Exp_11 0x3ff00000
228 #define Ebits 11
229 #define Frac_mask 0xfffff
230 #define Frac_mask1 0xfffff
231 #define Ten_pmax 22
232 #define Bletch 0x10
233 #define Bndry_mask 0xfffff
234 #define Bndry_mask1 0xfffff
235 #define LSB 1
236 #define Sign_bit 0x80000000
237 #define Log2P 1
238 #define Tiny0 0
239 #define Tiny1 1
240 #define Quick_max 14
241 #define Int_max 14
242 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
243 #else
244 #undef Sudden_Underflow
245 #define Sudden_Underflow
246 #ifdef IBM
247 #define Exp_shift 24
248 #define Exp_shift1 24
249 #define Exp_msk1 0x1000000
250 #define Exp_msk11 0x1000000
251 #define Exp_mask 0x7f000000
252 #define P 14
253 #define Bias 65
254 #define Exp_1 0x41000000
255 #define Exp_11 0x41000000
256 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
257 #define Frac_mask 0xffffff
258 #define Frac_mask1 0xffffff
259 #define Bletch 4
260 #define Ten_pmax 22
261 #define Bndry_mask 0xefffff
262 #define Bndry_mask1 0xffffff
263 #define LSB 1
264 #define Sign_bit 0x80000000
265 #define Log2P 4
266 #define Tiny0 0x100000
267 #define Tiny1 0
268 #define Quick_max 14
269 #define Int_max 15
270 #else /* VAX */
271 #define Exp_shift 23
272 #define Exp_shift1 7
273 #define Exp_msk1 0x80
274 #define Exp_msk11 0x800000
275 #define Exp_mask 0x7f80
276 #define P 56
277 #define Bias 129
278 #define Exp_1 0x40800000
279 #define Exp_11 0x4080
280 #define Ebits 8
281 #define Frac_mask 0x7fffff
282 #define Frac_mask1 0xffff007f
283 #define Ten_pmax 24
284 #define Bletch 2
285 #define Bndry_mask 0xffff007f
286 #define Bndry_mask1 0xffff007f
287 #define LSB 0x10000
288 #define Sign_bit 0x8000
289 #define Log2P 1
290 #define Tiny0 0x80
291 #define Tiny1 0
292 #define Quick_max 15
293 #define Int_max 15
294 #endif
295 #endif
296
297 #ifndef IEEE_Arith
298 #define ROUND_BIASED
299 #endif
300
301 #ifdef RND_PRODQUOT
302 #define rounded_product(a,b) a = rnd_prod(a, b)
303 #define rounded_quotient(a,b) a = rnd_quot(a, b)
304 #ifdef KR_headers
305 extern double rnd_prod(), rnd_quot();
306 #else
307 extern double rnd_prod(double, double), rnd_quot(double, double);
308 #endif
309 #else
310 #define rounded_product(a,b) a *= b
311 #define rounded_quotient(a,b) a /= b
312 #endif
313
314 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
315 #define Big1 0xffffffff
316
317 #ifndef Just_16
318 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
319 * This makes some inner loops simpler and sometimes saves work
320 * during multiplications, but it often seems to make things slightly
321 * slower. Hence the default is now to store 32 bits per long.
322 */
323 #ifndef Pack_32
324 #define Pack_32
325 #endif
326 #endif
327
328 #define Kmax 15
329
330 #ifdef __cplusplus
331 extern "C" double strtod(const char *s00, char **se);
332 extern "C" char *dtoa(double d, int mode, int ndigits,
333 int *decpt, int *sign, char **rve);
334 #endif
335
336 struct
337 Bigint {
338 struct Bigint *next;
339 int k, maxwds, sign, wds;
340 unsigned long x[1];
341 };
342
343 typedef struct Bigint Bigint;
344
345 static Bigint *freelist[Kmax+1];
346
347 static Bigint *
Balloc(k)348 Balloc
349 #ifdef KR_headers
350 (k) int k;
351 #else
352 (int k)
353 #endif
354 {
355 int x;
356 Bigint *rv;
357
358 if (rv = freelist[k]) {
359 freelist[k] = rv->next;
360 } else {
361 x = 1 << k;
362 rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
363 rv->k = k;
364 rv->maxwds = x;
365 }
366 rv->sign = rv->wds = 0;
367 return rv;
368 }
369
370 static void
Bfree(v)371 Bfree
372 #ifdef KR_headers
373 (v) Bigint *v;
374 #else
375 (Bigint *v)
376 #endif
377 {
378 if (v) {
379 v->next = freelist[v->k];
380 freelist[v->k] = v;
381 }
382 }
383
384 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
385 y->wds*sizeof(long) + 2*sizeof(int))
386
387 static Bigint *
multadd(b,m,a)388 multadd
389 #ifdef KR_headers
390 (b, m, a) Bigint *b; int m, a;
391 #else
392 (Bigint *b, int m, int a) /* multiply by m and add a */
393 #endif
394 {
395 int i, wds;
396 unsigned long *x, y;
397 #ifdef Pack_32
398 unsigned long xi, z;
399 #endif
400 Bigint *b1;
401
402 wds = b->wds;
403 x = b->x;
404 i = 0;
405 do {
406 #ifdef Pack_32
407 xi = *x;
408 y = (xi & 0xffff) * m + a;
409 z = (xi >> 16) * m + (y >> 16);
410 a = (int)(z >> 16);
411 *x++ = (z << 16) + (y & 0xffff);
412 #else
413 y = *x * m + a;
414 a = (int)(y >> 16);
415 *x++ = y & 0xffff;
416 #endif
417 } while (++i < wds);
418 if (a) {
419 if (wds >= b->maxwds) {
420 b1 = Balloc(b->k+1);
421 Bcopy(b1, b);
422 Bfree(b);
423 b = b1;
424 }
425 b->x[wds++] = a;
426 b->wds = wds;
427 }
428 return b;
429 }
430
431 static Bigint *
s2b(s,nd0,nd,y9)432 s2b
433 #ifdef KR_headers
434 (s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
435 #else
436 (CONST char *s, int nd0, int nd, unsigned long y9)
437 #endif
438 {
439 Bigint *b;
440 int i, k;
441 long x, y;
442
443 x = (nd + 8) / 9;
444 for (k = 0, y = 1; x > y; y <<= 1, k++) ;
445 #ifdef Pack_32
446 b = Balloc(k);
447 b->x[0] = y9;
448 b->wds = 1;
449 #else
450 b = Balloc(k+1);
451 b->x[0] = y9 & 0xffff;
452 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
453 #endif
454
455 i = 9;
456 if (9 < nd0) {
457 s += 9;
458 do
459 b = multadd(b, 10, *s++ - '0');
460 while (++i < nd0);
461 s++;
462 } else
463 s += 10;
464 for (; i < nd; i++)
465 b = multadd(b, 10, *s++ - '0');
466 return b;
467 }
468
469 static int
hi0bits(x)470 hi0bits
471 #ifdef KR_headers
472 (x) register unsigned long x;
473 #else
474 (register unsigned long x)
475 #endif
476 {
477 register int k = 0;
478
479 if (!(x & 0xffff0000)) {
480 k = 16;
481 x <<= 16;
482 }
483 if (!(x & 0xff000000)) {
484 k += 8;
485 x <<= 8;
486 }
487 if (!(x & 0xf0000000)) {
488 k += 4;
489 x <<= 4;
490 }
491 if (!(x & 0xc0000000)) {
492 k += 2;
493 x <<= 2;
494 }
495 if (!(x & 0x80000000)) {
496 k++;
497 if (!(x & 0x40000000))
498 return 32;
499 }
500 return k;
501 }
502
503 static int
lo0bits(y)504 lo0bits
505 #ifdef KR_headers
506 (y) unsigned long *y;
507 #else
508 (unsigned long *y)
509 #endif
510 {
511 register int k;
512 register unsigned long x = *y;
513
514 if (x & 7) {
515 if (x & 1)
516 return 0;
517 if (x & 2) {
518 *y = x >> 1;
519 return 1;
520 }
521 *y = x >> 2;
522 return 2;
523 }
524 k = 0;
525 if (!(x & 0xffff)) {
526 k = 16;
527 x >>= 16;
528 }
529 if (!(x & 0xff)) {
530 k += 8;
531 x >>= 8;
532 }
533 if (!(x & 0xf)) {
534 k += 4;
535 x >>= 4;
536 }
537 if (!(x & 0x3)) {
538 k += 2;
539 x >>= 2;
540 }
541 if (!(x & 1)) {
542 k++;
543 x >>= 1;
544 if (!x & 1)
545 return 32;
546 }
547 *y = x;
548 return k;
549 }
550
551 static Bigint *
i2b(i)552 i2b
553 #ifdef KR_headers
554 (i) int i;
555 #else
556 (int i)
557 #endif
558 {
559 Bigint *b;
560
561 b = Balloc(1);
562 b->x[0] = i;
563 b->wds = 1;
564 return b;
565 }
566
567 static Bigint *
mult(a,b)568 mult
569 #ifdef KR_headers
570 (a, b) Bigint *a, *b;
571 #else
572 (Bigint *a, Bigint *b)
573 #endif
574 {
575 Bigint *c;
576 int k, wa, wb, wc;
577 unsigned long carry, y, z;
578 unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
579 #ifdef Pack_32
580 unsigned long z2;
581 #endif
582
583 if (a->wds < b->wds) {
584 c = a;
585 a = b;
586 b = c;
587 }
588 k = a->k;
589 wa = a->wds;
590 wb = b->wds;
591 wc = wa + wb;
592 if (wc > a->maxwds)
593 k++;
594 c = Balloc(k);
595 for (x = c->x, xa = x + wc; x < xa; x++)
596 *x = 0;
597 xa = a->x;
598 xae = xa + wa;
599 xb = b->x;
600 xbe = xb + wb;
601 xc0 = c->x;
602 #ifdef Pack_32
603 for (; xb < xbe; xb++, xc0++) {
604 if (y = *xb & 0xffff) {
605 x = xa;
606 xc = xc0;
607 carry = 0;
608 do {
609 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
610 carry = z >> 16;
611 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
612 carry = z2 >> 16;
613 Storeinc(xc, z2, z);
614 } while (x < xae);
615 *xc = carry;
616 }
617 if (y = *xb >> 16) {
618 x = xa;
619 xc = xc0;
620 carry = 0;
621 z2 = *xc;
622 do {
623 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
624 carry = z >> 16;
625 Storeinc(xc, z, z2);
626 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
627 carry = z2 >> 16;
628 } while (x < xae);
629 *xc = z2;
630 }
631 }
632 #else
633 for (; xb < xbe; xc0++) {
634 if (y = *xb++) {
635 x = xa;
636 xc = xc0;
637 carry = 0;
638 do {
639 z = *x++ * y + *xc + carry;
640 carry = z >> 16;
641 *xc++ = z & 0xffff;
642 } while (x < xae);
643 *xc = carry;
644 }
645 }
646 #endif
647 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
648 c->wds = wc;
649 return c;
650 }
651
652 static Bigint *p5s;
653
654 static Bigint *
pow5mult(b,k)655 pow5mult
656 #ifdef KR_headers
657 (b, k) Bigint *b; int k;
658 #else
659 (Bigint *b, int k)
660 #endif
661 {
662 Bigint *b1, *p5, *p51;
663 int i;
664 static int p05[3] = { 5, 25, 125 };
665
666 if (i = k & 3)
667 b = multadd(b, p05[i-1], 0);
668
669 if (!(k >>= 2))
670 return b;
671 if (!(p5 = p5s)) {
672 /* first time */
673 p5 = p5s = i2b(625);
674 p5->next = 0;
675 }
676 for (;;) {
677 if (k & 1) {
678 b1 = mult(b, p5);
679 Bfree(b);
680 b = b1;
681 }
682 if (!(k >>= 1))
683 break;
684 if (!(p51 = p5->next)) {
685 p51 = p5->next = mult(p5,p5);
686 p51->next = 0;
687 }
688 p5 = p51;
689 }
690 return b;
691 }
692
693 static Bigint *
lshift(b,k)694 lshift
695 #ifdef KR_headers
696 (b, k) Bigint *b; int k;
697 #else
698 (Bigint *b, int k)
699 #endif
700 {
701 int i, k1, n, n1;
702 Bigint *b1;
703 unsigned long *x, *x1, *xe, z;
704
705 #ifdef Pack_32
706 n = k >> 5;
707 #else
708 n = k >> 4;
709 #endif
710 k1 = b->k;
711 n1 = n + b->wds + 1;
712 for (i = b->maxwds; n1 > i; i <<= 1)
713 k1++;
714 b1 = Balloc(k1);
715 x1 = b1->x;
716 for (i = 0; i < n; i++)
717 *x1++ = 0;
718 x = b->x;
719 xe = x + b->wds;
720 #ifdef Pack_32
721 if (k &= 0x1f) {
722 k1 = 32 - k;
723 z = 0;
724 do {
725 *x1++ = *x << k | z;
726 z = *x++ >> k1;
727 } while (x < xe);
728 if (*x1 = z)
729 ++n1;
730 }
731 #else
732 if (k &= 0xf) {
733 k1 = 16 - k;
734 z = 0;
735 do {
736 *x1++ = *x << k & 0xffff | z;
737 z = *x++ >> k1;
738 } while (x < xe);
739 if (*x1 = z)
740 ++n1;
741 }
742 #endif
743 else
744 do
745 *x1++ = *x++;
746 while (x < xe);
747 b1->wds = n1 - 1;
748 Bfree(b);
749 return b1;
750 }
751
752 static int
cmp(a,b)753 cmp
754 #ifdef KR_headers
755 (a, b) Bigint *a, *b;
756 #else
757 (Bigint *a, Bigint *b)
758 #endif
759 {
760 unsigned long *xa, *xa0, *xb, *xb0;
761 int i, j;
762
763 i = a->wds;
764 j = b->wds;
765 #ifdef DEBUG
766 if (i > 1 && !a->x[i-1])
767 Bug("cmp called with a->x[a->wds-1] == 0");
768 if (j > 1 && !b->x[j-1])
769 Bug("cmp called with b->x[b->wds-1] == 0");
770 #endif
771 if (i -= j)
772 return i;
773 xa0 = a->x;
774 xa = xa0 + j;
775 xb0 = b->x;
776 xb = xb0 + j;
777 for (;;) {
778 if (*--xa != *--xb)
779 return *xa < *xb ? -1 : 1;
780 if (xa <= xa0)
781 break;
782 }
783 return 0;
784 }
785
786 static Bigint *
diff(a,b)787 diff
788 #ifdef KR_headers
789 (a, b) Bigint *a, *b;
790 #else
791 (Bigint *a, Bigint *b)
792 #endif
793 {
794 Bigint *c;
795 int i, wa, wb;
796 long borrow, y; /* We need signed shifts here. */
797 unsigned long *xa, *xae, *xb, *xbe, *xc;
798 #ifdef Pack_32
799 long z;
800 #endif
801
802 i = cmp(a,b);
803 if (!i) {
804 c = Balloc(0);
805 c->wds = 1;
806 c->x[0] = 0;
807 return c;
808 }
809 if (i < 0) {
810 c = a;
811 a = b;
812 b = c;
813 i = 1;
814 } else
815 i = 0;
816 c = Balloc(a->k);
817 c->sign = i;
818 wa = a->wds;
819 xa = a->x;
820 xae = xa + wa;
821 wb = b->wds;
822 xb = b->x;
823 xbe = xb + wb;
824 xc = c->x;
825 borrow = 0;
826 #ifdef Pack_32
827 do {
828 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
829 borrow = y >> 16;
830 Sign_Extend(borrow, y);
831 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
832 borrow = z >> 16;
833 Sign_Extend(borrow, z);
834 Storeinc(xc, z, y);
835 } while (xb < xbe);
836 while (xa < xae) {
837 y = (*xa & 0xffff) + borrow;
838 borrow = y >> 16;
839 Sign_Extend(borrow, y);
840 z = (*xa++ >> 16) + borrow;
841 borrow = z >> 16;
842 Sign_Extend(borrow, z);
843 Storeinc(xc, z, y);
844 }
845 #else
846 do {
847 y = *xa++ - *xb++ + borrow;
848 borrow = y >> 16;
849 Sign_Extend(borrow, y);
850 *xc++ = y & 0xffff;
851 } while (xb < xbe);
852 while (xa < xae) {
853 y = *xa++ + borrow;
854 borrow = y >> 16;
855 Sign_Extend(borrow, y);
856 *xc++ = y & 0xffff;
857 }
858 #endif
859 while (!*--xc)
860 wa--;
861 c->wds = wa;
862 return c;
863 }
864
865 static double
ulp(x)866 ulp
867 #ifdef KR_headers
868 (x) double x;
869 #else
870 (double x)
871 #endif
872 {
873 register long L;
874 double a;
875
876 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
877 #ifndef Sudden_Underflow
878 if (L > 0) {
879 #endif
880 #ifdef IBM
881 L |= Exp_msk1 >> 4;
882 #endif
883 word0(a) = L;
884 word1(a) = 0;
885 #ifndef Sudden_Underflow
886 } else {
887 L = -L >> Exp_shift;
888 if (L < Exp_shift) {
889 word0(a) = 0x80000 >> L;
890 word1(a) = 0;
891 } else {
892 word0(a) = 0;
893 L -= Exp_shift;
894 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
895 }
896 }
897 #endif
898 return a;
899 }
900
901 static double
b2d(a,e)902 b2d
903 #ifdef KR_headers
904 (a, e) Bigint *a; int *e;
905 #else
906 (Bigint *a, int *e)
907 #endif
908 {
909 unsigned long *xa, *xa0, w, y, z;
910 int k;
911 double d;
912 #ifdef VAX
913 unsigned long d0, d1;
914 #else
915 #define d0 word0(d)
916 #define d1 word1(d)
917 #endif
918
919 xa0 = a->x;
920 xa = xa0 + a->wds;
921 y = *--xa;
922 #ifdef DEBUG
923 if (!y) Bug("zero y in b2d");
924 #endif
925 k = hi0bits(y);
926 *e = 32 - k;
927 #ifdef Pack_32
928 if (k < Ebits) {
929 d0 = Exp_1 | y >> Ebits - k;
930 w = xa > xa0 ? *--xa : 0;
931 d1 = y << (32-Ebits) + k | w >> Ebits - k;
932 goto ret_d;
933 }
934 z = xa > xa0 ? *--xa : 0;
935 if (k -= Ebits) {
936 d0 = Exp_1 | y << k | z >> 32 - k;
937 y = xa > xa0 ? *--xa : 0;
938 d1 = z << k | y >> 32 - k;
939 } else {
940 d0 = Exp_1 | y;
941 d1 = z;
942 }
943 #else
944 if (k < Ebits + 16) {
945 z = xa > xa0 ? *--xa : 0;
946 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
947 w = xa > xa0 ? *--xa : 0;
948 y = xa > xa0 ? *--xa : 0;
949 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
950 goto ret_d;
951 }
952 z = xa > xa0 ? *--xa : 0;
953 w = xa > xa0 ? *--xa : 0;
954 k -= Ebits + 16;
955 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
956 y = xa > xa0 ? *--xa : 0;
957 d1 = w << k + 16 | y << k;
958 #endif
959 ret_d:
960 #ifdef VAX
961 word0(d) = d0 >> 16 | d0 << 16;
962 word1(d) = d1 >> 16 | d1 << 16;
963 #else
964 #undef d0
965 #undef d1
966 #endif
967 return d;
968 }
969
970 static Bigint *
d2b(d,e,bits)971 d2b
972 #ifdef KR_headers
973 (d, e, bits) double d; int *e, *bits;
974 #else
975 (double d, int *e, int *bits)
976 #endif
977 {
978 Bigint *b;
979 int de, i, k;
980 unsigned long *x, y, z;
981 #ifdef VAX
982 unsigned long d0, d1;
983 d0 = word0(d) >> 16 | word0(d) << 16;
984 d1 = word1(d) >> 16 | word1(d) << 16;
985 #else
986 #define d0 word0(d)
987 #define d1 word1(d)
988 #endif
989
990 #ifdef Pack_32
991 b = Balloc(1);
992 #else
993 b = Balloc(2);
994 #endif
995 x = b->x;
996
997 z = d0 & Frac_mask;
998 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
999 #ifdef Sudden_Underflow
1000 de = (int)(d0 >> Exp_shift);
1001 #ifndef IBM
1002 z |= Exp_msk11;
1003 #endif
1004 #else
1005 if (de = (int)(d0 >> Exp_shift))
1006 z |= Exp_msk1;
1007 #endif
1008 #ifdef Pack_32
1009 if (y = d1) {
1010 if (k = lo0bits(&y)) {
1011 x[0] = y | z << 32 - k;
1012 z >>= k;
1013 }
1014 else
1015 x[0] = y;
1016 i = b->wds = (x[1] = z) ? 2 : 1;
1017 } else {
1018 #ifdef DEBUG
1019 if (!z)
1020 Bug("Zero passed to d2b");
1021 #endif
1022 k = lo0bits(&z);
1023 x[0] = z;
1024 i = b->wds = 1;
1025 k += 32;
1026 }
1027 #else
1028 if (y = d1) {
1029 if (k = lo0bits(&y))
1030 if (k >= 16) {
1031 x[0] = y | z << 32 - k & 0xffff;
1032 x[1] = z >> k - 16 & 0xffff;
1033 x[2] = z >> k;
1034 i = 2;
1035 } else {
1036 x[0] = y & 0xffff;
1037 x[1] = y >> 16 | z << 16 - k & 0xffff;
1038 x[2] = z >> k & 0xffff;
1039 x[3] = z >> k+16;
1040 i = 3;
1041 }
1042 else {
1043 x[0] = y & 0xffff;
1044 x[1] = y >> 16;
1045 x[2] = z & 0xffff;
1046 x[3] = z >> 16;
1047 i = 3;
1048 }
1049 } else {
1050 #ifdef DEBUG
1051 if (!z)
1052 Bug("Zero passed to d2b");
1053 #endif
1054 k = lo0bits(&z);
1055 if (k >= 16) {
1056 x[0] = z;
1057 i = 0;
1058 } else {
1059 x[0] = z & 0xffff;
1060 x[1] = z >> 16;
1061 i = 1;
1062 }
1063 k += 32;
1064 }
1065 while (!x[i])
1066 --i;
1067 b->wds = i + 1;
1068 #endif
1069 #ifndef Sudden_Underflow
1070 if (de) {
1071 #endif
1072 #ifdef IBM
1073 *e = (de - Bias - (P-1) << 2) + k;
1074 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1075 #else
1076 *e = de - Bias - (P-1) + k;
1077 *bits = P - k;
1078 #endif
1079 #ifndef Sudden_Underflow
1080 } else {
1081 *e = de - Bias - (P-1) + 1 + k;
1082 #ifdef Pack_32
1083 *bits = 32*i - hi0bits(x[i-1]);
1084 #else
1085 *bits = (i+2)*16 - hi0bits(x[i]);
1086 #endif
1087 }
1088 #endif
1089 return b;
1090 }
1091 #undef d0
1092 #undef d1
1093
1094 static double
ratio(a,b)1095 ratio
1096 #ifdef KR_headers
1097 (a, b) Bigint *a, *b;
1098 #else
1099 (Bigint *a, Bigint *b)
1100 #endif
1101 {
1102 double da, db;
1103 int k, ka, kb;
1104
1105 da = b2d(a, &ka);
1106 db = b2d(b, &kb);
1107 #ifdef Pack_32
1108 k = ka - kb + 32*(a->wds - b->wds);
1109 #else
1110 k = ka - kb + 16*(a->wds - b->wds);
1111 #endif
1112 #ifdef IBM
1113 if (k > 0) {
1114 word0(da) += (k >> 2)*Exp_msk1;
1115 if (k &= 3)
1116 da *= 1 << k;
1117 } else {
1118 k = -k;
1119 word0(db) += (k >> 2)*Exp_msk1;
1120 if (k &= 3)
1121 db *= 1 << k;
1122 }
1123 #else
1124 if (k > 0)
1125 word0(da) += k*Exp_msk1;
1126 else {
1127 k = -k;
1128 word0(db) += k*Exp_msk1;
1129 }
1130 #endif
1131 return da / db;
1132 }
1133
1134 static double
1135 tens[] = {
1136 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1137 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1138 1e20, 1e21, 1e22
1139 #ifdef VAX
1140 , 1e23, 1e24
1141 #endif
1142 };
1143
1144 static double
1145 #ifdef IEEE_Arith
1146 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1147 static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1148 #define n_bigtens 5
1149 #else
1150 #ifdef IBM
1151 bigtens[] = { 1e16, 1e32, 1e64 };
1152 static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1153 #define n_bigtens 3
1154 #else
1155 bigtens[] = { 1e16, 1e32 };
1156 static double tinytens[] = { 1e-16, 1e-32 };
1157 #define n_bigtens 2
1158 #endif
1159 #endif
1160
1161 double
strtod(s00,se)1162 strtod
1163 #ifdef KR_headers
1164 (s00, se) CONST char *s00; char **se;
1165 #else
1166 (CONST char *s00, char **se)
1167 #endif
1168 {
1169 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1170 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1171 CONST char *s, *s0, *s1;
1172 double aadj, aadj1, adj, rv, rv0;
1173 long L;
1174 unsigned long y, z;
1175 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1176 sign = nz0 = nz = 0;
1177 rv = 0.;
1178 for (s = s00;;s++) switch(*s) {
1179 case '-':
1180 sign = 1;
1181 /* no break */
1182 case '+':
1183 if (*++s)
1184 goto break2;
1185 /* no break */
1186 case 0:
1187 s = s00;
1188 goto ret;
1189 case '\t':
1190 case '\n':
1191 case '\v':
1192 case '\f':
1193 case '\r':
1194 case ' ':
1195 continue;
1196 default:
1197 goto break2;
1198 }
1199 break2:
1200 if (*s == '0') {
1201 nz0 = 1;
1202 while (*++s == '0') ;
1203 if (!*s)
1204 goto ret;
1205 }
1206 s0 = s;
1207 y = z = 0;
1208 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1209 if (nd < 9)
1210 y = 10*y + c - '0';
1211 else if (nd < 16)
1212 z = 10*z + c - '0';
1213 nd0 = nd;
1214 if (c == '.') {
1215 c = *++s;
1216 if (!nd) {
1217 for (; c == '0'; c = *++s)
1218 nz++;
1219 if (c > '0' && c <= '9') {
1220 s0 = s;
1221 nf += nz;
1222 nz = 0;
1223 goto have_dig;
1224 }
1225 goto dig_done;
1226 }
1227 for (; c >= '0' && c <= '9'; c = *++s) {
1228 have_dig:
1229 nz++;
1230 if (c -= '0') {
1231 nf += nz;
1232 for (i = 1; i < nz; i++)
1233 if (nd++ < 9)
1234 y *= 10;
1235 else if (nd <= DBL_DIG + 1)
1236 z *= 10;
1237 if (nd++ < 9)
1238 y = 10*y + c;
1239 else if (nd <= DBL_DIG + 1)
1240 z = 10*z + c;
1241 nz = 0;
1242 }
1243 }
1244 }
1245 dig_done:
1246 e = 0;
1247 if (c == 'e' || c == 'E') {
1248 if (!nd && !nz && !nz0) {
1249 s = s00;
1250 goto ret;
1251 }
1252 s00 = s;
1253 esign = 0;
1254 switch(c = *++s) {
1255 case '-':
1256 esign = 1;
1257 case '+':
1258 c = *++s;
1259 }
1260 if (c >= '0' && c <= '9') {
1261 while (c == '0')
1262 c = *++s;
1263 if (c > '0' && c <= '9') {
1264 L = c - '0';
1265 s1 = s;
1266 while ((c = *++s) >= '0' && c <= '9')
1267 L = 10*L + c - '0';
1268 if (s - s1 > 8 || L > 19999)
1269 /* Avoid confusion from exponents
1270 * so large that e might overflow.
1271 */
1272 e = 19999; /* safe for 16 bit ints */
1273 else
1274 e = (int)L;
1275 if (esign)
1276 e = -e;
1277 } else
1278 e = 0;
1279 } else
1280 s = s00;
1281 }
1282 if (!nd) {
1283 if (!nz && !nz0)
1284 s = s00;
1285 goto ret;
1286 }
1287 e1 = e -= nf;
1288
1289 /* Now we have nd0 digits, starting at s0, followed by a
1290 * decimal point, followed by nd-nd0 digits. The number we're
1291 * after is the integer represented by those digits times
1292 * 10**e */
1293
1294 if (!nd0)
1295 nd0 = nd;
1296 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1297 rv = y;
1298 if (k > 9)
1299 rv = tens[k - 9] * rv + z;
1300 if (nd <= DBL_DIG
1301 #ifndef RND_PRODQUOT
1302 && FLT_ROUNDS == 1
1303 #endif
1304 ) {
1305 if (!e)
1306 goto ret;
1307 if (e > 0) {
1308 if (e <= Ten_pmax) {
1309 #ifdef VAX
1310 goto vax_ovfl_check;
1311 #else
1312 /* rv = */ rounded_product(rv, tens[e]);
1313 goto ret;
1314 #endif
1315 }
1316 i = DBL_DIG - nd;
1317 if (e <= Ten_pmax + i) {
1318 /* A fancier test would sometimes let us do
1319 * this for larger i values.
1320 */
1321 e -= i;
1322 rv *= tens[i];
1323 #ifdef VAX
1324 /* VAX exponent range is so narrow we must
1325 * worry about overflow here...
1326 */
1327 vax_ovfl_check:
1328 word0(rv) -= P*Exp_msk1;
1329 /* rv = */ rounded_product(rv, tens[e]);
1330 if ((word0(rv) & Exp_mask)
1331 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1332 goto ovfl;
1333 word0(rv) += P*Exp_msk1;
1334 #else
1335 /* rv = */ rounded_product(rv, tens[e]);
1336 #endif
1337 goto ret;
1338 }
1339 }
1340 #ifndef Inaccurate_Divide
1341 else if (e >= -Ten_pmax) {
1342 /* rv = */ rounded_quotient(rv, tens[-e]);
1343 goto ret;
1344 }
1345 #endif
1346 }
1347 e1 += nd - k;
1348
1349 /* Get starting approximation = rv * 10**e1 */
1350
1351 if (e1 > 0) {
1352 if (i = e1 & 15)
1353 rv *= tens[i];
1354 if (e1 &= ~15) {
1355 if (e1 > DBL_MAX_10_EXP) {
1356 ovfl:
1357 errno = ERANGE;
1358 #ifdef __STDC__
1359 rv = HUGE_VAL;
1360 #else
1361 /* Can't trust HUGE_VAL */
1362 #ifdef IEEE_Arith
1363 word0(rv) = Exp_mask;
1364 word1(rv) = 0;
1365 #else
1366 word0(rv) = Big0;
1367 word1(rv) = Big1;
1368 #endif
1369 #endif
1370 goto ret;
1371 }
1372 if (e1 >>= 4) {
1373 for (j = 0; e1 > 1; j++, e1 >>= 1)
1374 if (e1 & 1)
1375 rv *= bigtens[j];
1376 /* The last multiplication could overflow. */
1377 word0(rv) -= P*Exp_msk1;
1378 rv *= bigtens[j];
1379 if ((z = word0(rv) & Exp_mask)
1380 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1381 goto ovfl;
1382 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1383 /* set to largest number */
1384 /* (Can't trust DBL_MAX) */
1385 word0(rv) = Big0;
1386 word1(rv) = Big1;
1387 }
1388 else
1389 word0(rv) += P*Exp_msk1;
1390 }
1391 }
1392 } else if (e1 < 0) {
1393 e1 = -e1;
1394 if (i = e1 & 15)
1395 rv /= tens[i];
1396 if (e1 &= ~15) {
1397 e1 >>= 4;
1398 for (j = 0; e1 > 1; j++, e1 >>= 1)
1399 if (e1 & 1)
1400 rv *= tinytens[j];
1401 /* The last multiplication could underflow. */
1402 rv0 = rv;
1403 rv *= tinytens[j];
1404 if (!rv) {
1405 rv = 2.*rv0;
1406 rv *= tinytens[j];
1407 if (!rv) {
1408 undfl:
1409 rv = 0.;
1410 errno = ERANGE;
1411 goto ret;
1412 }
1413 word0(rv) = Tiny0;
1414 word1(rv) = Tiny1;
1415 /* The refinement below will clean
1416 * this approximation up.
1417 */
1418 }
1419 }
1420 }
1421
1422 /* Now the hard part -- adjusting rv to the correct value.*/
1423
1424 /* Put digits into bd: true value = bd * 10^e */
1425
1426 bd0 = s2b(s0, nd0, nd, y);
1427
1428 for (;;) {
1429 bd = Balloc(bd0->k);
1430 Bcopy(bd, bd0);
1431 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1432 bs = i2b(1);
1433
1434 if (e >= 0) {
1435 bb2 = bb5 = 0;
1436 bd2 = bd5 = e;
1437 } else {
1438 bb2 = bb5 = -e;
1439 bd2 = bd5 = 0;
1440 }
1441 if (bbe >= 0)
1442 bb2 += bbe;
1443 else
1444 bd2 -= bbe;
1445 bs2 = bb2;
1446 #ifdef Sudden_Underflow
1447 #ifdef IBM
1448 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1449 #else
1450 j = P + 1 - bbbits;
1451 #endif
1452 #else
1453 i = bbe + bbbits - 1; /* logb(rv) */
1454 if (i < Emin) /* denormal */
1455 j = bbe + (P-Emin);
1456 else
1457 j = P + 1 - bbbits;
1458 #endif
1459 bb2 += j;
1460 bd2 += j;
1461 i = bb2 < bd2 ? bb2 : bd2;
1462 if (i > bs2)
1463 i = bs2;
1464 if (i > 0) {
1465 bb2 -= i;
1466 bd2 -= i;
1467 bs2 -= i;
1468 }
1469 if (bb5 > 0) {
1470 bs = pow5mult(bs, bb5);
1471 bb1 = mult(bs, bb);
1472 Bfree(bb);
1473 bb = bb1;
1474 }
1475 if (bb2 > 0)
1476 bb = lshift(bb, bb2);
1477 if (bd5 > 0)
1478 bd = pow5mult(bd, bd5);
1479 if (bd2 > 0)
1480 bd = lshift(bd, bd2);
1481 if (bs2 > 0)
1482 bs = lshift(bs, bs2);
1483 delta = diff(bb, bd);
1484 dsign = delta->sign;
1485 delta->sign = 0;
1486 i = cmp(delta, bs);
1487 if (i < 0) {
1488 /* Error is less than half an ulp -- check for
1489 * special case of mantissa a power of two.
1490 */
1491 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1492 break;
1493 delta = lshift(delta,Log2P);
1494 if (cmp(delta, bs) > 0)
1495 goto drop_down;
1496 break;
1497 }
1498 if (i == 0) {
1499 /* exactly half-way between */
1500 if (dsign) {
1501 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1502 && word1(rv) == 0xffffffff) {
1503 /*boundary case -- increment exponent*/
1504 word0(rv) = (word0(rv) & Exp_mask)
1505 + Exp_msk1
1506 #ifdef IBM
1507 | Exp_msk1 >> 4
1508 #endif
1509 ;
1510 word1(rv) = 0;
1511 break;
1512 }
1513 } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1514 drop_down:
1515 /* boundary case -- decrement exponent */
1516 #ifdef Sudden_Underflow
1517 L = word0(rv) & Exp_mask;
1518 #ifdef IBM
1519 if (L < Exp_msk1)
1520 #else
1521 if (L <= Exp_msk1)
1522 #endif
1523 goto undfl;
1524 L -= Exp_msk1;
1525 #else
1526 L = (word0(rv) & Exp_mask) - Exp_msk1;
1527 #endif
1528 word0(rv) = L | Bndry_mask1;
1529 word1(rv) = 0xffffffff;
1530 #ifdef IBM
1531 goto cont;
1532 #else
1533 break;
1534 #endif
1535 }
1536 #ifndef ROUND_BIASED
1537 if (!(word1(rv) & LSB))
1538 break;
1539 #endif
1540 if (dsign)
1541 rv += ulp(rv);
1542 #ifndef ROUND_BIASED
1543 else {
1544 rv -= ulp(rv);
1545 #ifndef Sudden_Underflow
1546 if (!rv)
1547 goto undfl;
1548 #endif
1549 }
1550 #endif
1551 break;
1552 }
1553 if ((aadj = ratio(delta, bs)) <= 2.) {
1554 if (dsign)
1555 aadj = aadj1 = 1.;
1556 else if (word1(rv) || word0(rv) & Bndry_mask) {
1557 #ifndef Sudden_Underflow
1558 if (word1(rv) == Tiny1 && !word0(rv))
1559 goto undfl;
1560 #endif
1561 aadj = 1.;
1562 aadj1 = -1.;
1563 } else {
1564 /* special case -- power of FLT_RADIX to be */
1565 /* rounded down... */
1566
1567 if (aadj < 2./FLT_RADIX)
1568 aadj = 1./FLT_RADIX;
1569 else
1570 aadj *= 0.5;
1571 aadj1 = -aadj;
1572 }
1573 } else {
1574 aadj *= 0.5;
1575 aadj1 = dsign ? aadj : -aadj;
1576 #ifdef Check_FLT_ROUNDS
1577 switch(FLT_ROUNDS) {
1578 case 2: /* towards +infinity */
1579 aadj1 -= 0.5;
1580 break;
1581 case 0: /* towards 0 */
1582 case 3: /* towards -infinity */
1583 aadj1 += 0.5;
1584 }
1585 #else
1586 if (FLT_ROUNDS == 0)
1587 aadj1 += 0.5;
1588 #endif
1589 }
1590 y = word0(rv) & Exp_mask;
1591
1592 /* Check for overflow */
1593
1594 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1595 rv0 = rv;
1596 word0(rv) -= P*Exp_msk1;
1597 adj = aadj1 * ulp(rv);
1598 rv += adj;
1599 if ((word0(rv) & Exp_mask) >=
1600 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1601 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1602 goto ovfl;
1603 word0(rv) = Big0;
1604 word1(rv) = Big1;
1605 goto cont;
1606 } else
1607 word0(rv) += P*Exp_msk1;
1608 } else {
1609 #ifdef Sudden_Underflow
1610 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1611 rv0 = rv;
1612 word0(rv) += P*Exp_msk1;
1613 adj = aadj1 * ulp(rv);
1614 rv += adj;
1615 #ifdef IBM
1616 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1617 #else
1618 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1619 #endif
1620 {
1621 if (word0(rv0) == Tiny0
1622 && word1(rv0) == Tiny1)
1623 goto undfl;
1624 word0(rv) = Tiny0;
1625 word1(rv) = Tiny1;
1626 goto cont;
1627 } else
1628 word0(rv) -= P*Exp_msk1;
1629 } else {
1630 adj = aadj1 * ulp(rv);
1631 rv += adj;
1632 }
1633 #else
1634 /* Compute adj so that the IEEE rounding rules will
1635 * correctly round rv + adj in some half-way cases.
1636 * If rv * ulp(rv) is denormalized (i.e.,
1637 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1638 * trouble from bits lost to denormalization;
1639 * example: 1.2e-307 .
1640 */
1641 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1642 aadj1 = (double)(int)(aadj + 0.5);
1643 if (!dsign)
1644 aadj1 = -aadj1;
1645 }
1646 adj = aadj1 * ulp(rv);
1647 rv += adj;
1648 #endif
1649 }
1650 z = word0(rv) & Exp_mask;
1651 if (y == z) {
1652 /* Can we stop now? */
1653 L = aadj;
1654 aadj -= L;
1655 /* The tolerances below are conservative. */
1656 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1657 if (aadj < .4999999 || aadj > .5000001)
1658 break;
1659 } else if (aadj < .4999999/FLT_RADIX)
1660 break;
1661 }
1662 cont:
1663 Bfree(bb);
1664 Bfree(bd);
1665 Bfree(bs);
1666 Bfree(delta);
1667 }
1668 Bfree(bb);
1669 Bfree(bd);
1670 Bfree(bs);
1671 Bfree(bd0);
1672 Bfree(delta);
1673 ret:
1674 if (se)
1675 *se = (char *)s;
1676 return sign ? -rv : rv;
1677 }
1678
1679 static int
quorem(b,S)1680 quorem
1681 #ifdef KR_headers
1682 (b, S) Bigint *b, *S;
1683 #else
1684 (Bigint *b, Bigint *S)
1685 #endif
1686 {
1687 int n;
1688 long borrow, y;
1689 unsigned long carry, q, ys;
1690 unsigned long *bx, *bxe, *sx, *sxe;
1691 #ifdef Pack_32
1692 long z;
1693 unsigned long si, zs;
1694 #endif
1695
1696 n = S->wds;
1697 #ifdef DEBUG
1698 /*debug*/ if (b->wds > n)
1699 /*debug*/ Bug("oversize b in quorem");
1700 #endif
1701 if (b->wds < n)
1702 return 0;
1703 sx = S->x;
1704 sxe = sx + --n;
1705 bx = b->x;
1706 bxe = bx + n;
1707 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1708 #ifdef DEBUG
1709 /*debug*/ if (q > 9)
1710 /*debug*/ Bug("oversized quotient in quorem");
1711 #endif
1712 if (q) {
1713 borrow = 0;
1714 carry = 0;
1715 do {
1716 #ifdef Pack_32
1717 si = *sx++;
1718 ys = (si & 0xffff) * q + carry;
1719 zs = (si >> 16) * q + (ys >> 16);
1720 carry = zs >> 16;
1721 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1722 borrow = y >> 16;
1723 Sign_Extend(borrow, y);
1724 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1725 borrow = z >> 16;
1726 Sign_Extend(borrow, z);
1727 Storeinc(bx, z, y);
1728 #else
1729 ys = *sx++ * q + carry;
1730 carry = ys >> 16;
1731 y = *bx - (ys & 0xffff) + borrow;
1732 borrow = y >> 16;
1733 Sign_Extend(borrow, y);
1734 *bx++ = y & 0xffff;
1735 #endif
1736 } while (sx <= sxe);
1737 if (!*bxe) {
1738 bx = b->x;
1739 while (--bxe > bx && !*bxe)
1740 --n;
1741 b->wds = n;
1742 }
1743 }
1744 if (cmp(b, S) >= 0) {
1745 q++;
1746 borrow = 0;
1747 carry = 0;
1748 bx = b->x;
1749 sx = S->x;
1750 do {
1751 #ifdef Pack_32
1752 si = *sx++;
1753 ys = (si & 0xffff) + carry;
1754 zs = (si >> 16) + (ys >> 16);
1755 carry = zs >> 16;
1756 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1757 borrow = y >> 16;
1758 Sign_Extend(borrow, y);
1759 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1760 borrow = z >> 16;
1761 Sign_Extend(borrow, z);
1762 Storeinc(bx, z, y);
1763 #else
1764 ys = *sx++ + carry;
1765 carry = ys >> 16;
1766 y = *bx - (ys & 0xffff) + borrow;
1767 borrow = y >> 16;
1768 Sign_Extend(borrow, y);
1769 *bx++ = y & 0xffff;
1770 #endif
1771 } while (sx <= sxe);
1772 bx = b->x;
1773 bxe = bx + n;
1774 if (!*bxe) {
1775 while (--bxe > bx && !*bxe)
1776 --n;
1777 b->wds = n;
1778 }
1779 }
1780 return q;
1781 }
1782
1783 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1784 *
1785 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1786 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1787 *
1788 * Modifications:
1789 * 1. Rather than iterating, we use a simple numeric overestimate
1790 * to determine k = floor(log10(d)). We scale relevant
1791 * quantities using O(log2(k)) rather than O(k) multiplications.
1792 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1793 * try to generate digits strictly left to right. Instead, we
1794 * compute with fewer bits and propagate the carry if necessary
1795 * when rounding the final digit up. This is often faster.
1796 * 3. Under the assumption that input will be rounded nearest,
1797 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1798 * That is, we allow equality in stopping tests when the
1799 * round-nearest rule will give the same floating-point value
1800 * as would satisfaction of the stopping test with strict
1801 * inequality.
1802 * 4. We remove common factors of powers of 2 from relevant
1803 * quantities.
1804 * 5. When converting floating-point integers less than 1e16,
1805 * we use floating-point arithmetic rather than resorting
1806 * to multiple-precision integers.
1807 * 6. When asked to produce fewer than 15 digits, we first try
1808 * to get by with floating-point arithmetic; we resort to
1809 * multiple-precision integer arithmetic only if we cannot
1810 * guarantee that the floating-point calculation has given
1811 * the correctly rounded result. For k requested digits and
1812 * "uniformly" distributed input, the probability is
1813 * something like 10^(k-15) that we must resort to the long
1814 * calculation.
1815 */
1816
1817 char *
__dtoa(d,mode,ndigits,decpt,sign,rve)1818 __dtoa
1819 #ifdef KR_headers
1820 (d, mode, ndigits, decpt, sign, rve)
1821 double d; int mode, ndigits, *decpt, *sign; char **rve;
1822 #else
1823 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1824 #endif
1825 {
1826 /* Arguments ndigits, decpt, sign are similar to those
1827 of ecvt and fcvt; trailing zeros are suppressed from
1828 the returned string. If not null, *rve is set to point
1829 to the end of the return value. If d is +-Infinity or NaN,
1830 then *decpt is set to 9999.
1831
1832 mode:
1833 0 ==> shortest string that yields d when read in
1834 and rounded to nearest.
1835 1 ==> like 0, but with Steele & White stopping rule;
1836 e.g. with IEEE P754 arithmetic , mode 0 gives
1837 1e23 whereas mode 1 gives 9.999999999999999e22.
1838 2 ==> max(1,ndigits) significant digits. This gives a
1839 return value similar to that of ecvt, except
1840 that trailing zeros are suppressed.
1841 3 ==> through ndigits past the decimal point. This
1842 gives a return value similar to that from fcvt,
1843 except that trailing zeros are suppressed, and
1844 ndigits can be negative.
1845 4-9 should give the same return values as 2-3, i.e.,
1846 4 <= mode <= 9 ==> same return as mode
1847 2 + (mode & 1). These modes are mainly for
1848 debugging; often they run slower but sometimes
1849 faster than modes 2-3.
1850 4,5,8,9 ==> left-to-right digit generation.
1851 6-9 ==> don't try fast floating-point estimate
1852 (if applicable).
1853
1854 Values of mode other than 0-9 are treated as mode 0.
1855
1856 Sufficient space is allocated to the return value
1857 to hold the suppressed trailing zeros.
1858 */
1859
1860 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1861 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1862 spec_case, try_quick;
1863 long L;
1864 #ifndef Sudden_Underflow
1865 int denorm;
1866 unsigned long x;
1867 #endif
1868 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1869 double d2, ds, eps;
1870 char *s, *s0;
1871 static Bigint *result;
1872 static int result_k;
1873
1874 if (result) {
1875 result->k = result_k;
1876 result->maxwds = 1 << result_k;
1877 Bfree(result);
1878 result = 0;
1879 }
1880
1881 if (word0(d) & Sign_bit) {
1882 /* set sign for everything, including 0's and NaNs */
1883 *sign = 1;
1884 word0(d) &= ~Sign_bit; /* clear sign bit */
1885 }
1886 else
1887 *sign = 0;
1888
1889 #if defined(IEEE_Arith) + defined(VAX)
1890 #ifdef IEEE_Arith
1891 if ((word0(d) & Exp_mask) == Exp_mask)
1892 #else
1893 if (word0(d) == 0x8000)
1894 #endif
1895 {
1896 /* Infinity or NaN */
1897 *decpt = 9999;
1898 s =
1899 #ifdef IEEE_Arith
1900 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1901 #endif
1902 "NaN";
1903 if (rve)
1904 *rve =
1905 #ifdef IEEE_Arith
1906 s[3] ? s + 8 :
1907 #endif
1908 s + 3;
1909 return s;
1910 }
1911 #endif
1912 #ifdef IBM
1913 d += 0; /* normalize */
1914 #endif
1915 if (!d) {
1916 *decpt = 1;
1917 s = "0";
1918 if (rve)
1919 *rve = s + 1;
1920 return s;
1921 }
1922
1923 b = d2b(d, &be, &bbits);
1924 #ifdef Sudden_Underflow
1925 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1926 #else
1927 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
1928 #endif
1929 d2 = d;
1930 word0(d2) &= Frac_mask1;
1931 word0(d2) |= Exp_11;
1932 #ifdef IBM
1933 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
1934 d2 /= 1 << j;
1935 #endif
1936
1937 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1938 * log10(x) = log(x) / log(10)
1939 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1940 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1941 *
1942 * This suggests computing an approximation k to log10(d) by
1943 *
1944 * k = (i - Bias)*0.301029995663981
1945 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1946 *
1947 * We want k to be too large rather than too small.
1948 * The error in the first-order Taylor series approximation
1949 * is in our favor, so we just round up the constant enough
1950 * to compensate for any error in the multiplication of
1951 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1952 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1953 * adding 1e-13 to the constant term more than suffices.
1954 * Hence we adjust the constant term to 0.1760912590558.
1955 * (We could get a more accurate k by invoking log10,
1956 * but this is probably not worthwhile.)
1957 */
1958
1959 i -= Bias;
1960 #ifdef IBM
1961 i <<= 2;
1962 i += j;
1963 #endif
1964 #ifndef Sudden_Underflow
1965 denorm = 0;
1966 } else {
1967 /* d is denormalized */
1968
1969 i = bbits + be + (Bias + (P-1) - 1);
1970 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
1971 : word1(d) << 32 - i;
1972 d2 = x;
1973 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
1974 i -= (Bias + (P-1) - 1) + 1;
1975 denorm = 1;
1976 }
1977 #endif
1978 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
1979 k = (int)ds;
1980 if (ds < 0. && ds != k)
1981 k--; /* want k = floor(ds) */
1982 k_check = 1;
1983 if (k >= 0 && k <= Ten_pmax) {
1984 if (d < tens[k])
1985 k--;
1986 k_check = 0;
1987 }
1988 j = bbits - i - 1;
1989 if (j >= 0) {
1990 b2 = 0;
1991 s2 = j;
1992 } else {
1993 b2 = -j;
1994 s2 = 0;
1995 }
1996 if (k >= 0) {
1997 b5 = 0;
1998 s5 = k;
1999 s2 += k;
2000 } else {
2001 b2 -= k;
2002 b5 = -k;
2003 s5 = 0;
2004 }
2005 if (mode < 0 || mode > 9)
2006 mode = 0;
2007 try_quick = 1;
2008 if (mode > 5) {
2009 mode -= 4;
2010 try_quick = 0;
2011 }
2012 leftright = 1;
2013 switch(mode) {
2014 case 0:
2015 case 1:
2016 ilim = ilim1 = -1;
2017 i = 18;
2018 ndigits = 0;
2019 break;
2020 case 2:
2021 leftright = 0;
2022 /* no break */
2023 case 4:
2024 if (ndigits <= 0)
2025 ndigits = 1;
2026 ilim = ilim1 = i = ndigits;
2027 break;
2028 case 3:
2029 leftright = 0;
2030 /* no break */
2031 case 5:
2032 i = ndigits + k + 1;
2033 ilim = i;
2034 ilim1 = i - 1;
2035 if (i <= 0)
2036 i = 1;
2037 }
2038 j = sizeof(unsigned long);
2039 for (result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i;
2040 j <<= 1) result_k++;
2041 result = Balloc(result_k);
2042 s = s0 = (char *)result;
2043
2044 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2045
2046 /* Try to get by with floating-point arithmetic. */
2047
2048 i = 0;
2049 d2 = d;
2050 k0 = k;
2051 ilim0 = ilim;
2052 ieps = 2; /* conservative */
2053 if (k > 0) {
2054 ds = tens[k&0xf];
2055 j = k >> 4;
2056 if (j & Bletch) {
2057 /* prevent overflows */
2058 j &= Bletch - 1;
2059 d /= bigtens[n_bigtens-1];
2060 ieps++;
2061 }
2062 for (; j; j >>= 1, i++)
2063 if (j & 1) {
2064 ieps++;
2065 ds *= bigtens[i];
2066 }
2067 d /= ds;
2068 } else if (j1 = -k) {
2069 d *= tens[j1 & 0xf];
2070 for (j = j1 >> 4; j; j >>= 1, i++)
2071 if (j & 1) {
2072 ieps++;
2073 d *= bigtens[i];
2074 }
2075 }
2076 if (k_check && d < 1. && ilim > 0) {
2077 if (ilim1 <= 0)
2078 goto fast_failed;
2079 ilim = ilim1;
2080 k--;
2081 d *= 10.;
2082 ieps++;
2083 }
2084 eps = ieps*d + 7.;
2085 word0(eps) -= (P-1)*Exp_msk1;
2086 if (ilim == 0) {
2087 S = mhi = 0;
2088 d -= 5.;
2089 if (d > eps)
2090 goto one_digit;
2091 if (d < -eps)
2092 goto no_digits;
2093 goto fast_failed;
2094 }
2095 #ifndef No_leftright
2096 if (leftright) {
2097 /* Use Steele & White method of only
2098 * generating digits needed.
2099 */
2100 eps = 0.5/tens[ilim-1] - eps;
2101 for (i = 0;;) {
2102 L = d;
2103 d -= L;
2104 *s++ = '0' + (int)L;
2105 if (d < eps)
2106 goto ret1;
2107 if (1. - d < eps)
2108 goto bump_up;
2109 if (++i >= ilim)
2110 break;
2111 eps *= 10.;
2112 d *= 10.;
2113 }
2114 } else {
2115 #endif
2116 /* Generate ilim digits, then fix them up. */
2117 eps *= tens[ilim-1];
2118 for (i = 1;; i++, d *= 10.) {
2119 L = d;
2120 d -= L;
2121 *s++ = '0' + (int)L;
2122 if (i == ilim) {
2123 if (d > 0.5 + eps)
2124 goto bump_up;
2125 else if (d < 0.5 - eps) {
2126 while (*--s == '0');
2127 s++;
2128 goto ret1;
2129 }
2130 break;
2131 }
2132 }
2133 #ifndef No_leftright
2134 }
2135 #endif
2136 fast_failed:
2137 s = s0;
2138 d = d2;
2139 k = k0;
2140 ilim = ilim0;
2141 }
2142
2143 /* Do we have a "small" integer? */
2144
2145 if (be >= 0 && k <= Int_max) {
2146 /* Yes. */
2147 ds = tens[k];
2148 if (ndigits < 0 && ilim <= 0) {
2149 S = mhi = 0;
2150 if (ilim < 0 || d <= 5*ds)
2151 goto no_digits;
2152 goto one_digit;
2153 }
2154 for (i = 1;; i++) {
2155 L = d / ds;
2156 d -= L*ds;
2157 #ifdef Check_FLT_ROUNDS
2158 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2159 if (d < 0) {
2160 L--;
2161 d += ds;
2162 }
2163 #endif
2164 *s++ = '0' + (int)L;
2165 if (i == ilim) {
2166 d += d;
2167 if (d > ds || d == ds && L & 1) {
2168 bump_up:
2169 while (*--s == '9')
2170 if (s == s0) {
2171 k++;
2172 *s = '0';
2173 break;
2174 }
2175 ++*s++;
2176 }
2177 break;
2178 }
2179 if (!(d *= 10.))
2180 break;
2181 }
2182 goto ret1;
2183 }
2184
2185 m2 = b2;
2186 m5 = b5;
2187 mhi = mlo = 0;
2188 if (leftright) {
2189 if (mode < 2) {
2190 i =
2191 #ifndef Sudden_Underflow
2192 denorm ? be + (Bias + (P-1) - 1 + 1) :
2193 #endif
2194 #ifdef IBM
2195 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2196 #else
2197 1 + P - bbits;
2198 #endif
2199 } else {
2200 j = ilim - 1;
2201 if (m5 >= j)
2202 m5 -= j;
2203 else {
2204 s5 += j -= m5;
2205 b5 += j;
2206 m5 = 0;
2207 }
2208 if ((i = ilim) < 0) {
2209 m2 -= i;
2210 i = 0;
2211 }
2212 }
2213 b2 += i;
2214 s2 += i;
2215 mhi = i2b(1);
2216 }
2217 if (m2 > 0 && s2 > 0) {
2218 i = m2 < s2 ? m2 : s2;
2219 b2 -= i;
2220 m2 -= i;
2221 s2 -= i;
2222 }
2223 if (b5 > 0) {
2224 if (leftright) {
2225 if (m5 > 0) {
2226 mhi = pow5mult(mhi, m5);
2227 b1 = mult(mhi, b);
2228 Bfree(b);
2229 b = b1;
2230 }
2231 if (j = b5 - m5)
2232 b = pow5mult(b, j);
2233 } else
2234 b = pow5mult(b, b5);
2235 }
2236 S = i2b(1);
2237 if (s5 > 0)
2238 S = pow5mult(S, s5);
2239
2240 /* Check for special case that d is a normalized power of 2. */
2241
2242 if (mode < 2) {
2243 if (!word1(d) && !(word0(d) & Bndry_mask)
2244 #ifndef Sudden_Underflow
2245 && word0(d) & Exp_mask
2246 #endif
2247 ) {
2248 /* The special case */
2249 b2 += Log2P;
2250 s2 += Log2P;
2251 spec_case = 1;
2252 } else
2253 spec_case = 0;
2254 }
2255
2256 /* Arrange for convenient computation of quotients:
2257 * shift left if necessary so divisor has 4 leading 0 bits.
2258 *
2259 * Perhaps we should just compute leading 28 bits of S once
2260 * and for all and pass them and a shift to quorem, so it
2261 * can do shifts and ors to compute the numerator for q.
2262 */
2263 #ifdef Pack_32
2264 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
2265 i = 32 - i;
2266 #else
2267 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2268 i = 16 - i;
2269 #endif
2270 if (i > 4) {
2271 i -= 4;
2272 b2 += i;
2273 m2 += i;
2274 s2 += i;
2275 } else if (i < 4) {
2276 i += 28;
2277 b2 += i;
2278 m2 += i;
2279 s2 += i;
2280 }
2281 if (b2 > 0)
2282 b = lshift(b, b2);
2283 if (s2 > 0)
2284 S = lshift(S, s2);
2285 if (k_check) {
2286 if (cmp(b,S) < 0) {
2287 k--;
2288 b = multadd(b, 10, 0); /* we botched the k estimate */
2289 if (leftright)
2290 mhi = multadd(mhi, 10, 0);
2291 ilim = ilim1;
2292 }
2293 }
2294 if (ilim <= 0 && mode > 2) {
2295 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2296 /* no digits, fcvt style */
2297 no_digits:
2298 k = -1 - ndigits;
2299 goto ret;
2300 }
2301 one_digit:
2302 *s++ = '1';
2303 k++;
2304 goto ret;
2305 }
2306 if (leftright) {
2307 if (m2 > 0)
2308 mhi = lshift(mhi, m2);
2309
2310 /* Compute mlo -- check for special case
2311 * that d is a normalized power of 2.
2312 */
2313
2314 mlo = mhi;
2315 if (spec_case) {
2316 mhi = Balloc(mhi->k);
2317 Bcopy(mhi, mlo);
2318 mhi = lshift(mhi, Log2P);
2319 }
2320
2321 for (i = 1;;i++) {
2322 dig = quorem(b,S) + '0';
2323 /* Do we yet have the shortest decimal string
2324 * that will round to d?
2325 */
2326 j = cmp(b, mlo);
2327 delta = diff(S, mhi);
2328 j1 = delta->sign ? 1 : cmp(b, delta);
2329 Bfree(delta);
2330 #ifndef ROUND_BIASED
2331 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2332 if (dig == '9')
2333 goto round_9_up;
2334 if (j > 0)
2335 dig++;
2336 *s++ = dig;
2337 goto ret;
2338 }
2339 #endif
2340 if (j < 0 || j == 0 && !mode
2341 #ifndef ROUND_BIASED
2342 && !(word1(d) & 1)
2343 #endif
2344 ) {
2345 if (j1 > 0) {
2346 b = lshift(b, 1);
2347 j1 = cmp(b, S);
2348 if ((j1 > 0 || j1 == 0 && dig & 1)
2349 && dig++ == '9')
2350 goto round_9_up;
2351 }
2352 *s++ = dig;
2353 goto ret;
2354 }
2355 if (j1 > 0) {
2356 if (dig == '9') { /* possible if i == 1 */
2357 round_9_up:
2358 *s++ = '9';
2359 goto roundoff;
2360 }
2361 *s++ = dig + 1;
2362 goto ret;
2363 }
2364 *s++ = dig;
2365 if (i == ilim)
2366 break;
2367 b = multadd(b, 10, 0);
2368 if (mlo == mhi)
2369 mlo = mhi = multadd(mhi, 10, 0);
2370 else {
2371 mlo = multadd(mlo, 10, 0);
2372 mhi = multadd(mhi, 10, 0);
2373 }
2374 }
2375 } else
2376 for (i = 1;; i++) {
2377 *s++ = dig = quorem(b,S) + '0';
2378 if (i >= ilim)
2379 break;
2380 b = multadd(b, 10, 0);
2381 }
2382
2383 /* Round off last digit */
2384
2385 b = lshift(b, 1);
2386 j = cmp(b, S);
2387 if (j > 0 || j == 0 && dig & 1) {
2388 roundoff:
2389 while (*--s == '9')
2390 if (s == s0) {
2391 k++;
2392 *s++ = '1';
2393 goto ret;
2394 }
2395 ++*s++;
2396 } else {
2397 while (*--s == '0');
2398 s++;
2399 }
2400 ret:
2401 Bfree(S);
2402 if (mhi) {
2403 if (mlo && mlo != mhi)
2404 Bfree(mlo);
2405 Bfree(mhi);
2406 }
2407 ret1:
2408 Bfree(b);
2409 if (s == s0) { /* don't return empty string */
2410 *s++ = '0';
2411 k = 0;
2412 }
2413 *s = 0;
2414 *decpt = k + 1;
2415 if (rve)
2416 *rve = s;
2417 return s0;
2418 }
2419 #ifdef __cplusplus
2420 }
2421 #endif
2422