1 /*
2 * Copyright © 1988-2004 Keith Packard and Bart Massey.
3 * All Rights Reserved. See the file COPYING in this directory
4 * for licensing information.
5 */
6
7 /*
8 * natural.c
9 *
10 * arithmetic for natural numbers
11 */
12
13 #include <math.h>
14 #include <stdio.h>
15 #include "nickle.h"
16
17 # define length(n) ((n)->length)
18 # define data(n) NaturalDigits(n)
19 # define max(a,b) ((a) > (b) ? (a) : (b))
20 # define zerop(n) (length(n) == 0)
21
22 #define BIGDOUBLE ((double)1.79769313486231470e+308)
23
24 #ifndef MAXDOUBLE
25 # define MAXDOUBLE ((double)1.79769313486231470e+308)
26 #endif
27
28 Natural *zero_natural;
29 Natural *one_natural;
30 Natural *two_natural;
31 Natural *max_int_natural;
32 Natural *max_signed_digit_natural;
33 #ifndef LBASE10
34 static Natural *max_tenpow_natural;
35 static int tenpow_digits;
36 #endif
37 DataCachePtr naturalCache;
38
39 int
NaturalToInt(Natural * n)40 NaturalToInt (Natural *n)
41 {
42 int i;
43 digit *d;
44 int index;
45
46 d = data(n) + length (n);
47 i = 0;
48 for (index = 0; index < length(n); index++)
49 i = i * BASE + (int) *--d;
50 return i;
51 }
52
53 double_digit
NaturalToDoubleDigit(Natural * n)54 NaturalToDoubleDigit(Natural *n)
55 {
56 double_digit i;
57 digit *d;
58 int index;
59
60 d = data(n) + length (n);
61 i = 0;
62 for (index = 0; index < length(n); index++)
63 i = i * BASE + (double_digit) *--d;
64 return i;
65 }
66
67 void
NaturalCopy(Natural * a,Natural * b)68 NaturalCopy (Natural *a, Natural *b)
69 {
70 digit *ad, *bd;
71 int index;
72
73 length(b) = length(a);
74 ad = data(a);
75 bd = data(b);
76 for (index = 0; index < length(a); index++)
77 *bd++ = *ad++;
78 }
79
80 #if 0
81 static void
82 NaturalClear (Natural *n)
83 {
84 int i;
85
86 for (i = 0; i < length(n); i++)
87 data(n)[i] = 0;
88 }
89 #endif
90
91 Bool
NaturalEven(Natural * n)92 NaturalEven (Natural *n)
93 {
94 if (!length (n) || (data(n)[0] & 1) == 0)
95 return True;
96 return False;
97 }
98
99 Bool
NaturalZero(Natural * n)100 NaturalZero (Natural *n)
101 {
102 return length (n) == 0;
103 }
104
105 #if 0
106 static int
107 NaturalOne (Natural *n)
108 {
109 return length (n) == 1 && data(n)[0] == 1;
110 }
111 #endif
112
113 Bool
NaturalLess(Natural * a,Natural * b)114 NaturalLess (Natural *a, Natural *b)
115 {
116 int index;
117 digit *at, *bt;
118
119 if (length (a) < length (b))
120 return True;
121 else if (length (b) < length (a))
122 return False;
123 else {
124 at = data(a) + length(a) - 1;
125 bt = data(b) + length(b) - 1;
126 for (index = 0; index < length(a); index++) {
127 if (*at < *bt)
128 return True;
129 else if (*bt < *at)
130 return False;
131 at--; bt--;
132 }
133 return False;
134 }
135 }
136
137 Bool
NaturalEqual(Natural * a,Natural * b)138 NaturalEqual (Natural *a, Natural *b)
139 {
140 int index;
141 digit *at, *bt;
142
143 if (length (a) == length (b)) {
144 at = data(a);
145 bt = data(b);
146 for (index = 0; index < length(a); index++)
147 if (*at++ != *bt++)
148 return False;
149 return True;
150 }
151 return False;
152 }
153
154 /*
155 * Primitive functions that operate on sequences
156 * of digits
157 */
158
159 static int
DigitsLen(digit * x,int len)160 DigitsLen (digit *x, int len)
161 {
162 x += (len - 1);
163 while (len && *x == 0)
164 {
165 len--;
166 x--;
167 }
168 return len;
169 }
170
171 static int
DigitsAdd(digit * x,int xlen,digit * y,int ylen,digit * r_orig)172 DigitsAdd (digit *x, int xlen, digit *y, int ylen, digit *r_orig)
173 {
174 digit *r = r_orig;
175 int rlen;
176 digit carry = 0;
177 digit xv, yv, rv;
178
179 while (xlen && ylen)
180 {
181 xlen--;
182 ylen--;
183 rv = xv = *x++;
184 yv = *y++ + carry;
185 if (yv)
186 {
187 carry = 0;
188 if ((rv = xv + yv) < xv)
189 carry = 1;
190 }
191 *r++ = rv;
192 }
193 while (ylen)
194 {
195 ylen--;
196 yv = *y++ + carry;
197 if (yv)
198 carry = 0;
199 *r++ = yv;
200 }
201 while (xlen)
202 {
203 xlen--;
204 rv = xv = *x++;
205 if (carry)
206 {
207 yv = carry;
208 carry = 0;
209 if ((rv = xv + yv) < xv)
210 carry = 1;
211 }
212 *r++ = rv;
213 }
214 if (carry)
215 *r++ = carry;
216 rlen = r - r_orig;
217 r--;
218 while (rlen && *r == 0)
219 {
220 r--;
221 rlen--;
222 }
223 return rlen;
224 }
225
226 static int
DigitsAddInPlace(digit * x_orig,int xlen,digit * y,int ylen,int off)227 DigitsAddInPlace (digit *x_orig, int xlen, digit *y, int ylen, int off)
228 {
229 digit *x = x_orig;
230 digit carry = 0;
231 digit xv, yv;
232
233 x += off;
234 xlen -= off;
235 if (xlen < 0)
236 {
237 x += xlen;
238 while (xlen++)
239 *x++ = 0;
240 }
241 while (xlen && ylen)
242 {
243 xlen--;
244 ylen--;
245 yv = *y++ + carry;
246 if (yv)
247 {
248 carry = 0;
249 xv = *x;
250 if ((*x = xv + yv) < xv)
251 carry = 1;
252 }
253 x++;
254 }
255 while (ylen)
256 {
257 ylen--;
258 yv = *y++ + carry;
259 if (yv)
260 carry = 0;
261 *x++ = yv;
262 }
263 while (xlen && carry)
264 {
265 xlen--;
266 xv = *x;
267 yv = carry;
268 carry = 0;
269 if ((*x = xv + yv) < xv)
270 carry = 1;
271 x++;
272 }
273 if (carry)
274 *x++ = carry;
275 return xlen + (x - x_orig);
276 }
277
278 static int
DigitsSubInPlace(digit * x_orig,int xlen,digit * y,int ylen,int off)279 DigitsSubInPlace (digit *x_orig, int xlen, digit *y, int ylen, int off)
280 {
281 digit *x = x_orig;
282 digit carry = 0;
283 digit xv, yv;
284
285 x += off;
286 xlen -= off;
287 while (ylen--)
288 {
289 xlen--;
290 xv = *x;
291 yv = *y++ + carry;
292 if (yv)
293 {
294 carry = 0;
295 if ((*x = xv - yv) > xv)
296 carry = 1;
297 }
298 x++;
299 }
300 while (carry)
301 {
302 xlen--;
303 xv = *x;
304 yv = carry;
305 carry = 0;
306 if ((*x = xv - yv) > xv)
307 carry = 1;
308 x++;
309 }
310 return xlen + (x - x_orig);
311 }
312
313 static int
DigitTimes(digit * x,int xlen,digit y,digit * result)314 DigitTimes (digit *x, int xlen, digit y, digit *result)
315 {
316 double_digit q;
317 digit carry;
318 int rlen = xlen;
319
320 if (y == 1)
321 {
322 memcpy (x, result, xlen * sizeof (digit));
323 return xlen;
324 }
325 carry = 0;
326 while (xlen--)
327 {
328 q = (double_digit) y * (double_digit) *x++ + (double_digit) carry;
329 carry = DivBase (q);
330 *result++ = ModBase (q);
331 }
332 if (carry)
333 {
334 *result++ = carry;
335 rlen++;
336 }
337 return rlen;
338 }
339
340 static int
DigitsGradeSchool(digit * x_orig,int xlen,digit * y_orig,int ylen,digit * result)341 DigitsGradeSchool (digit *x_orig, int xlen, digit *y_orig, int ylen, digit *result)
342 {
343 digit *x, *y, *r, *rbase, *rloop;
344 double_digit temp;
345 digit carry;
346 digit xd;
347 int xindex, yindex;
348 int rlen;
349
350 if (xlen == 0 || ylen == 0)
351 return 0;
352 if (xlen == 1)
353 return DigitTimes (y_orig, ylen, *x_orig, result);
354 if (ylen == 1)
355 return DigitTimes (x_orig, xlen, *y_orig, result);
356 memset (result, 0, (xlen + ylen + 1) * sizeof (digit));
357 rbase = 0;
358 x = x_orig;
359 xindex = xlen;
360 rbase = result;
361 while (xindex--)
362 {
363 carry = 0;
364 rloop = rbase++;
365 xd = *x++;
366
367 y = y_orig;
368 yindex = ylen;
369 while (yindex--)
370 {
371 temp = (double_digit) xd * (double_digit) *y++ + (double_digit) carry;
372 carry = DivBase (temp);
373 temp = ModBase (temp);
374 r = rloop++;
375 while (temp)
376 {
377 temp += (double_digit) *r;
378 *r++ = ModBase (temp);
379 temp = DivBase (temp);
380 }
381 }
382 if (carry)
383 {
384 r = rloop;
385 temp = carry;
386 while (temp)
387 {
388 temp += (double_digit) *r;
389 *r++ = ModBase (temp);
390 temp = DivBase (temp);
391 }
392 }
393 }
394 rlen = xlen + ylen + 1;
395 r = result + (rlen - 1);
396 while (rlen && *r == 0)
397 {
398 rlen--;
399 r--;
400 }
401 return rlen;
402 }
403
404 #define KARATSUBA_LIMIT 100
405
406 /*
407 * Karatsuba multiplication as found in
408 @article{ karatsuba62multiplication,
409 author = "A. Karatsuba and Yu Ofman",
410 title = "Multiplication of multidigit numbers on automata",
411 journal = "Doklady Akademii Nauk SSSR",
412 volume = "145",
413 number = "2",
414 pages = "293--294",
415 year = "1962"
416 }
417 */
418 static int
DigitsKaratsuba(digit * x,int xlen,digit * y,int ylen,digit * result,digit * tmp)419 DigitsKaratsuba (digit *x, int xlen, digit *y, int ylen, digit *result, digit *tmp)
420 {
421 /*
422 * x * y = (x1 * b + x0) * (y1 * b + y0);
423 * = b^2 x1 y1 + b (x1 y0 + x0 y1) + x0 y0
424 * = b^2 x1 y1 + b (x1 y0 + x0 y1 + x1 y1 + x0 y0) + x0 y0 - b x1 y1 - b x0 y0
425 * = (b^2 - b) x1 y1 + b (x1 + x0) (y0 + y1) + (1 - b) x0 y0
426 */
427 int off;
428 int off2;
429 digit *x1, *x0, *y1, *y0;
430 digit *f, *m1, *m2;
431 digit *next_tmp;
432 int x1len, x0len, y1len, y0len;
433 int flen, m1len, m2len;
434 int rlen;
435
436 if (aborting)
437 return 0;
438
439 if (xlen < KARATSUBA_LIMIT || ylen < KARATSUBA_LIMIT)
440 return DigitsGradeSchool (x, xlen, y, ylen, result);
441
442 off = xlen > ylen ? (xlen >> 1) : (ylen >> 1);
443 off2 = off << 1;
444 /*
445 * Normalize partial quotients
446 */
447 x0 = x;
448 x0len = xlen;
449 if (x0len > off)
450 x0len = DigitsLen (x0, off);
451 if (off < xlen)
452 {
453 x1 = x + off;
454 x1len = DigitsLen (x1, xlen - off);
455 }
456 else
457 {
458 x1 = x0;
459 x1len = 0;
460 }
461
462 y0 = y;
463 y0len = ylen;
464 if (y0len > off)
465 y0len = DigitsLen (y0, off);
466 if (off < ylen)
467 {
468 y1 = y + off;
469 y1len = DigitsLen (y1, ylen - off);
470 }
471 else
472 {
473 y1 = y0;
474 y1len = 0;
475 }
476
477 /*
478 * Allocate temp space
479 */
480 m1 = tmp;
481 m2 = m1 + off + 1;
482 f = tmp; /* overlay first factor on minuends */
483 next_tmp = m2 + off + 1;
484
485 /*
486 * Generate middle factor first
487 */
488 m1len = DigitsAdd (x0, x0len, x1, x1len, m1);
489 m2len = DigitsAdd (y0, y0len, y1, y1len, m2);
490
491 /*
492 * Compute middle factor
493 */
494 rlen = 0;
495 if (m1len && m2len)
496 {
497 memset (result, 0, off * sizeof (digit));
498 rlen = DigitsKaratsuba (m1, m1len, m2, m2len, result + off, next_tmp) + off;
499 if (aborting)
500 return rlen;
501 }
502
503 /*
504 * Compute first factor
505 */
506 if (x1len && y1len)
507 {
508 flen = DigitsKaratsuba (x1, x1len, y1, y1len, f, next_tmp);
509 if (aborting)
510 return rlen;
511 rlen = DigitsAddInPlace (result, rlen, f, flen, off2);
512 rlen = DigitsSubInPlace (result, rlen, f, flen, off);
513 }
514
515 /*
516 * Compute third factor
517 */
518
519 if (x0len && y0len)
520 {
521 flen = DigitsKaratsuba (x0, x0len, y0, y0len, f, next_tmp);
522 if (aborting)
523 return rlen;
524 rlen = DigitsAddInPlace (result, rlen, f, flen, 0);
525 rlen = DigitsSubInPlace (result, rlen, f, flen, off);
526 }
527 return rlen;
528 }
529
530 Natural *
NaturalPlus(Natural * a,Natural * b)531 NaturalPlus (Natural *a, Natural *b)
532 {
533 ENTER ();
534 Natural *result;
535
536 result = AllocNatural (max(length(a), length(b)) + 1);
537 result->length = DigitsAdd (data(a), length(a),
538 data(b), length(b),
539 data(result));
540 RETURN (result);
541 }
542
543 Natural *
NaturalMinus(Natural * a,Natural * b)544 NaturalMinus (Natural *a, Natural *b)
545 {
546 ENTER ();
547 int resultlen;
548 Natural *result;
549 signed_digit temp, carry;
550 digit *at, *bt, *rt;
551 int index;
552 int len;
553
554 resultlen = length(a);
555 result = AllocNatural (resultlen);
556 at = data(a);
557 bt = data(b);
558 rt = data(result);
559 carry = 0;
560 len = -1;
561 for (index = 0; index < resultlen; index++) {
562 temp = ((signed_digit) (index < length(a) ? *at++ : 0) -
563 (signed_digit) (index < length(b) ? *bt++ : 0) -
564 (signed_digit) carry);
565 carry = 0;
566 if (temp < 0) {
567 temp += BASE;
568 carry = 1;
569 }
570 if (temp > 0)
571 len = index;
572 *rt++ = temp;
573 }
574 length(result) = len + 1;
575 RETURN(result);
576 }
577
578 Natural *
NaturalTimes(Natural * a,Natural * b)579 NaturalTimes (Natural *a, Natural *b)
580 {
581 ENTER ();
582 Natural *result;
583 int rlen;
584 digit *tmp;
585 int tmp_len;
586
587 if (length (a) < KARATSUBA_LIMIT || length (b) < KARATSUBA_LIMIT)
588 {
589 if (zeroNp (a) || zeroNp (b))
590 RETURN (zero_natural);
591 if (oneNp (a))
592 RETURN(b);
593 if (oneNp (b))
594 RETURN (a);
595 result = AllocNatural (length(a) + length (b) + 1);
596 result->length = DigitsGradeSchool (data(a), length(a), data(b), length (b), data(result));
597 }
598 else
599 {
600 if (length (a) > length (b))
601 rlen = length (a) << 1;
602 else
603 rlen = length (b) << 1;
604 result = AllocNatural (rlen);
605 tmp_len = rlen << 3;
606 tmp = AllocateTemp (tmp_len * sizeof (digit));
607 rlen = DigitsKaratsuba (data(a), length (a), data(b), length (b), data(result), tmp);
608 if (aborting)
609 RETURN(zero_natural);
610 tmp = data(result) + (rlen - 1);
611 while (rlen && *tmp == 0)
612 {
613 rlen--;
614 tmp--;
615 }
616 result->length = rlen;
617 }
618
619 RETURN (result);
620 }
621
622 Natural *
NaturalLand(Natural * a,Natural * b)623 NaturalLand (Natural *a, Natural *b)
624 {
625 ENTER ();
626 digit *at, *bt, *rt;
627 Natural *result;
628 int resultlen;
629
630 resultlen = length (a);
631 if (resultlen > length(b))
632 resultlen = length(b);
633 at = data(a) + (resultlen-1);
634 bt = data(b) + (resultlen-1);
635 while (resultlen > 0 && (*at & *bt) == 0)
636 {
637 resultlen--;
638 at--;
639 bt--;
640 }
641 if (resultlen == 0)
642 RETURN (zero_natural);
643 result = AllocNatural (resultlen);
644 rt = data(result) + (resultlen-1);
645 while (resultlen-- > 0)
646 *rt-- = *at-- & *bt--;
647 RETURN (result);
648 }
649
650 Natural *
NaturalLor(Natural * a,Natural * b)651 NaturalLor (Natural *a, Natural *b)
652 {
653 ENTER ();
654 digit *at, *bt, *rt;
655 Natural *result;
656 int alength;
657 int blength;
658
659 alength = length(a);
660 blength = length(b);
661 if (alength < blength)
662 {
663 result = a;
664 a = b;
665 b = result;
666 alength = length(a);
667 blength = length(b);
668 }
669 if (alength == 0)
670 RETURN (zero_natural);
671 result = AllocNatural (alength);
672 at = data(a);
673 bt = data(b);
674 rt = data(result);
675 alength -= blength;
676 while (blength--)
677 *rt++ = *at++ | *bt++;
678 while (alength--)
679 *rt++ = *at++;
680 RETURN (result);
681 }
682
683 Natural *
NaturalCompliment(Natural * a,int len)684 NaturalCompliment (Natural *a, int len)
685 {
686 ENTER ();
687 digit *at, *rt;
688 Natural *result;
689 int resultlen;
690
691 resultlen = length (a);
692 at = data(a) + (resultlen-1);
693 while (resultlen > len && ~*at == 0)
694 {
695 resultlen--;
696 at--;
697 }
698 if (resultlen == 0)
699 RETURN (zero_natural);
700 if (resultlen > len)
701 len = resultlen;
702 result = AllocNatural (len);
703 rt = data(result) + (len-1);
704 while (len > resultlen)
705 {
706 *rt-- = ~0;
707 len--;
708 }
709 while (resultlen-- > 0)
710 *rt-- = ~*at--;
711 RETURN (result);
712 }
713
714 Natural *
NaturalNegate(Natural * n,int len)715 NaturalNegate (Natural *n, int len)
716 {
717 ENTER ();
718 RETURN (NaturalPlus (NaturalCompliment (n, len), one_natural));
719 }
720
721 Natural *
NaturalSqrt(Natural * n)722 NaturalSqrt (Natural *n)
723 {
724 ENTER ();
725 Natural *l, *h, *m, *rem;
726
727 l = two_natural;
728 h = NaturalDivide (n, two_natural, &rem);
729 while (NaturalLess (one_natural,
730 NaturalMinus (h, l)))
731 {
732 m = NaturalDivide (NaturalPlus (l, h), two_natural, &rem);
733 if (NaturalLess (NaturalTimes (m, m), n))
734 l = m;
735 else
736 h = m;
737 }
738 RETURN (h);
739 }
740
741 Natural *
NaturalFactor(Natural * n,Natural * max)742 NaturalFactor (Natural *n, Natural *max)
743 {
744 ENTER ();
745 Natural *v, *lim, *rem;
746
747 if (zerop (n))
748 RETURN(zero_natural);
749 if ((data(n)[0] & 1) == 0)
750 RETURN(two_natural);
751 lim = NaturalSqrt (n);
752 for (v = NewNatural (3);
753 !NaturalLess (lim, v);
754 v = NaturalPlus (v, two_natural))
755 {
756 (void) NaturalDivide (n, v, &rem);
757 if (zerop (rem))
758 RETURN (v);
759 if (aborting)
760 break;
761 if (max && NaturalLess (max, v))
762 RETURN (0);
763 }
764 RETURN (n);
765 }
766
767 Natural *
NaturalIntPow(Natural * n,int p)768 NaturalIntPow (Natural *n, int p)
769 {
770 ENTER ();
771 Natural *result;
772
773 result = one_natural;
774 while (p)
775 {
776 if (p & 1)
777 result = NaturalTimes (result, n);
778 p >>= 1;
779 if (p)
780 n = NaturalTimes (n, n);
781 if (aborting)
782 break;
783 }
784 RETURN (result);
785 }
786
787 Natural *
NaturalPow(Natural * n,Natural * p)788 NaturalPow (Natural *n, Natural *p)
789 {
790 ENTER ();
791 Natural *result;
792
793 result = one_natural;
794 while (!zerop (p))
795 {
796 if (data(p)[0] & 1)
797 result = NaturalTimes (result, n);
798 p = NaturalRsl(p, 1);
799 if (!zerop (p))
800 n = NaturalTimes (n, n);
801 if (aborting)
802 break;
803 }
804 RETURN (result);
805 }
806
807 #define evenp(n) ((zerop (n) || ((data(n)[0] & 1) == 0)))
808
809 Natural *
NaturalPowMod(Natural * n,Natural * p,Natural * m)810 NaturalPowMod (Natural *n, Natural *p, Natural *m)
811 {
812 ENTER ();
813 Natural *result;
814 Natural *rem;
815
816 result = one_natural;
817 while (!zerop (p))
818 {
819 if (!evenp (p))
820 (void) NaturalDivide (NaturalTimes (result, n), m, &result);
821 p = NaturalDivide (p, two_natural, &rem);
822 if (!zerop(p))
823 (void) NaturalDivide (NaturalTimes (n, n), m, &n);
824 if (aborting)
825 break;
826 }
827 RETURN (result);
828 }
829
830 static int
digit_width(digit d,int base)831 digit_width (digit d, int base)
832 {
833 int width = 1;
834 while (d >= base)
835 {
836 width++;
837 d /= base;
838 }
839 return width;
840 }
841
842 int
NaturalEstimateLength(Natural * a,int base)843 NaturalEstimateLength (Natural *a, int base)
844 {
845 if (length (a) == 0)
846 return 2;
847 return length(a) * digit_width (MAXDIGIT, base) + 1;
848 }
849
850 char *naturalBuffer;
851 int naturalBufferSize;
852
853 static char *
NaturalBottom(char * result,digit partial,int base,int digits,Bool fill)854 NaturalBottom (char *result, digit partial, int base, int digits, Bool fill)
855 {
856 digit dig;
857
858 do
859 {
860 dig = partial % base;
861 if (dig < 10)
862 dig = '0' + dig;
863 else
864 dig = 'a' + dig - 10;
865 *--result = dig;
866 digits--;
867 partial = partial / base;
868 } while (partial);
869 if (fill)
870 while (digits-- > 0)
871 *--result = '0';
872 return result;
873 }
874
875 static char *
NaturalSplit(char * result,Natural * a,Natural ** divisors,int base,int digits,Bool fill)876 NaturalSplit (char *result, Natural *a, Natural **divisors, int base, int digits, Bool fill)
877 {
878 ENTER ();
879 Natural *q, *r;
880 Bool rfill;
881
882 if (aborting)
883 return 0;
884 if (zerop (a))
885 {
886 if (fill)
887 while (digits--)
888 *--result = '0';
889 }
890 else if (!divisors[0])
891 {
892 result = NaturalBottom (result, data(a)[0], base, digits, fill);
893 }
894 else
895 {
896 q = NaturalDivide (a, divisors[0], &r);
897 digits = digits / 2;
898 divisors--;
899 rfill = True;
900 if (zerop (q))
901 rfill = fill;
902 result = NaturalSplit (result, r, divisors,
903 base, digits, rfill);
904 if (rfill)
905 result = NaturalSplit (result, q, divisors,
906 base, digits, fill);
907 }
908 EXIT ();
909 return result;
910 }
911
912 char *
NaturalSprint(char * result,Natural * a,int base,int * width)913 NaturalSprint (char *result, Natural *a, int base, int *width)
914 {
915 ENTER ();
916 int len;
917 double_digit max_base;
918 int digits;
919 digit *t;
920 Natural *divisor;
921 char *r;
922 digit partial;
923 int print_width;
924 Natural **divisors;
925 int ndivisors;
926 int idivisor;
927
928 if (!result)
929 {
930 /*
931 * Allocate temporary space for the string of digits
932 */
933 print_width = NaturalEstimateLength (a, base);
934 if (naturalBufferSize < print_width)
935 {
936 if (naturalBuffer)
937 free (naturalBuffer);
938 naturalBuffer = malloc (print_width);
939 if (!naturalBuffer)
940 {
941 naturalBufferSize = 0;
942 EXIT ();
943 return 0;
944 }
945 naturalBufferSize = print_width;
946 }
947 result = naturalBuffer + naturalBufferSize;
948 }
949 r = result;
950 *--r = '\0';
951 len = length (a);
952 if (len == 0)
953 {
954 *--r = '0';
955 if (width)
956 *width = 1;
957 EXIT ();
958 return r;
959 }
960 /*
961 * Compute the number of base digits which can be
962 * held in BASE
963 */
964 max_base = base;
965 digits = 0;
966 while (max_base <= BASE)
967 {
968 max_base *= base;
969 digits++;
970 }
971 max_base /= base;
972 t = 0;
973 divisor = 0;
974 if (max_base == BASE)
975 {
976 t = data(a);
977 while (len)
978 {
979 if (aborting)
980 {
981 r = 0;
982 break;
983 }
984 partial = *t++;
985 len--;
986 r = NaturalBottom (r, partial, base, digits, len != 0);
987 }
988 }
989 else
990 {
991 divisor = NewNatural ((unsigned) max_base);
992 divisors = 0;
993 ndivisors = 0;
994 idivisor = 0;
995 do
996 {
997 if (idivisor >= ndivisors - 1)
998 {
999 ndivisors += 128;
1000 if (divisors)
1001 divisors = realloc (divisors, ndivisors * sizeof (Natural *));
1002 else
1003 divisors = malloc (ndivisors * sizeof (Natural *));
1004 if (!divisors)
1005 return 0;
1006 }
1007 if (!idivisor)
1008 divisors[idivisor++] = 0;
1009 divisors[idivisor++] = divisor;
1010 divisor = NaturalTimes (divisor, divisor);
1011 digits = digits * 2;
1012 if (aborting)
1013 {
1014 r = 0;
1015 break;
1016 }
1017 } while (NaturalLess (divisor, a));
1018 if (!aborting)
1019 r = NaturalSplit (r, a, divisors + idivisor - 1, base, digits, False);
1020 free (divisors);
1021 }
1022 if (width && r)
1023 *width = (result - 1) - r;
1024 EXIT ();
1025 return r;
1026 }
1027
1028 DataType NaturalType = { 0, 0, "NaturalType" };
1029
1030 Natural *
AllocNatural(int size)1031 AllocNatural (int size)
1032 {
1033 Natural *result;
1034
1035 result = ALLOCATE (&NaturalType, sizeof (Natural) + size * sizeof (digit));
1036 result->length = size;
1037 return result;
1038 }
1039
1040 static Natural *
NewDoubleDigitNaturalReal(double_digit dd)1041 NewDoubleDigitNaturalReal (double_digit dd)
1042 {
1043 Natural *result;
1044 int len;
1045 double_digit temp;
1046 digit *d;
1047
1048 len = 0;
1049 temp = dd;
1050 while (temp) {
1051 len++;
1052 temp = DivBase (temp);
1053 }
1054 result = AllocNatural (len);
1055 temp = dd;
1056 d = data(result);
1057 while (temp) {
1058 *d++ = ModBase (temp);
1059 temp = DivBase (temp);
1060 }
1061 return result;
1062 }
1063
1064 #define NATURAL_CACHE_SIZE 8191
1065
1066 Natural *
NewDoubleDigitNatural(double_digit dd)1067 NewDoubleDigitNatural (double_digit dd)
1068 {
1069
1070 switch (dd) {
1071 case 0:
1072 return zero_natural;
1073 case 1:
1074 return one_natural;
1075 case 2:
1076 return two_natural;
1077 case MAX_NICKLE_INT:
1078 return max_int_natural;
1079 default:
1080 {
1081 digit l = ModBase(dd), u = DivBase(dd);
1082 unsigned c = l % NATURAL_CACHE_SIZE;
1083 Natural **re = (Natural **) DataCacheValues(naturalCache) + c;
1084 Natural *ret = *re;
1085 digit *d;
1086
1087 if (ret)
1088 {
1089 d = data(ret);
1090 if (l == d[0] && u == (ret->length == 1 ? 0 : d[1]))
1091 {
1092 REFERENCE (ret);
1093 return ret;
1094 }
1095 }
1096 ret = NewDoubleDigitNaturalReal (dd);
1097 *re = ret;
1098 return ret;
1099 }
1100 }
1101 }
1102
1103 Natural *
NewNatural(unsigned value)1104 NewNatural (unsigned value)
1105 {
1106 return NewDoubleDigitNatural ((double_digit) value);
1107 }
1108
1109 Natural *
NaturalRsl(Natural * v,int shift)1110 NaturalRsl (Natural *v, int shift)
1111 {
1112 ENTER ();
1113 Natural *r;
1114 digit *vt, *rt;
1115 digit d1, d2;
1116 int length;
1117 int dshift;
1118 int index, last;
1119
1120 if (v->length == 0)
1121 RETURN (zero_natural);
1122 #ifdef LLBASE2
1123 dshift = (shift >> LLBASE2);
1124 shift = (shift & (LBASE2 - 1));
1125 #else
1126 dshift = shift / LBASE2;
1127 shift = shift % LBASE2;
1128 #endif
1129 length = v->length - dshift;
1130 index = length;
1131 last = 1;
1132 if ((NaturalDigits(v)[v->length - 1] >> shift) == 0)
1133 {
1134 length--;
1135 last = 0;
1136 }
1137 if (length <= 0)
1138 RETURN (zero_natural);
1139 r = AllocNatural (length);
1140 rt = NaturalDigits (r);
1141 vt = NaturalDigits (v) + dshift;
1142 if (shift)
1143 {
1144 d2 = *vt++;
1145 while (--index)
1146 {
1147 d1 = d2;
1148 d2 = *vt++;
1149 *rt++ = (d1 >> shift) | (d2 << (LBASE2 - shift));
1150 }
1151 if (last)
1152 *rt++ = (d2 >> shift);
1153 }
1154 else
1155 {
1156 while (length--)
1157 {
1158 *rt++ = *vt++;
1159 }
1160 }
1161 RETURN (r);
1162 }
1163
1164 Natural *
NaturalLsl(Natural * v,int shift)1165 NaturalLsl (Natural *v, int shift)
1166 {
1167 ENTER ();
1168 Natural *r;
1169 digit *vt, *rt;
1170 digit d1, d2;
1171 int length;
1172 int dshift;
1173 int index;
1174 int last;
1175
1176 if (v->length == 0)
1177 RETURN (zero_natural);
1178 #ifdef LLBASE2
1179 dshift = (shift >> LLBASE2);
1180 shift = (shift & (LBASE2 - 1));
1181 #else
1182 dshift = shift / LBASE2;
1183 shift = shift % LBASE2;
1184 #endif
1185 length = v->length + dshift;
1186 index = v->length;
1187 last = 0;
1188 if (shift)
1189 {
1190 if ((NaturalDigits(v)[v->length - 1] >> (LBASE2 - shift)) != 0)
1191 {
1192 length++;
1193 last = 1;
1194 }
1195 }
1196 r = AllocNatural (length);
1197 rt = NaturalDigits (r);
1198 while (dshift--)
1199 *rt++ = 0;
1200 vt = NaturalDigits (v);
1201 if (shift)
1202 {
1203 d2 = *vt++;
1204 *rt++ = d2 << shift;
1205 while (--index)
1206 {
1207 d1 = d2;
1208 d2 = *vt++;
1209 *rt++ = (d1 >> (LBASE2 - shift)) | (d2 << shift);
1210 }
1211 if (last)
1212 *rt++ = (d2 >> (LBASE2 - shift));
1213 }
1214 else
1215 {
1216 while (index--)
1217 *rt++ = *vt++;
1218 }
1219 RETURN (r);
1220 }
1221
1222 Natural *
NaturalMask(Natural * v,int bits)1223 NaturalMask (Natural *v, int bits)
1224 {
1225 ENTER ();
1226 Natural *r;
1227 digit *vt, *rt;
1228 digit mask;
1229 int length;
1230
1231 #ifdef LLBASE2
1232 length = (bits + LBASE2) >> LLBASE2;
1233 mask = bits & (LBASE2 - 1);
1234 #else
1235 length = (bits + LBASE2) / LBASE2;
1236 mask = bits % LBASE2;
1237 #endif
1238 mask = (1 << mask) - 1;
1239 if (length > v->length)
1240 {
1241 length = v->length;
1242 mask = (digit) ~0;
1243 }
1244 while (length && (NaturalDigits(v)[length - 1] & mask) == 0)
1245 {
1246 length--;
1247 mask = (digit) ~0;
1248 }
1249 r = AllocNatural (length);
1250 rt = NaturalDigits (r);
1251 vt = NaturalDigits (v);
1252 if (length)
1253 {
1254 length--;
1255 while (length--)
1256 *rt++ = *vt++;
1257 *rt = *vt & mask;
1258 }
1259 RETURN (r);
1260 }
1261
1262 int
NaturalPowerOfTwo(Natural * v)1263 NaturalPowerOfTwo (Natural *v)
1264 {
1265 int bit;
1266 int l;
1267 digit *vt, last;
1268
1269 if (!v->length)
1270 return -1;
1271 vt = NaturalDigits(v);
1272 l = v->length - 1;
1273 while (l--)
1274 {
1275 if (*vt++ != 0)
1276 return -1;
1277 }
1278 last = *vt;
1279 if (last & (last - 1))
1280 return -1;
1281 bit = (v->length - 1) * LBASE2;
1282 while (!(last & 1))
1283 {
1284 bit++;
1285 last >>= 1;
1286 }
1287 return bit;
1288 }
1289
1290 void
NaturalDigitMultiply(Natural * a,digit i,Natural * result)1291 NaturalDigitMultiply (Natural *a, digit i, Natural *result)
1292 {
1293 result->length = DigitTimes (data(a), length(a), i,
1294 data(result));
1295 }
1296
1297 /*
1298 * subtract b from a in place with offset implied zeros to the
1299 * right of b. Return if a carry out occured
1300 */
1301
1302 digit
NaturalSubtractOffset(Natural * a,Natural * b,int offset)1303 NaturalSubtractOffset (Natural *a, Natural *b, int offset)
1304 {
1305 int index;
1306 digit carry;
1307 digit *at, *bt;
1308 digit av, bv;
1309 int len;
1310
1311 carry = 0;
1312 at = NaturalDigits(a) + offset;
1313 bt = NaturalDigits(b);
1314 index = a->length - offset;
1315 if (index > b->length)
1316 index = b->length;
1317 while (index--)
1318 {
1319 av = *at;
1320 bv = *bt++ + carry;
1321 if (bv)
1322 {
1323 carry = 0;
1324 if ((*at = av - bv) > av)
1325 carry = 1;
1326 }
1327 at++;
1328 }
1329 if (carry && a->length > b->length + offset)
1330 {
1331 *at = *at - carry;
1332 carry = 0;
1333 }
1334 len = a->length;
1335 at = NaturalDigits(a) + len;
1336 while (len > 0 && *--at == 0)
1337 len--;
1338 a->length = len;
1339 return carry;
1340 }
1341
1342 digit
NaturalSubtractOffsetReverse(Natural * a,Natural * b,int offset)1343 NaturalSubtractOffsetReverse (Natural *a, Natural *b, int offset)
1344 {
1345 int index;
1346 digit carry;
1347 digit *at, *bt;
1348 digit av, bv;
1349 int len;
1350
1351 carry = 0;
1352 at = NaturalDigits(a) + offset;
1353 bt = NaturalDigits(b);
1354 index = a->length - offset;
1355 if (index > b->length)
1356 index = b->length;
1357 while (index--)
1358 {
1359 av = *at + carry;
1360 bv = *bt++;
1361 if (bv)
1362 {
1363 carry = 0;
1364 if ((*at = bv - av) > bv)
1365 carry = 1;
1366 }
1367 at++;
1368 }
1369 if (carry && a->length > b->length + offset)
1370 {
1371 *at = carry;
1372 carry = 0;
1373 }
1374 len = a->length;
1375 at = NaturalDigits(a) + len;
1376 while (len > 0 && *--at == 0)
1377 len--;
1378 a->length = len;
1379 return carry;
1380 }
1381
1382 void
NaturalAddOffset(Natural * a,Natural * b,int offset)1383 NaturalAddOffset (Natural *a, Natural *b, int offset)
1384 {
1385 int index;
1386 digit carry;
1387 digit *at, *bt;
1388 digit av, bv;
1389
1390 carry = 0;
1391 at = NaturalDigits(a) + offset;
1392 bt = NaturalDigits(b);
1393 index = b->length;
1394 while (index--)
1395 {
1396 bv = *bt++ + carry;
1397 if (bv)
1398 {
1399 carry = 0;
1400 av = *at;
1401 if ((*at = av + bv) < av)
1402 carry = 1;
1403 }
1404 at++;
1405 }
1406 if (carry)
1407 *at = *at + carry;
1408 if (at == NaturalDigits(a) + a->length - 1)
1409 {
1410 while (a->length && *at == 0)
1411 {
1412 at--;
1413 a->length--;
1414 }
1415 }
1416 }
1417
1418 Bool
NaturalGreaterEqualOffset(Natural * a,Natural * b,int offset)1419 NaturalGreaterEqualOffset (Natural *a, Natural *b, int offset)
1420 {
1421 digit *ad, *bd;
1422 int index;
1423
1424 if (a->length > b->length + offset)
1425 return True;
1426 if (a->length < b->length + offset)
1427 return False;
1428 ad = NaturalDigits(a) + a->length - 1;
1429 bd = NaturalDigits(b) + b->length - 1;
1430 index = b->length;
1431 while (index--)
1432 {
1433 if (*ad > *bd)
1434 return True;
1435 if (*ad < *bd)
1436 return False;
1437 --ad;
1438 --bd;
1439 }
1440 return True;
1441 }
1442
1443 HashValue
NaturalHash(Natural * a)1444 NaturalHash (Natural *a)
1445 {
1446 return HashCrc32 ((unsigned char *) &a->length,
1447 sizeof (int) + sizeof (digit) * a->length);
1448 }
1449
1450 int
NaturalInit(void)1451 NaturalInit (void)
1452 {
1453 ENTER ();
1454 #ifndef LBASE10
1455 int max_tenpow, i;
1456 #endif
1457
1458 naturalCache = NewDataCache (NATURAL_CACHE_SIZE);
1459 zero_natural = NewDoubleDigitNaturalReal (0);
1460 MemAddRoot (zero_natural);
1461 one_natural = NewDoubleDigitNaturalReal (1);
1462 MemAddRoot (one_natural);
1463 two_natural = NewDoubleDigitNaturalReal (2);
1464 MemAddRoot (two_natural);
1465 max_int_natural = NewDoubleDigitNaturalReal (MAX_NICKLE_INT);
1466 MemAddRoot (max_int_natural);
1467 max_signed_digit_natural = NewDoubleDigitNaturalReal (MAX_NICKLE_SIGNED_DIGIT);
1468 MemAddRoot(max_signed_digit_natural);
1469 #ifndef LBASE10
1470 tenpow_digits = (int) floor (log10 ((double) MAX_NICKLE_INT));
1471 max_tenpow = 1;
1472 for (i = 0; i < tenpow_digits; i++)
1473 max_tenpow *= 10;
1474 #ifdef DEBUG
1475 printf ("max_tenpow: %d\n", max_tenpow);
1476 #endif
1477 max_tenpow_natural = NewNatural (max_tenpow);
1478 MemAddRoot (max_tenpow_natural);
1479 #endif
1480 EXIT ();
1481 return 1;
1482 }
1483
1484