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