9  * is included in all copies of any software which is or includes a copy
22 	AT&T Bell Laboratories, Room 2C-463
24 	Murray Hill, NJ 07974-2070
25 	U.S.A.
29 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
31  * This strtod returns a nearest machine number to the input decimal
33  * broken by the IEEE round-even rule.  Otherwise ties are broken by
37  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
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
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
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
60  * #define IEEE_8087 for IEEE-arithmetic machines where the least
62  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
64  * #define Sudden_Underflow for IEEE-format machines without gradual
66  * #define IBM for IBM mainframe-style floating-point arithmetic.
67  * #define VAX for VAX-style floating-point arithmetic.
69  * #define No_leftright to omit left-right logic in fast floating-point
73  *	that use extended-precision instructions to compute rounded
75  * #define ROUND_BIASED for IEEE-format with biased rounding.
76  * #define Inaccurate_Divide for IEEE-format with correctly rounded
78  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
84 #pragma GCC diagnostic ignored "-Wpragmas"
85 #pragma GCC diagnostic ignored "-Wunknown-warning-option"
86 #pragma GCC diagnostic ignored "-Wanalyzer-malloc-leak"
87 #pragma GCC diagnostic ignored "-Wanalyzer-use-of-uninitialized-value"
88 #pragma GCC diagnostic ignored "-Wanalyzer-out-of-bounds"
89 #pragma GCC diagnostic ignored "-Warray-bounds"
90 #pragma GCC diagnostic ignored "-Wanalyzer-null-dereference"
101    overrun attacks (CVE-2009-0689), same as other implementations of gdtoa
159 			  (x) * sizeof(rv->_x[0]));  in Balloc()
161   rv->_k = k;  in Balloc()
162   rv->_maxwds = x;  in Balloc()
163   rv->_sign = rv->_wds = 0;  in Balloc()
177 	int a)  in multadd()  argument
189   wds = b->_wds;  in multadd()
190   x = b->_x;  in multadd()
196       y = (xi & 0xffff) * m + a;  in multadd()
198       a = (int) (z >> 16);  in multadd()
201       y = *x * m + a;  in multadd()
202       a = (int) (y >> 16);  in multadd()
207   if (a)  in multadd()
209       if (wds >= b->_maxwds)  in multadd()
211 	  b1 = Balloc (b->_k + 1);  in multadd()
220       b->_x[wds++] = a;  in multadd()
221       b->_wds = wds;  in multadd()
243   b->_x[0] = y9;  in s2b()
244   b->_wds = 1;  in s2b()
249   b->_x[0] = y9 & 0xffff;  in s2b()
250   b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;  in s2b()
258 	b = multadd (b, 10, *s++ - '0');  in s2b()
265     b = multadd (b, 10, *s++ - '0');  in s2b()
360     b->_x[0] = i;  in i2b()
361     b->_wds = 1;  in i2b()
367 mult (_Bigint * a, _Bigint * b)  in mult()  argument
377   if (!a || !b)  in mult()
380   if (a->_wds < b->_wds)  in mult()
382       c = a;  in mult()
383       a = b;  in mult()
386   k = a->_k;  in mult()
387   wa = a->_wds;  in mult()
388   wb = b->_wds;  in mult()
390   if (wc > a->_maxwds)  in mult()
395   for (x = c->_x, xa = x + wc; x < xa; x++)  in mult()
397   xa = a->_x;  in mult()
399   xb = b->_x;  in mult()
401   xc0 = c->_x;  in mult()
458   for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);  in mult()
459   c->_wds = wc;  in mult()
471     b = multadd (b, p05[i - 1], 0);  in pow5mult()
509   k1 = b->_k;  in lshift()
510   n1 = n + b->_wds + 1;  in lshift()
511   for (i = b->_maxwds; n1 > i; i <<= 1)  in lshift()
517   x1 = b1->_x;  in lshift()
520   x = b->_x;  in lshift()
521   xe = x + b->_wds;  in lshift()
525       k1 = 32 - k;  in lshift()
539       k1 = 16 - k;  in lshift()
555   b1->_wds = n1 - 1;  in lshift()
562 cmp (_Bigint * a, _Bigint * b)  in cmp()  argument
567   if (!a || !b)  in cmp()
570   i = a->_wds;  in cmp()
571   j = b->_wds;  in cmp()
573   if (i > 1 && !a->_x[i - 1])  in cmp()
574     Bug ("cmp called with a->_x[a->_wds-1] == 0");  in cmp()
575   if (j > 1 && !b->_x[j - 1])  in cmp()
576     Bug ("cmp called with b->_x[b->_wds-1] == 0");  in cmp()
578   if (i -= j)  in cmp()
580   xa0 = a->_x;  in cmp()
582   xb0 = b->_x;  in cmp()
586       if (*--xa != *--xb)  in cmp()
587 	return *xa < *xb ? -1 : 1;  in cmp()
596 	_Bigint * a, _Bigint * b)  in diff()  argument
606   i = cmp (a, b);  in diff()
612       c->_wds = 1;  in diff()
613       c->_x[0] = 0;  in diff()
618       c = a;  in diff()
619       a = b;  in diff()
625   c = Balloc (a->_k);  in diff()
628   c->_sign = i;  in diff()
629   wa = a->_wds;  in diff()
630   xa = a->_x;  in diff()
632   wb = b->_wds;  in diff()
633   xb = b->_x;  in diff()
635   xc = c->_x;  in diff()
640       y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;  in diff()
643       z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;  in diff()
662       y = *xa++ - *xb++ + borrow;  in diff()
676   while (!*--xc)  in diff()
677     wa--;  in diff()
678   c->_wds = wa;  in diff()
685   union double_union x, a;  in ulp()  local
690   L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;  in ulp()
698       word0 (a) = L;  in ulp()
700       word1 (a) = 0;  in ulp()
707       L = -L >> Exp_shift;  in ulp()
710 	  word0 (a) = 0x80000 >> L;  in ulp()
712 	  word1 (a) = 0;  in ulp()
717 	  word0 (a) = 0;  in ulp()
718 	  L -= Exp_shift;  in ulp()
720           word1 (a) = L >= 31 ? 1 : (__uint32_t) 1 << (31 - L);  in ulp()
725   return a.d;  in ulp()
729 b2d (_Bigint * a, int *e)  in b2d()  argument
744   xa0 = a->_x;  in b2d()
745   xa = xa0 + a->_wds;  in b2d()
746   y = *--xa;  in b2d()
752   *e = 32 - k;  in b2d()
756       d0 = Exp_1 | y >> (Ebits - k);  in b2d()
758       w = xa > xa0 ? *--xa : 0;  in b2d()
759       d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);  in b2d()
763   z = xa > xa0 ? *--xa : 0;  in b2d()
764   if (k -= Ebits)  in b2d()
766       d0 = Exp_1 | y << k | z >> (32 - k);  in b2d()
767       y = xa > xa0 ? *--xa : 0;  in b2d()
769       d1 = z << k | y >> (32 - k);  in b2d()
782       z = xa > xa0 ? *--xa : 0;  in b2d()
783       d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;  in b2d()
784       w = xa > xa0 ? *--xa : 0;  in b2d()
785       y = xa > xa0 ? *--xa : 0;  in b2d()
786       d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;  in b2d()
789   z = xa > xa0 ? *--xa : 0;  in b2d()
790   w = xa > xa0 ? *--xa : 0;  in b2d()
791   k -= Ebits + 16;  in b2d()
792   d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;  in b2d()
793   y = xa > xa0 ? *--xa : 0;  in b2d()
841   x = b->_x;  in d2b()
862          x[0] = y | z << (32 - k);  in d2b()
867       i = b->_wds = (x[1] = z) ? 2 : 1;  in d2b()
878       i = b->_wds = 1;  in d2b()
891 	    x[0] = y | z << 32 - k & 0xffff;  in d2b()
892 	    x[1] = z >> k - 16 & 0xffff;  in d2b()
899 	    x[1] = y >> 16 | z << 16 - k & 0xffff;  in d2b()
934     --i;  in d2b()
935   b->_wds = i + 1;  in d2b()
942       *e = (de - Bias - (P - 1) << 2) + k;  in d2b()
943       *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);  in d2b()
945       *e = de - Bias - (P - 1) + k;  in d2b()
946       *bits = P - k;  in d2b()
952       *e = de - Bias - (P - 1) + 1 + k;  in d2b()
954       *bits = 32 * i - hi0bits (x[i - 1]);  in d2b()
956       *bits = (i + 2) * 16 - hi0bits (x[i]);  in d2b()
966 ratio (_Bigint * a, _Bigint * b)  in ratio()  argument
972   da.d = b2d (a, &ka);  in ratio()
975   k = ka - kb + 32 * (a->_wds - b->_wds);  in ratio()
977   k = ka - kb + 16 * (a->_wds - b->_wds);  in ratio()
988       k = -k;  in ratio()
998       k = -k;  in ratio()
1020 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
1026 {1e-16, 1e-32};
1039       dig--;  in _mprec_log10()
1054 	ce = c + ((n-1) >> kshift) + 1;  in copybits()
1055 	x = b->_x;  in copybits()
1057 	xe = x + b->_wds;  in copybits()
1061 	nw = b->_wds;  in copybits()
1063 	for(xe = x + (nw - nw1); x < xe; x += 2)  in copybits()
1079 	x = b->_x;  in any_on()
1080 	nwds = b->_wds;  in any_on()
1094 		if (*--x)  in any_on()