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