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_8087 for IEEE-arithmetic machines where the least
61 * significant byte has the lowest address.
62 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
63 * significant byte has the lowest address.
64 * #define Sudden_Underflow for IEEE-format machines without gradual
65 * underflow (i.e., that flush to zero on underflow).
66 * #define IBM for IBM mainframe-style floating-point arithmetic.
67 * #define VAX for VAX-style floating-point arithmetic.
68 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
69 * #define No_leftright to omit left-right logic in fast floating-point
70 * computation of dtoa.
71 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
72 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
73 * that use extended-precision instructions to compute rounded
74 * products and quotients) with IBM.
75 * #define ROUND_BIASED for IEEE-format with biased rounding.
76 * #define Inaccurate_Divide for IEEE-format with correctly rounded
77 * products but inaccurate quotients, e.g., for Intel i860.
78 * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
79 * integer arithmetic. Whether this speeds things up or slows things
80 * down depends on the machine and the number being converted.
81 */
82
83 #ifdef __GNUC__
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"
91 #endif
92
93 #define _DEFAULT_SOURCE
94 #include <stdlib.h>
95 #include <string.h>
96 #include "mprec.h"
97 #include "atexit.h"
98
99 /* This is defined in sys/reent.h as (sizeof (size_t) << 3) now, as in NetBSD.
100 The old value of 15 was wrong and made newlib vulnerable against buffer
101 overrun attacks (CVE-2009-0689), same as other implementations of gdtoa
102 based on BSD code.
103 #define _Kmax 15
104 */
105 /* This value is used in stdlib/misc.c. reent/reent.c has to know it
106 as well to make sure the freelist is correctly free'd. Therefore
107 we define it here, rather than in stdlib/misc.c, as before. */
108 #define _Kmax (sizeof (size_t) << 3)
109
110 static NEWLIB_THREAD_LOCAL char *__dtoa_result;
111 static NEWLIB_THREAD_LOCAL int __dtoa_result_len;
112
__alloc_dtoa_result(int len)113 char *__alloc_dtoa_result(int len)
114 {
115 if (len > __dtoa_result_len) {
116 if (__mprec_register_exit() != 0)
117 return NULL;
118 if (__dtoa_result)
119 free(__dtoa_result);
120 __dtoa_result = malloc(len);
121 if (__dtoa_result)
122 __dtoa_result_len = len;
123 }
124 return __dtoa_result;
125 }
126
127 static NEWLIB_THREAD_LOCAL int _mprec_exit_registered;
128
129 static void
__mprec_exit(void)130 __mprec_exit(void)
131 {
132 if (__dtoa_result) {
133 free(__dtoa_result);
134 __dtoa_result = NULL;
135 __dtoa_result_len = 0;
136 }
137 }
138
139 int
__mprec_register_exit(void)140 __mprec_register_exit(void)
141 {
142 if (!_mprec_exit_registered) {
143 _mprec_exit_registered = 1;
144 return atexit(__mprec_exit);
145 }
146 return 0;
147 }
148
149 _Bigint *
Balloc(int k)150 Balloc (int k)
151 {
152 int x;
153 _Bigint *rv ;
154
155 x = 1 << k;
156 /* Allocate an mprec Bigint */
157 rv = (_Bigint *) calloc(1,
158 sizeof (_Bigint) +
159 (x) * sizeof(rv->_x[0]));
160 if (rv == NULL) return NULL;
161 rv->_k = k;
162 rv->_maxwds = x;
163 rv->_sign = rv->_wds = 0;
164 return rv;
165 }
166
167 void
Bfree(_Bigint * v)168 Bfree (_Bigint * v)
169 {
170 free(v);
171 }
172
173 _Bigint *
multadd(_Bigint * b,int m,int a)174 multadd (
175 _Bigint * b,
176 int m,
177 int a)
178 {
179 int i, wds;
180 __ULong *x, y;
181 #ifdef Pack_32
182 __ULong xi, z;
183 #endif
184 _Bigint *b1;
185
186 if (!b)
187 return NULL;
188
189 wds = b->_wds;
190 x = b->_x;
191 i = 0;
192 do
193 {
194 #ifdef Pack_32
195 xi = *x;
196 y = (xi & 0xffff) * m + a;
197 z = (xi >> 16) * m + (y >> 16);
198 a = (int) (z >> 16);
199 *x++ = (z << 16) + (y & 0xffff);
200 #else
201 y = *x * m + a;
202 a = (int) (y >> 16);
203 *x++ = y & 0xffff;
204 #endif
205 }
206 while (++i < wds);
207 if (a)
208 {
209 if (wds >= b->_maxwds)
210 {
211 b1 = Balloc (b->_k + 1);
212 if (!b1) {
213 Bfree(b);
214 return NULL;
215 }
216 Bcopy (b1, b);
217 Bfree (b);
218 b = b1;
219 }
220 b->_x[wds++] = a;
221 b->_wds = wds;
222 }
223 return b;
224 }
225
226 _Bigint *
s2b(const char * s,int nd0,int nd,__ULong y9)227 s2b (
228 const char *s,
229 int nd0,
230 int nd,
231 __ULong y9)
232 {
233 _Bigint *b;
234 int i, k;
235 __Long x, y;
236
237 x = (nd + 8) / 9;
238 for (k = 0, y = 1; x > y; y <<= 1, k++);
239 #ifdef Pack_32
240 b = Balloc (k);
241 if (!b)
242 return NULL;
243 b->_x[0] = y9;
244 b->_wds = 1;
245 #else
246 b = Balloc (k + 1);
247 if (!b)
248 return NULL;
249 b->_x[0] = y9 & 0xffff;
250 b->_wds = (b->_x[1] = y9 >> 16) ? 2 : 1;
251 #endif
252
253 i = 9;
254 if (9 < nd0)
255 {
256 s += 9;
257 do
258 b = multadd (b, 10, *s++ - '0');
259 while (++i < nd0);
260 s++;
261 }
262 else
263 s += 10;
264 for (; i < nd; i++)
265 b = multadd (b, 10, *s++ - '0');
266 return b;
267 }
268
269 int
hi0bits(register __ULong x)270 hi0bits (register __ULong x)
271 {
272 register int k = 0;
273
274 if (!(x & 0xffff0000))
275 {
276 k = 16;
277 x <<= 16;
278 }
279 if (!(x & 0xff000000))
280 {
281 k += 8;
282 x <<= 8;
283 }
284 if (!(x & 0xf0000000))
285 {
286 k += 4;
287 x <<= 4;
288 }
289 if (!(x & 0xc0000000))
290 {
291 k += 2;
292 x <<= 2;
293 }
294 if (!(x & 0x80000000))
295 {
296 k++;
297 if (!(x & 0x40000000))
298 return 32;
299 }
300 return k;
301 }
302
303 int
lo0bits(__ULong * y)304 lo0bits (__ULong *y)
305 {
306 register int k;
307 register __ULong x = *y;
308
309 if (x & 7)
310 {
311 if (x & 1)
312 return 0;
313 if (x & 2)
314 {
315 *y = x >> 1;
316 return 1;
317 }
318 *y = x >> 2;
319 return 2;
320 }
321 k = 0;
322 if (!(x & 0xffff))
323 {
324 k = 16;
325 x >>= 16;
326 }
327 if (!(x & 0xff))
328 {
329 k += 8;
330 x >>= 8;
331 }
332 if (!(x & 0xf))
333 {
334 k += 4;
335 x >>= 4;
336 }
337 if (!(x & 0x3))
338 {
339 k += 2;
340 x >>= 2;
341 }
342 if (!(x & 1))
343 {
344 k++;
345 x >>= 1;
346 if (!(x & 1))
347 return 32;
348 }
349 *y = x;
350 return k;
351 }
352
353 _Bigint *
i2b(int i)354 i2b (int i)
355 {
356 _Bigint *b;
357
358 b = Balloc (1);
359 if (b) {
360 b->_x[0] = i;
361 b->_wds = 1;
362 }
363 return b;
364 }
365
366 _Bigint *
mult(_Bigint * a,_Bigint * b)367 mult (_Bigint * a, _Bigint * b)
368 {
369 _Bigint *c;
370 int k, wa, wb, wc;
371 __ULong carry, y, z;
372 __ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
373 #ifdef Pack_32
374 __ULong z2;
375 #endif
376
377 if (!a || !b)
378 return NULL;
379
380 if (a->_wds < b->_wds)
381 {
382 c = a;
383 a = b;
384 b = c;
385 }
386 k = a->_k;
387 wa = a->_wds;
388 wb = b->_wds;
389 wc = wa + wb;
390 if (wc > a->_maxwds)
391 k++;
392 c = Balloc (k);
393 if (!c)
394 return NULL;
395 for (x = c->_x, xa = x + wc; x < xa; x++)
396 *x = 0;
397 xa = a->_x;
398 xae = xa + wa;
399 xb = b->_x;
400 xbe = xb + wb;
401 xc0 = c->_x;
402 #ifdef Pack_32
403 for (; xb < xbe; xb++, xc0++)
404 {
405 if ((y = *xb & 0xffff) != 0)
406 {
407 x = xa;
408 xc = xc0;
409 carry = 0;
410 do
411 {
412 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
413 carry = z >> 16;
414 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
415 carry = z2 >> 16;
416 Storeinc (xc, z2, z);
417 }
418 while (x < xae);
419 *xc = carry;
420 }
421 if ((y = *xb >> 16) != 0)
422 {
423 x = xa;
424 xc = xc0;
425 carry = 0;
426 z2 = *xc;
427 do
428 {
429 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
430 carry = z >> 16;
431 Storeinc (xc, z, z2);
432 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
433 carry = z2 >> 16;
434 }
435 while (x < xae);
436 *xc = z2;
437 }
438 }
439 #else
440 for (; xb < xbe; xc0++)
441 {
442 if (y = *xb++)
443 {
444 x = xa;
445 xc = xc0;
446 carry = 0;
447 do
448 {
449 z = *x++ * y + *xc + carry;
450 carry = z >> 16;
451 *xc++ = z & 0xffff;
452 }
453 while (x < xae);
454 *xc = carry;
455 }
456 }
457 #endif
458 for (xc0 = c->_x, xc = xc0 + wc; wc > 0 && !*--xc; --wc);
459 c->_wds = wc;
460 return c;
461 }
462
463 _Bigint *
pow5mult(_Bigint * b,int k)464 pow5mult (_Bigint * b, int k)
465 {
466 _Bigint *b1, *p5, *p51;
467 int i;
468 static const int p05[3] = {5, 25, 125};
469
470 if ((i = k & 3) != 0)
471 b = multadd (b, p05[i - 1], 0);
472
473 if (!(k >>= 2))
474 return b;
475 p5 = i2b (625);
476 for (;;)
477 {
478 if (k & 1)
479 {
480 b1 = mult (b, p5);
481 Bfree (b);
482 b = b1;
483 }
484 if (!(k >>= 1))
485 break;
486 p51 = mult (p5, p5);
487 Bfree(p5);
488 p5 = p51;
489 }
490 Bfree(p5);
491 return b;
492 }
493
494 _Bigint *
lshift(_Bigint * b,int k)495 lshift (_Bigint * b, int k)
496 {
497 int i, k1, n, n1;
498 _Bigint *b1;
499 __ULong *x, *x1, *xe, z;
500
501 if (!b)
502 return NULL;
503
504 #ifdef Pack_32
505 n = k >> 5;
506 #else
507 n = k >> 4;
508 #endif
509 k1 = b->_k;
510 n1 = n + b->_wds + 1;
511 for (i = b->_maxwds; n1 > i; i <<= 1)
512 k1++;
513 b1 = Balloc (k1);
514 if (!b1)
515 goto bail;
516
517 x1 = b1->_x;
518 for (i = 0; i < n; i++)
519 *x1++ = 0;
520 x = b->_x;
521 xe = x + b->_wds;
522 #ifdef Pack_32
523 if (k &= 0x1f)
524 {
525 k1 = 32 - k;
526 z = 0;
527 do
528 {
529 *x1++ = *x << k | z;
530 z = *x++ >> k1;
531 }
532 while (x < xe);
533 if ((*x1 = z) != 0)
534 ++n1;
535 }
536 #else
537 if (k &= 0xf)
538 {
539 k1 = 16 - k;
540 z = 0;
541 do
542 {
543 *x1++ = *x << k & 0xffff | z;
544 z = *x++ >> k1;
545 }
546 while (x < xe);
547 if (*x1 = z)
548 ++n1;
549 }
550 #endif
551 else
552 do
553 *x1++ = *x++;
554 while (x < xe);
555 b1->_wds = n1 - 1;
556 bail:
557 Bfree (b);
558 return b1;
559 }
560
561 int
cmp(_Bigint * a,_Bigint * b)562 cmp (_Bigint * a, _Bigint * b)
563 {
564 __ULong *xa, *xa0, *xb, *xb0;
565 int i, j;
566
567 if (!a || !b)
568 return 0;
569
570 i = a->_wds;
571 j = b->_wds;
572 #ifdef DEBUG
573 if (i > 1 && !a->_x[i - 1])
574 Bug ("cmp called with a->_x[a->_wds-1] == 0");
575 if (j > 1 && !b->_x[j - 1])
576 Bug ("cmp called with b->_x[b->_wds-1] == 0");
577 #endif
578 if (i -= j)
579 return i;
580 xa0 = a->_x;
581 xa = xa0 + j;
582 xb0 = b->_x;
583 xb = xb0 + j;
584 for (;;)
585 {
586 if (*--xa != *--xb)
587 return *xa < *xb ? -1 : 1;
588 if (xa <= xa0)
589 break;
590 }
591 return 0;
592 }
593
594 _Bigint *
diff(_Bigint * a,_Bigint * b)595 diff (
596 _Bigint * a, _Bigint * b)
597 {
598 _Bigint *c;
599 int i, wa, wb;
600 __Long borrow, y; /* We need signed shifts here. */
601 __ULong *xa, *xae, *xb, *xbe, *xc;
602 #ifdef Pack_32
603 __Long z;
604 #endif
605
606 i = cmp (a, b);
607 if (!i)
608 {
609 c = Balloc (0);
610 if (!c)
611 return NULL;
612 c->_wds = 1;
613 c->_x[0] = 0;
614 return c;
615 }
616 if (i < 0)
617 {
618 c = a;
619 a = b;
620 b = c;
621 i = 1;
622 }
623 else
624 i = 0;
625 c = Balloc (a->_k);
626 if (!c)
627 return NULL;
628 c->_sign = i;
629 wa = a->_wds;
630 xa = a->_x;
631 xae = xa + wa;
632 wb = b->_wds;
633 xb = b->_x;
634 xbe = xb + wb;
635 xc = c->_x;
636 borrow = 0;
637 #ifdef Pack_32
638 do
639 {
640 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
641 borrow = y >> 16;
642 Sign_Extend (borrow, y);
643 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
644 borrow = z >> 16;
645 Sign_Extend (borrow, z);
646 Storeinc (xc, z, y);
647 }
648 while (xb < xbe);
649 while (xa < xae)
650 {
651 y = (*xa & 0xffff) + borrow;
652 borrow = y >> 16;
653 Sign_Extend (borrow, y);
654 z = (*xa++ >> 16) + borrow;
655 borrow = z >> 16;
656 Sign_Extend (borrow, z);
657 Storeinc (xc, z, y);
658 }
659 #else
660 do
661 {
662 y = *xa++ - *xb++ + borrow;
663 borrow = y >> 16;
664 Sign_Extend (borrow, y);
665 *xc++ = y & 0xffff;
666 }
667 while (xb < xbe);
668 while (xa < xae)
669 {
670 y = *xa++ + borrow;
671 borrow = y >> 16;
672 Sign_Extend (borrow, y);
673 *xc++ = y & 0xffff;
674 }
675 #endif
676 while (!*--xc)
677 wa--;
678 c->_wds = wa;
679 return c;
680 }
681
682 double
ulp(double _x)683 ulp (double _x)
684 {
685 union double_union x, a;
686 register __Long L;
687
688 x.d = _x;
689
690 L = (word0 (x) & Exp_mask) - (P - 1) * Exp_msk1;
691 #ifndef Sudden_Underflow
692 if (L > 0)
693 {
694 #endif
695 #ifdef IBM
696 L |= Exp_msk1 >> 4;
697 #endif
698 word0 (a) = L;
699 #ifndef _DOUBLE_IS_32BITS
700 word1 (a) = 0;
701 #endif
702
703 #ifndef Sudden_Underflow
704 }
705 else
706 {
707 L = -L >> Exp_shift;
708 if (L < Exp_shift)
709 {
710 word0 (a) = 0x80000 >> L;
711 #ifndef _DOUBLE_IS_32BITS
712 word1 (a) = 0;
713 #endif
714 }
715 else
716 {
717 word0 (a) = 0;
718 L -= Exp_shift;
719 #ifndef _DOUBLE_IS_32BITS
720 word1 (a) = L >= 31 ? 1 : (__uint32_t) 1 << (31 - L);
721 #endif
722 }
723 }
724 #endif
725 return a.d;
726 }
727
728 double
b2d(_Bigint * a,int * e)729 b2d (_Bigint * a, int *e)
730 {
731 __ULong *xa, *xa0, y, z;
732 #if !defined(Pack_32) || !defined(_DOUBLE_IS_32BITS)
733 __ULong w;
734 #endif
735 int k;
736 union double_union d;
737 #ifdef VAX
738 __ULong d0, d1;
739 #else
740 #define d0 word0(d)
741 #define d1 word1(d)
742 #endif
743
744 xa0 = a->_x;
745 xa = xa0 + a->_wds;
746 y = *--xa;
747 #ifdef DEBUG
748 if (!y)
749 Bug ("zero y in b2d");
750 #endif
751 k = hi0bits (y);
752 *e = 32 - k;
753 #ifdef Pack_32
754 if (k < Ebits)
755 {
756 d0 = Exp_1 | y >> (Ebits - k);
757 #ifndef _DOUBLE_IS_32BITS
758 w = xa > xa0 ? *--xa : 0;
759 d1 = y << ((32 - Ebits) + k) | w >> (Ebits - k);
760 #endif
761 goto ret_d;
762 }
763 z = xa > xa0 ? *--xa : 0;
764 if (k -= Ebits)
765 {
766 d0 = Exp_1 | y << k | z >> (32 - k);
767 y = xa > xa0 ? *--xa : 0;
768 #ifndef _DOUBLE_IS_32BITS
769 d1 = z << k | y >> (32 - k);
770 #endif
771 }
772 else
773 {
774 d0 = Exp_1 | y;
775 #ifndef _DOUBLE_IS_32BITS
776 d1 = z;
777 #endif
778 }
779 #else
780 if (k < Ebits + 16)
781 {
782 z = xa > xa0 ? *--xa : 0;
783 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
784 w = xa > xa0 ? *--xa : 0;
785 y = xa > xa0 ? *--xa : 0;
786 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
787 goto ret_d;
788 }
789 z = xa > xa0 ? *--xa : 0;
790 w = xa > xa0 ? *--xa : 0;
791 k -= Ebits + 16;
792 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
793 y = xa > xa0 ? *--xa : 0;
794 d1 = w << k + 16 | y << k;
795 #endif
796 ret_d:
797 #ifdef VAX
798 word0 (d) = d0 >> 16 | d0 << 16;
799 word1 (d) = d1 >> 16 | d1 << 16;
800 #else
801 #undef d0
802 #undef d1
803 #endif
804 return d.d;
805 }
806
807 _Bigint *
d2b(double _d,int * e,int * bits)808 d2b (
809 double _d,
810 int *e,
811 int *bits)
812
813 {
814 union double_union d;
815 _Bigint *b;
816 int de, i, k;
817 __ULong *x, z;
818 #if !defined(Pack_32) || !defined(_DOUBLE_IS_32BITS)
819 __ULong y;
820 #endif
821 #ifdef VAX
822 __ULong d0, d1;
823 #endif
824 d.d = _d;
825 #ifdef VAX
826 d0 = word0 (d) >> 16 | word0 (d) << 16;
827 d1 = word1 (d) >> 16 | word1 (d) << 16;
828 #else
829 #define d0 word0(d)
830 #define d1 word1(d)
831 d.d = _d;
832 #endif
833
834 #ifdef Pack_32
835 b = Balloc (1);
836 #else
837 b = Balloc (2);
838 #endif
839 if (!b)
840 return NULL;
841 x = b->_x;
842
843 z = d0 & Frac_mask;
844 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
845 #ifdef Sudden_Underflow
846 de = (int) (d0 >> Exp_shift);
847 #ifndef IBM
848 z |= Exp_msk11;
849 #endif
850 #else
851 if ((de = (int) (d0 >> Exp_shift)) != 0)
852 z |= Exp_msk1;
853 #endif
854 #ifdef Pack_32
855 #ifndef _DOUBLE_IS_32BITS
856 if (d1)
857 {
858 y = d1;
859 k = lo0bits (&y);
860 if (k)
861 {
862 x[0] = y | z << (32 - k);
863 z >>= k;
864 }
865 else
866 x[0] = y;
867 i = b->_wds = (x[1] = z) ? 2 : 1;
868 }
869 else
870 #endif
871 {
872 #ifdef DEBUG
873 if (!z)
874 Bug ("Zero passed to d2b");
875 #endif
876 k = lo0bits (&z);
877 x[0] = z;
878 i = b->_wds = 1;
879 #ifndef _DOUBLE_IS_32BITS
880 k += 32;
881 #endif
882 }
883 #else
884 if (d1)
885 {
886 y = d1;
887 k = lo0bits (&y);
888 if (k)
889 if (k >= 16)
890 {
891 x[0] = y | z << 32 - k & 0xffff;
892 x[1] = z >> k - 16 & 0xffff;
893 x[2] = z >> k;
894 i = 2;
895 }
896 else
897 {
898 x[0] = y & 0xffff;
899 x[1] = y >> 16 | z << 16 - k & 0xffff;
900 x[2] = z >> k & 0xffff;
901 x[3] = z >> k + 16;
902 i = 3;
903 }
904 else
905 {
906 x[0] = y & 0xffff;
907 x[1] = y >> 16;
908 x[2] = z & 0xffff;
909 x[3] = z >> 16;
910 i = 3;
911 }
912 }
913 else
914 {
915 #ifdef DEBUG
916 if (!z)
917 Bug ("Zero passed to d2b");
918 #endif
919 k = lo0bits (&z);
920 if (k >= 16)
921 {
922 x[0] = z;
923 i = 0;
924 }
925 else
926 {
927 x[0] = z & 0xffff;
928 x[1] = z >> 16;
929 i = 1;
930 }
931 k += 32;
932 }
933 while (!x[i])
934 --i;
935 b->_wds = i + 1;
936 #endif
937 #ifndef Sudden_Underflow
938 if (de)
939 {
940 #endif
941 #ifdef IBM
942 *e = (de - Bias - (P - 1) << 2) + k;
943 *bits = 4 * P + 8 - k - hi0bits (word0 (d) & Frac_mask);
944 #else
945 *e = de - Bias - (P - 1) + k;
946 *bits = P - k;
947 #endif
948 #ifndef Sudden_Underflow
949 }
950 else
951 {
952 *e = de - Bias - (P - 1) + 1 + k;
953 #ifdef Pack_32
954 *bits = 32 * i - hi0bits (x[i - 1]);
955 #else
956 *bits = (i + 2) * 16 - hi0bits (x[i]);
957 #endif
958 }
959 #endif
960 return b;
961 }
962 #undef d0
963 #undef d1
964
965 double
ratio(_Bigint * a,_Bigint * b)966 ratio (_Bigint * a, _Bigint * b)
967
968 {
969 union double_union da, db;
970 int k, ka, kb;
971
972 da.d = b2d (a, &ka);
973 db.d = b2d (b, &kb);
974 #ifdef Pack_32
975 k = ka - kb + 32 * (a->_wds - b->_wds);
976 #else
977 k = ka - kb + 16 * (a->_wds - b->_wds);
978 #endif
979 #ifdef IBM
980 if (k > 0)
981 {
982 word0 (da) += (k >> 2) * Exp_msk1;
983 if (k &= 3)
984 da.d *= 1 << k;
985 }
986 else
987 {
988 k = -k;
989 word0 (db) += (k >> 2) * Exp_msk1;
990 if (k &= 3)
991 db.d *= 1 << k;
992 }
993 #else
994 if (k > 0)
995 word0 (da) += k * Exp_msk1;
996 else
997 {
998 k = -k;
999 word0 (db) += k * Exp_msk1;
1000 }
1001 #endif
1002 return da.d / db.d;
1003 }
1004
1005
1006 const double
1007 tens[] =
1008 {
1009 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1010 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1011 1e20, 1e21, 1e22, 1e23, 1e24
1012
1013 };
1014
1015 #if !defined(_DOUBLE_IS_32BITS) && !defined(__v800)
1016 const double bigtens[] =
1017 {1e16, 1e32, 1e64, 1e128, 1e256};
1018
1019 const double tinytens[] =
1020 {1e-16, 1e-32, 1e-64, 1e-128, 1e-256};
1021 #else
1022 const double bigtens[] =
1023 {1e16, 1e32};
1024
1025 const double tinytens[] =
1026 {1e-16, 1e-32};
1027 #endif
1028
1029
1030 double
_mprec_log10(int dig)1031 _mprec_log10 (int dig)
1032 {
1033 double v = 1.0;
1034 if (dig < 24)
1035 return tens[dig];
1036 while (dig > 0)
1037 {
1038 v *= 10;
1039 dig--;
1040 }
1041 return v;
1042 }
1043
1044 void
copybits(__ULong * c,int n,_Bigint * b)1045 copybits (__ULong *c,
1046 int n,
1047 _Bigint *b)
1048 {
1049 __ULong *ce, *x, *xe;
1050 #ifdef Pack_16
1051 int nw, nw1;
1052 #endif
1053
1054 ce = c + ((n-1) >> kshift) + 1;
1055 x = b->_x;
1056 #ifdef Pack_32
1057 xe = x + b->_wds;
1058 while(x < xe)
1059 *c++ = *x++;
1060 #else
1061 nw = b->_wds;
1062 nw1 = nw & 1;
1063 for(xe = x + (nw - nw1); x < xe; x += 2)
1064 Storeinc(c, x[1], x[0]);
1065 if (nw1)
1066 *c++ = *x;
1067 #endif
1068 while(c < ce)
1069 *c++ = 0;
1070 }
1071
1072 __ULong
any_on(_Bigint * b,int k)1073 any_on (_Bigint *b,
1074 int k)
1075 {
1076 int n, nwds;
1077 __ULong *x, *x0, x1, x2;
1078
1079 x = b->_x;
1080 nwds = b->_wds;
1081 n = k >> kshift;
1082 if (n > nwds)
1083 n = nwds;
1084 else if (n < nwds && (k &= kmask)) {
1085 x1 = x2 = x[n];
1086 x1 >>= k;
1087 x1 <<= k;
1088 if (x1 != x2)
1089 return 1;
1090 }
1091 x0 = x;
1092 x += n;
1093 while(x > x0)
1094 if (*--x)
1095 return 1;
1096 return 0;
1097 }
1098
1099
