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