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