Comment out alien.remote-control tests for now
[factor/jcg.git] / vm / bignum.c
blob1f4bc3ce7693f0435c41792bf884422eb5b6cf89
1 /* :tabSize=2:indentSize=2:noTabs=true:
3 Copyright (C) 1989-94 Massachusetts Institute of Technology
4 Portions copyright (C) 2004-2008 Slava Pestov
6 This material was developed by the Scheme project at the Massachusetts
7 Institute of Technology, Department of Electrical Engineering and
8 Computer Science. Permission to copy and modify this software, to
9 redistribute either the original software or a modified version, and
10 to use this software for any purpose is granted, subject to the
11 following restrictions and understandings.
13 1. Any copy made of this software must include this copyright notice
14 in full.
16 2. Users of this software agree to make their best efforts (a) to
17 return to the MIT Scheme project any improvements or extensions that
18 they make, so that these may be included in future releases; and (b)
19 to inform MIT of noteworthy uses of this software.
21 3. All materials developed as a consequence of the use of this
22 software shall duly acknowledge such use, in accordance with the usual
23 standards of acknowledging credit in academic research.
25 4. MIT has made no warrantee or representation that the operation of
26 this software will be error-free, and MIT is under no obligation to
27 provide any services, by way of maintenance, update, or otherwise.
29 5. In conjunction with products arising from the use of this material,
30 there shall be no use of the name of the Massachusetts Institute of
31 Technology nor of any adaptation thereof in any advertising,
32 promotional, or sales literature without prior written consent from
33 MIT in each case. */
35 /* Changes for Scheme 48:
36 * - Converted to ANSI.
37 * - Added bitwise operations.
38 * - Added s48 to the beginning of all externally visible names.
39 * - Cached the bignum representations of -1, 0, and 1.
42 /* Changes for Factor:
43 * - Adapt bignumint.h for Factor memory manager
44 * - Add more bignum <-> C type conversions
45 * - Remove unused functions
46 * - Add local variable GC root recording
47 * - Remove s48 prefix from function names
48 * - Various fixes for Win64
51 #include "master.h"
52 #include <limits.h>
53 #include <stdio.h>
54 #include <stdlib.h> /* abort */
55 #include <math.h>
57 /* Exports */
59 int
60 bignum_equal_p(bignum_type x, bignum_type y)
62 return
63 ((BIGNUM_ZERO_P (x))
64 ? (BIGNUM_ZERO_P (y))
65 : ((! (BIGNUM_ZERO_P (y)))
66 && ((BIGNUM_NEGATIVE_P (x))
67 ? (BIGNUM_NEGATIVE_P (y))
68 : (! (BIGNUM_NEGATIVE_P (y))))
69 && (bignum_equal_p_unsigned (x, y))));
72 enum bignum_comparison
73 bignum_compare(bignum_type x, bignum_type y)
75 return
76 ((BIGNUM_ZERO_P (x))
77 ? ((BIGNUM_ZERO_P (y))
78 ? bignum_comparison_equal
79 : (BIGNUM_NEGATIVE_P (y))
80 ? bignum_comparison_greater
81 : bignum_comparison_less)
82 : (BIGNUM_ZERO_P (y))
83 ? ((BIGNUM_NEGATIVE_P (x))
84 ? bignum_comparison_less
85 : bignum_comparison_greater)
86 : (BIGNUM_NEGATIVE_P (x))
87 ? ((BIGNUM_NEGATIVE_P (y))
88 ? (bignum_compare_unsigned (y, x))
89 : (bignum_comparison_less))
90 : ((BIGNUM_NEGATIVE_P (y))
91 ? (bignum_comparison_greater)
92 : (bignum_compare_unsigned (x, y))));
95 /* allocates memory */
96 bignum_type
97 bignum_add(bignum_type x, bignum_type y)
99 return
100 ((BIGNUM_ZERO_P (x))
101 ? (y)
102 : (BIGNUM_ZERO_P (y))
103 ? (x)
104 : ((BIGNUM_NEGATIVE_P (x))
105 ? ((BIGNUM_NEGATIVE_P (y))
106 ? (bignum_add_unsigned (x, y, 1))
107 : (bignum_subtract_unsigned (y, x)))
108 : ((BIGNUM_NEGATIVE_P (y))
109 ? (bignum_subtract_unsigned (x, y))
110 : (bignum_add_unsigned (x, y, 0)))));
113 /* allocates memory */
114 bignum_type
115 bignum_subtract(bignum_type x, bignum_type y)
117 return
118 ((BIGNUM_ZERO_P (x))
119 ? ((BIGNUM_ZERO_P (y))
120 ? (y)
121 : (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y))))))
122 : ((BIGNUM_ZERO_P (y))
123 ? (x)
124 : ((BIGNUM_NEGATIVE_P (x))
125 ? ((BIGNUM_NEGATIVE_P (y))
126 ? (bignum_subtract_unsigned (y, x))
127 : (bignum_add_unsigned (x, y, 1)))
128 : ((BIGNUM_NEGATIVE_P (y))
129 ? (bignum_add_unsigned (x, y, 0))
130 : (bignum_subtract_unsigned (x, y))))));
133 /* allocates memory */
134 bignum_type
135 bignum_multiply(bignum_type x, bignum_type y)
137 bignum_length_type x_length = (BIGNUM_LENGTH (x));
138 bignum_length_type y_length = (BIGNUM_LENGTH (y));
139 int negative_p =
140 ((BIGNUM_NEGATIVE_P (x))
141 ? (! (BIGNUM_NEGATIVE_P (y)))
142 : (BIGNUM_NEGATIVE_P (y)));
143 if (BIGNUM_ZERO_P (x))
144 return (x);
145 if (BIGNUM_ZERO_P (y))
146 return (y);
147 if (x_length == 1)
149 bignum_digit_type digit = (BIGNUM_REF (x, 0));
150 if (digit == 1)
151 return (bignum_maybe_new_sign (y, negative_p));
152 if (digit < BIGNUM_RADIX_ROOT)
153 return (bignum_multiply_unsigned_small_factor (y, digit, negative_p));
155 if (y_length == 1)
157 bignum_digit_type digit = (BIGNUM_REF (y, 0));
158 if (digit == 1)
159 return (bignum_maybe_new_sign (x, negative_p));
160 if (digit < BIGNUM_RADIX_ROOT)
161 return (bignum_multiply_unsigned_small_factor (x, digit, negative_p));
163 return (bignum_multiply_unsigned (x, y, negative_p));
166 /* allocates memory */
167 void
168 bignum_divide(bignum_type numerator, bignum_type denominator,
169 bignum_type * quotient, bignum_type * remainder)
171 if (BIGNUM_ZERO_P (denominator))
173 divide_by_zero_error(NULL);
174 return;
176 if (BIGNUM_ZERO_P (numerator))
178 (*quotient) = numerator;
179 (*remainder) = numerator;
181 else
183 int r_negative_p = (BIGNUM_NEGATIVE_P (numerator));
184 int q_negative_p =
185 ((BIGNUM_NEGATIVE_P (denominator)) ? (! r_negative_p) : r_negative_p);
186 switch (bignum_compare_unsigned (numerator, denominator))
188 case bignum_comparison_equal:
190 (*quotient) = (BIGNUM_ONE (q_negative_p));
191 (*remainder) = (BIGNUM_ZERO ());
192 break;
194 case bignum_comparison_less:
196 (*quotient) = (BIGNUM_ZERO ());
197 (*remainder) = numerator;
198 break;
200 case bignum_comparison_greater:
202 if ((BIGNUM_LENGTH (denominator)) == 1)
204 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
205 if (digit == 1)
207 (*quotient) =
208 (bignum_maybe_new_sign (numerator, q_negative_p));
209 (*remainder) = (BIGNUM_ZERO ());
210 break;
212 else if (digit < BIGNUM_RADIX_ROOT)
214 bignum_divide_unsigned_small_denominator
215 (numerator, digit,
216 quotient, remainder,
217 q_negative_p, r_negative_p);
218 break;
220 else
222 bignum_divide_unsigned_medium_denominator
223 (numerator, digit,
224 quotient, remainder,
225 q_negative_p, r_negative_p);
226 break;
229 bignum_divide_unsigned_large_denominator
230 (numerator, denominator,
231 quotient, remainder,
232 q_negative_p, r_negative_p);
233 break;
239 /* allocates memory */
240 bignum_type
241 bignum_quotient(bignum_type numerator, bignum_type denominator)
243 if (BIGNUM_ZERO_P (denominator))
245 divide_by_zero_error(NULL);
246 return (BIGNUM_OUT_OF_BAND);
248 if (BIGNUM_ZERO_P (numerator))
249 return numerator;
251 int q_negative_p =
252 ((BIGNUM_NEGATIVE_P (denominator))
253 ? (! (BIGNUM_NEGATIVE_P (numerator)))
254 : (BIGNUM_NEGATIVE_P (numerator)));
255 switch (bignum_compare_unsigned (numerator, denominator))
257 case bignum_comparison_equal:
258 return (BIGNUM_ONE (q_negative_p));
259 case bignum_comparison_less:
260 return (BIGNUM_ZERO ());
261 case bignum_comparison_greater:
262 default: /* to appease gcc -Wall */
264 bignum_type quotient;
265 if ((BIGNUM_LENGTH (denominator)) == 1)
267 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
268 if (digit == 1)
269 return (bignum_maybe_new_sign (numerator, q_negative_p));
270 if (digit < BIGNUM_RADIX_ROOT)
271 bignum_divide_unsigned_small_denominator
272 (numerator, digit,
273 (&quotient), ((bignum_type *) 0),
274 q_negative_p, 0);
275 else
276 bignum_divide_unsigned_medium_denominator
277 (numerator, digit,
278 (&quotient), ((bignum_type *) 0),
279 q_negative_p, 0);
281 else
282 bignum_divide_unsigned_large_denominator
283 (numerator, denominator,
284 (&quotient), ((bignum_type *) 0),
285 q_negative_p, 0);
286 return (quotient);
292 /* allocates memory */
293 bignum_type
294 bignum_remainder(bignum_type numerator, bignum_type denominator)
296 if (BIGNUM_ZERO_P (denominator))
298 divide_by_zero_error(NULL);
299 return (BIGNUM_OUT_OF_BAND);
301 if (BIGNUM_ZERO_P (numerator))
302 return numerator;
303 switch (bignum_compare_unsigned (numerator, denominator))
305 case bignum_comparison_equal:
306 return (BIGNUM_ZERO ());
307 case bignum_comparison_less:
308 return numerator;
309 case bignum_comparison_greater:
310 default: /* to appease gcc -Wall */
312 bignum_type remainder;
313 if ((BIGNUM_LENGTH (denominator)) == 1)
315 bignum_digit_type digit = (BIGNUM_REF (denominator, 0));
316 if (digit == 1)
317 return (BIGNUM_ZERO ());
318 if (digit < BIGNUM_RADIX_ROOT)
319 return
320 (bignum_remainder_unsigned_small_denominator
321 (numerator, digit, (BIGNUM_NEGATIVE_P (numerator))));
322 bignum_divide_unsigned_medium_denominator
323 (numerator, digit,
324 ((bignum_type *) 0), (&remainder),
325 0, (BIGNUM_NEGATIVE_P (numerator)));
327 else
328 bignum_divide_unsigned_large_denominator
329 (numerator, denominator,
330 ((bignum_type *) 0), (&remainder),
331 0, (BIGNUM_NEGATIVE_P (numerator)));
332 return (remainder);
337 #define FOO_TO_BIGNUM(name,type,utype) \
338 bignum_type name##_to_bignum(type n) \
340 int negative_p; \
341 bignum_digit_type result_digits [BIGNUM_DIGITS_FOR(type)]; \
342 bignum_digit_type * end_digits = result_digits; \
343 /* Special cases win when these small constants are cached. */ \
344 if (n == 0) return (BIGNUM_ZERO ()); \
345 if (n == 1) return (BIGNUM_ONE (0)); \
346 if (n < 0 && n == -1) return (BIGNUM_ONE (1)); \
348 utype accumulator = ((negative_p = (n < 0)) ? (-n) : n); \
349 do \
351 (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); \
352 accumulator >>= BIGNUM_DIGIT_LENGTH; \
354 while (accumulator != 0); \
357 bignum_type result = \
358 (allot_bignum ((end_digits - result_digits), negative_p)); \
359 bignum_digit_type * scan_digits = result_digits; \
360 bignum_digit_type * scan_result = (BIGNUM_START_PTR (result)); \
361 while (scan_digits < end_digits) \
362 (*scan_result++) = (*scan_digits++); \
363 return (result); \
367 /* all below allocate memory */
368 FOO_TO_BIGNUM(cell,CELL,CELL)
369 FOO_TO_BIGNUM(fixnum,F_FIXNUM,CELL)
370 FOO_TO_BIGNUM(long_long,s64,u64)
371 FOO_TO_BIGNUM(ulong_long,u64,u64)
373 #define BIGNUM_TO_FOO(name,type,utype) \
374 type bignum_to_##name(bignum_type bignum) \
376 if (BIGNUM_ZERO_P (bignum)) \
377 return (0); \
379 utype accumulator = 0; \
380 bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); \
381 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); \
382 while (start < scan) \
383 accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); \
384 return ((BIGNUM_NEGATIVE_P (bignum)) ? (-((type)accumulator)) : accumulator); \
388 /* all of the below allocate memory */
389 BIGNUM_TO_FOO(cell,CELL,CELL);
390 BIGNUM_TO_FOO(fixnum,F_FIXNUM,CELL);
391 BIGNUM_TO_FOO(long_long,s64,u64)
392 BIGNUM_TO_FOO(ulong_long,u64,u64)
394 double
395 bignum_to_double(bignum_type bignum)
397 if (BIGNUM_ZERO_P (bignum))
398 return (0);
400 double accumulator = 0;
401 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
402 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
403 while (start < scan)
404 accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan));
405 return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator);
409 #define DTB_WRITE_DIGIT(factor) \
411 significand *= (factor); \
412 digit = ((bignum_digit_type) significand); \
413 (*--scan) = digit; \
414 significand -= ((double) digit); \
417 /* allocates memory */
418 bignum_type
419 double_to_bignum(double x)
421 if (x == 1.0/0.0 || x == -1.0/0.0 || x != x) return (BIGNUM_ZERO ());
422 int exponent;
423 double significand = (frexp (x, (&exponent)));
424 if (exponent <= 0) return (BIGNUM_ZERO ());
425 if (exponent == 1) return (BIGNUM_ONE (x < 0));
426 if (significand < 0) significand = (-significand);
428 bignum_length_type length = (BIGNUM_BITS_TO_DIGITS (exponent));
429 bignum_type result = (allot_bignum (length, (x < 0)));
430 bignum_digit_type * start = (BIGNUM_START_PTR (result));
431 bignum_digit_type * scan = (start + length);
432 bignum_digit_type digit;
433 int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
434 if (odd_bits > 0)
435 DTB_WRITE_DIGIT ((F_FIXNUM)1 << odd_bits);
436 while (start < scan)
438 if (significand == 0)
440 while (start < scan)
441 (*--scan) = 0;
442 break;
444 DTB_WRITE_DIGIT (BIGNUM_RADIX);
446 return (result);
450 #undef DTB_WRITE_DIGIT
452 /* Comparisons */
455 bignum_equal_p_unsigned(bignum_type x, bignum_type y)
457 bignum_length_type length = (BIGNUM_LENGTH (x));
458 if (length != (BIGNUM_LENGTH (y)))
459 return (0);
460 else
462 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
463 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
464 bignum_digit_type * end_x = (scan_x + length);
465 while (scan_x < end_x)
466 if ((*scan_x++) != (*scan_y++))
467 return (0);
468 return (1);
472 enum bignum_comparison
473 bignum_compare_unsigned(bignum_type x, bignum_type y)
475 bignum_length_type x_length = (BIGNUM_LENGTH (x));
476 bignum_length_type y_length = (BIGNUM_LENGTH (y));
477 if (x_length < y_length)
478 return (bignum_comparison_less);
479 if (x_length > y_length)
480 return (bignum_comparison_greater);
482 bignum_digit_type * start_x = (BIGNUM_START_PTR (x));
483 bignum_digit_type * scan_x = (start_x + x_length);
484 bignum_digit_type * scan_y = ((BIGNUM_START_PTR (y)) + y_length);
485 while (start_x < scan_x)
487 bignum_digit_type digit_x = (*--scan_x);
488 bignum_digit_type digit_y = (*--scan_y);
489 if (digit_x < digit_y)
490 return (bignum_comparison_less);
491 if (digit_x > digit_y)
492 return (bignum_comparison_greater);
495 return (bignum_comparison_equal);
498 /* Addition */
500 /* allocates memory */
501 bignum_type
502 bignum_add_unsigned(bignum_type x, bignum_type y, int negative_p)
504 if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
506 bignum_type z = x;
507 x = y;
508 y = z;
511 bignum_length_type x_length = (BIGNUM_LENGTH (x));
513 REGISTER_BIGNUM(x);
514 REGISTER_BIGNUM(y);
515 bignum_type r = (allot_bignum ((x_length + 1), negative_p));
516 UNREGISTER_BIGNUM(y);
517 UNREGISTER_BIGNUM(x);
519 bignum_digit_type sum;
520 bignum_digit_type carry = 0;
521 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
522 bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
524 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
525 bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
526 while (scan_y < end_y)
528 sum = ((*scan_x++) + (*scan_y++) + carry);
529 if (sum < BIGNUM_RADIX)
531 (*scan_r++) = sum;
532 carry = 0;
534 else
536 (*scan_r++) = (sum - BIGNUM_RADIX);
537 carry = 1;
542 bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
543 if (carry != 0)
544 while (scan_x < end_x)
546 sum = ((*scan_x++) + 1);
547 if (sum < BIGNUM_RADIX)
549 (*scan_r++) = sum;
550 carry = 0;
551 break;
553 else
554 (*scan_r++) = (sum - BIGNUM_RADIX);
556 while (scan_x < end_x)
557 (*scan_r++) = (*scan_x++);
559 if (carry != 0)
561 (*scan_r) = 1;
562 return (r);
564 return (bignum_shorten_length (r, x_length));
568 /* Subtraction */
570 /* allocates memory */
571 bignum_type
572 bignum_subtract_unsigned(bignum_type x, bignum_type y)
574 int negative_p = 0;
575 switch (bignum_compare_unsigned (x, y))
577 case bignum_comparison_equal:
578 return (BIGNUM_ZERO ());
579 case bignum_comparison_less:
581 bignum_type z = x;
582 x = y;
583 y = z;
585 negative_p = 1;
586 break;
587 case bignum_comparison_greater:
588 negative_p = 0;
589 break;
592 bignum_length_type x_length = (BIGNUM_LENGTH (x));
594 REGISTER_BIGNUM(x);
595 REGISTER_BIGNUM(y);
596 bignum_type r = (allot_bignum (x_length, negative_p));
597 UNREGISTER_BIGNUM(y);
598 UNREGISTER_BIGNUM(x);
600 bignum_digit_type difference;
601 bignum_digit_type borrow = 0;
602 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
603 bignum_digit_type * scan_r = (BIGNUM_START_PTR (r));
605 bignum_digit_type * scan_y = (BIGNUM_START_PTR (y));
606 bignum_digit_type * end_y = (scan_y + (BIGNUM_LENGTH (y)));
607 while (scan_y < end_y)
609 difference = (((*scan_x++) - (*scan_y++)) - borrow);
610 if (difference < 0)
612 (*scan_r++) = (difference + BIGNUM_RADIX);
613 borrow = 1;
615 else
617 (*scan_r++) = difference;
618 borrow = 0;
623 bignum_digit_type * end_x = ((BIGNUM_START_PTR (x)) + x_length);
624 if (borrow != 0)
625 while (scan_x < end_x)
627 difference = ((*scan_x++) - borrow);
628 if (difference < 0)
629 (*scan_r++) = (difference + BIGNUM_RADIX);
630 else
632 (*scan_r++) = difference;
633 borrow = 0;
634 break;
637 BIGNUM_ASSERT (borrow == 0);
638 while (scan_x < end_x)
639 (*scan_r++) = (*scan_x++);
641 return (bignum_trim (r));
645 /* Multiplication
646 Maximum value for product_low or product_high:
647 ((R * R) + (R * (R - 2)) + (R - 1))
648 Maximum value for carry: ((R * (R - 1)) + (R - 1))
649 where R == BIGNUM_RADIX_ROOT */
651 /* allocates memory */
652 bignum_type
653 bignum_multiply_unsigned(bignum_type x, bignum_type y, int negative_p)
655 if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x)))
657 bignum_type z = x;
658 x = y;
659 y = z;
662 bignum_digit_type carry;
663 bignum_digit_type y_digit_low;
664 bignum_digit_type y_digit_high;
665 bignum_digit_type x_digit_low;
666 bignum_digit_type x_digit_high;
667 bignum_digit_type product_low;
668 bignum_digit_type * scan_r;
669 bignum_digit_type * scan_y;
670 bignum_length_type x_length = (BIGNUM_LENGTH (x));
671 bignum_length_type y_length = (BIGNUM_LENGTH (y));
673 REGISTER_BIGNUM(x);
674 REGISTER_BIGNUM(y);
675 bignum_type r =
676 (allot_bignum_zeroed ((x_length + y_length), negative_p));
677 UNREGISTER_BIGNUM(y);
678 UNREGISTER_BIGNUM(x);
680 bignum_digit_type * scan_x = (BIGNUM_START_PTR (x));
681 bignum_digit_type * end_x = (scan_x + x_length);
682 bignum_digit_type * start_y = (BIGNUM_START_PTR (y));
683 bignum_digit_type * end_y = (start_y + y_length);
684 bignum_digit_type * start_r = (BIGNUM_START_PTR (r));
685 #define x_digit x_digit_high
686 #define y_digit y_digit_high
687 #define product_high carry
688 while (scan_x < end_x)
690 x_digit = (*scan_x++);
691 x_digit_low = (HD_LOW (x_digit));
692 x_digit_high = (HD_HIGH (x_digit));
693 carry = 0;
694 scan_y = start_y;
695 scan_r = (start_r++);
696 while (scan_y < end_y)
698 y_digit = (*scan_y++);
699 y_digit_low = (HD_LOW (y_digit));
700 y_digit_high = (HD_HIGH (y_digit));
701 product_low =
702 ((*scan_r) +
703 (x_digit_low * y_digit_low) +
704 (HD_LOW (carry)));
705 product_high =
706 ((x_digit_high * y_digit_low) +
707 (x_digit_low * y_digit_high) +
708 (HD_HIGH (product_low)) +
709 (HD_HIGH (carry)));
710 (*scan_r++) =
711 (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
712 carry =
713 ((x_digit_high * y_digit_high) +
714 (HD_HIGH (product_high)));
716 (*scan_r) += carry;
718 return (bignum_trim (r));
719 #undef x_digit
720 #undef y_digit
721 #undef product_high
725 /* allocates memory */
726 bignum_type
727 bignum_multiply_unsigned_small_factor(bignum_type x, bignum_digit_type y,
728 int negative_p)
730 bignum_length_type length_x = (BIGNUM_LENGTH (x));
732 REGISTER_BIGNUM(x);
733 bignum_type p = (allot_bignum ((length_x + 1), negative_p));
734 UNREGISTER_BIGNUM(x);
736 bignum_destructive_copy (x, p);
737 (BIGNUM_REF (p, length_x)) = 0;
738 bignum_destructive_scale_up (p, y);
739 return (bignum_trim (p));
742 void
743 bignum_destructive_add(bignum_type bignum, bignum_digit_type n)
745 bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
746 bignum_digit_type digit;
747 digit = ((*scan) + n);
748 if (digit < BIGNUM_RADIX)
750 (*scan) = digit;
751 return;
753 (*scan++) = (digit - BIGNUM_RADIX);
754 while (1)
756 digit = ((*scan) + 1);
757 if (digit < BIGNUM_RADIX)
759 (*scan) = digit;
760 return;
762 (*scan++) = (digit - BIGNUM_RADIX);
766 void
767 bignum_destructive_scale_up(bignum_type bignum, bignum_digit_type factor)
769 bignum_digit_type carry = 0;
770 bignum_digit_type * scan = (BIGNUM_START_PTR (bignum));
771 bignum_digit_type two_digits;
772 bignum_digit_type product_low;
773 #define product_high carry
774 bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum)));
775 BIGNUM_ASSERT ((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
776 while (scan < end)
778 two_digits = (*scan);
779 product_low = ((factor * (HD_LOW (two_digits))) + (HD_LOW (carry)));
780 product_high =
781 ((factor * (HD_HIGH (two_digits))) +
782 (HD_HIGH (product_low)) +
783 (HD_HIGH (carry)));
784 (*scan++) = (HD_CONS ((HD_LOW (product_high)), (HD_LOW (product_low))));
785 carry = (HD_HIGH (product_high));
787 /* A carry here would be an overflow, i.e. it would not fit.
788 Hopefully the callers allocate enough space that this will
789 never happen.
791 BIGNUM_ASSERT (carry == 0);
792 return;
793 #undef product_high
796 /* Division */
798 /* For help understanding this algorithm, see:
799 Knuth, Donald E., "The Art of Computer Programming",
800 volume 2, "Seminumerical Algorithms"
801 section 4.3.1, "Multiple-Precision Arithmetic". */
803 /* allocates memory */
804 void
805 bignum_divide_unsigned_large_denominator(bignum_type numerator,
806 bignum_type denominator,
807 bignum_type * quotient,
808 bignum_type * remainder,
809 int q_negative_p,
810 int r_negative_p)
812 bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1);
813 bignum_length_type length_d = (BIGNUM_LENGTH (denominator));
815 REGISTER_BIGNUM(numerator);
816 REGISTER_BIGNUM(denominator);
818 bignum_type q =
819 ((quotient != ((bignum_type *) 0))
820 ? (allot_bignum ((length_n - length_d), q_negative_p))
821 : BIGNUM_OUT_OF_BAND);
823 REGISTER_BIGNUM(q);
824 bignum_type u = (allot_bignum (length_n, r_negative_p));
825 UNREGISTER_BIGNUM(q);
827 UNREGISTER_BIGNUM(denominator);
828 UNREGISTER_BIGNUM(numerator);
830 int shift = 0;
831 BIGNUM_ASSERT (length_d > 1);
833 bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1)));
834 while (v1 < (BIGNUM_RADIX / 2))
836 v1 <<= 1;
837 shift += 1;
840 if (shift == 0)
842 bignum_destructive_copy (numerator, u);
843 (BIGNUM_REF (u, (length_n - 1))) = 0;
844 bignum_divide_unsigned_normalized (u, denominator, q);
846 else
848 REGISTER_BIGNUM(numerator);
849 REGISTER_BIGNUM(denominator);
850 REGISTER_BIGNUM(u);
851 REGISTER_BIGNUM(q);
852 bignum_type v = (allot_bignum (length_d, 0));
853 UNREGISTER_BIGNUM(q);
854 UNREGISTER_BIGNUM(u);
855 UNREGISTER_BIGNUM(denominator);
856 UNREGISTER_BIGNUM(numerator);
858 bignum_destructive_normalization (numerator, u, shift);
859 bignum_destructive_normalization (denominator, v, shift);
860 bignum_divide_unsigned_normalized (u, v, q);
861 if (remainder != ((bignum_type *) 0))
862 bignum_destructive_unnormalization (u, shift);
865 REGISTER_BIGNUM(u);
866 if(q)
867 q = bignum_trim (q);
868 UNREGISTER_BIGNUM(u);
870 REGISTER_BIGNUM(q);
871 u = bignum_trim (u);
872 UNREGISTER_BIGNUM(q);
874 if (quotient != ((bignum_type *) 0))
875 (*quotient) = q;
877 if (remainder != ((bignum_type *) 0))
878 (*remainder) = u;
880 return;
883 void
884 bignum_divide_unsigned_normalized(bignum_type u, bignum_type v, bignum_type q)
886 bignum_length_type u_length = (BIGNUM_LENGTH (u));
887 bignum_length_type v_length = (BIGNUM_LENGTH (v));
888 bignum_digit_type * u_start = (BIGNUM_START_PTR (u));
889 bignum_digit_type * u_scan = (u_start + u_length);
890 bignum_digit_type * u_scan_limit = (u_start + v_length);
891 bignum_digit_type * u_scan_start = (u_scan - v_length);
892 bignum_digit_type * v_start = (BIGNUM_START_PTR (v));
893 bignum_digit_type * v_end = (v_start + v_length);
894 bignum_digit_type * q_scan = NULL;
895 bignum_digit_type v1 = (v_end[-1]);
896 bignum_digit_type v2 = (v_end[-2]);
897 bignum_digit_type ph; /* high half of double-digit product */
898 bignum_digit_type pl; /* low half of double-digit product */
899 bignum_digit_type guess;
900 bignum_digit_type gh; /* high half-digit of guess */
901 bignum_digit_type ch; /* high half of double-digit comparand */
902 bignum_digit_type v2l = (HD_LOW (v2));
903 bignum_digit_type v2h = (HD_HIGH (v2));
904 bignum_digit_type cl; /* low half of double-digit comparand */
905 #define gl ph /* low half-digit of guess */
906 #define uj pl
907 #define qj ph
908 bignum_digit_type gm; /* memory loc for reference parameter */
909 if (q != BIGNUM_OUT_OF_BAND)
910 q_scan = ((BIGNUM_START_PTR (q)) + (BIGNUM_LENGTH (q)));
911 while (u_scan_limit < u_scan)
913 uj = (*--u_scan);
914 if (uj != v1)
916 /* comparand =
917 (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
918 guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
919 cl = (u_scan[-2]);
920 ch = (bignum_digit_divide (uj, (u_scan[-1]), v1, (&gm)));
921 guess = gm;
923 else
925 cl = (u_scan[-2]);
926 ch = ((u_scan[-1]) + v1);
927 guess = (BIGNUM_RADIX - 1);
929 while (1)
931 /* product = (guess * v2); */
932 gl = (HD_LOW (guess));
933 gh = (HD_HIGH (guess));
934 pl = (v2l * gl);
935 ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH (pl)));
936 pl = (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl))));
937 ph = ((v2h * gh) + (HD_HIGH (ph)));
938 /* if (comparand >= product) */
939 if ((ch > ph) || ((ch == ph) && (cl >= pl)))
940 break;
941 guess -= 1;
942 /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
943 ch += v1;
944 /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
945 if (ch >= BIGNUM_RADIX)
946 break;
948 qj = (bignum_divide_subtract (v_start, v_end, guess, (--u_scan_start)));
949 if (q != BIGNUM_OUT_OF_BAND)
950 (*--q_scan) = qj;
952 return;
953 #undef gl
954 #undef uj
955 #undef qj
958 bignum_digit_type
959 bignum_divide_subtract(bignum_digit_type * v_start,
960 bignum_digit_type * v_end,
961 bignum_digit_type guess,
962 bignum_digit_type * u_start)
964 bignum_digit_type * v_scan = v_start;
965 bignum_digit_type * u_scan = u_start;
966 bignum_digit_type carry = 0;
967 if (guess == 0) return (0);
969 bignum_digit_type gl = (HD_LOW (guess));
970 bignum_digit_type gh = (HD_HIGH (guess));
971 bignum_digit_type v;
972 bignum_digit_type pl;
973 bignum_digit_type vl;
974 #define vh v
975 #define ph carry
976 #define diff pl
977 while (v_scan < v_end)
979 v = (*v_scan++);
980 vl = (HD_LOW (v));
981 vh = (HD_HIGH (v));
982 pl = ((vl * gl) + (HD_LOW (carry)));
983 ph = ((vl * gh) + (vh * gl) + (HD_HIGH (pl)) + (HD_HIGH (carry)));
984 diff = ((*u_scan) - (HD_CONS ((HD_LOW (ph)), (HD_LOW (pl)))));
985 if (diff < 0)
987 (*u_scan++) = (diff + BIGNUM_RADIX);
988 carry = ((vh * gh) + (HD_HIGH (ph)) + 1);
990 else
992 (*u_scan++) = diff;
993 carry = ((vh * gh) + (HD_HIGH (ph)));
996 if (carry == 0)
997 return (guess);
998 diff = ((*u_scan) - carry);
999 if (diff < 0)
1000 (*u_scan) = (diff + BIGNUM_RADIX);
1001 else
1003 (*u_scan) = diff;
1004 return (guess);
1006 #undef vh
1007 #undef ph
1008 #undef diff
1010 /* Subtraction generated carry, implying guess is one too large.
1011 Add v back in to bring it back down. */
1012 v_scan = v_start;
1013 u_scan = u_start;
1014 carry = 0;
1015 while (v_scan < v_end)
1017 bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
1018 if (sum < BIGNUM_RADIX)
1020 (*u_scan++) = sum;
1021 carry = 0;
1023 else
1025 (*u_scan++) = (sum - BIGNUM_RADIX);
1026 carry = 1;
1029 if (carry == 1)
1031 bignum_digit_type sum = ((*u_scan) + carry);
1032 (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
1034 return (guess - 1);
1037 /* allocates memory */
1038 void
1039 bignum_divide_unsigned_medium_denominator(bignum_type numerator,
1040 bignum_digit_type denominator,
1041 bignum_type * quotient,
1042 bignum_type * remainder,
1043 int q_negative_p,
1044 int r_negative_p)
1046 bignum_length_type length_n = (BIGNUM_LENGTH (numerator));
1047 bignum_length_type length_q;
1048 bignum_type q;
1049 int shift = 0;
1050 /* Because `bignum_digit_divide' requires a normalized denominator. */
1051 while (denominator < (BIGNUM_RADIX / 2))
1053 denominator <<= 1;
1054 shift += 1;
1056 if (shift == 0)
1058 length_q = length_n;
1060 REGISTER_BIGNUM(numerator);
1061 q = (allot_bignum (length_q, q_negative_p));
1062 UNREGISTER_BIGNUM(numerator);
1064 bignum_destructive_copy (numerator, q);
1066 else
1068 length_q = (length_n + 1);
1070 REGISTER_BIGNUM(numerator);
1071 q = (allot_bignum (length_q, q_negative_p));
1072 UNREGISTER_BIGNUM(numerator);
1074 bignum_destructive_normalization (numerator, q, shift);
1077 bignum_digit_type r = 0;
1078 bignum_digit_type * start = (BIGNUM_START_PTR (q));
1079 bignum_digit_type * scan = (start + length_q);
1080 bignum_digit_type qj;
1082 while (start < scan)
1084 r = (bignum_digit_divide (r, (*--scan), denominator, (&qj)));
1085 (*scan) = qj;
1088 q = bignum_trim (q);
1090 if (remainder != ((bignum_type *) 0))
1092 if (shift != 0)
1093 r >>= shift;
1095 REGISTER_BIGNUM(q);
1096 (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
1097 UNREGISTER_BIGNUM(q);
1100 if (quotient != ((bignum_type *) 0))
1101 (*quotient) = q;
1103 return;
1106 void
1107 bignum_destructive_normalization(bignum_type source, bignum_type target,
1108 int shift_left)
1110 bignum_digit_type digit;
1111 bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
1112 bignum_digit_type carry = 0;
1113 bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
1114 bignum_digit_type * end_source = (scan_source + (BIGNUM_LENGTH (source)));
1115 bignum_digit_type * end_target = (scan_target + (BIGNUM_LENGTH (target)));
1116 int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
1117 bignum_digit_type mask = (((CELL)1 << shift_right) - 1);
1118 while (scan_source < end_source)
1120 digit = (*scan_source++);
1121 (*scan_target++) = (((digit & mask) << shift_left) | carry);
1122 carry = (digit >> shift_right);
1124 if (scan_target < end_target)
1125 (*scan_target) = carry;
1126 else
1127 BIGNUM_ASSERT (carry == 0);
1128 return;
1131 void
1132 bignum_destructive_unnormalization(bignum_type bignum, int shift_right)
1134 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
1135 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
1136 bignum_digit_type digit;
1137 bignum_digit_type carry = 0;
1138 int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1139 bignum_digit_type mask = (((F_FIXNUM)1 << shift_right) - 1);
1140 while (start < scan)
1142 digit = (*--scan);
1143 (*scan) = ((digit >> shift_right) | carry);
1144 carry = ((digit & mask) << shift_left);
1146 BIGNUM_ASSERT (carry == 0);
1147 return;
1150 /* This is a reduced version of the division algorithm, applied to the
1151 case of dividing two bignum digits by one bignum digit. It is
1152 assumed that the numerator, denominator are normalized. */
1154 #define BDD_STEP(qn, j) \
1156 uj = (u[j]); \
1157 if (uj != v1) \
1159 uj_uj1 = (HD_CONS (uj, (u[j + 1]))); \
1160 guess = (uj_uj1 / v1); \
1161 comparand = (HD_CONS ((uj_uj1 % v1), (u[j + 2]))); \
1163 else \
1165 guess = (BIGNUM_RADIX_ROOT - 1); \
1166 comparand = (HD_CONS (((u[j + 1]) + v1), (u[j + 2]))); \
1168 while ((guess * v2) > comparand) \
1170 guess -= 1; \
1171 comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH); \
1172 if (comparand >= BIGNUM_RADIX) \
1173 break; \
1175 qn = (bignum_digit_divide_subtract (v1, v2, guess, (&u[j]))); \
1178 bignum_digit_type
1179 bignum_digit_divide(bignum_digit_type uh, bignum_digit_type ul,
1180 bignum_digit_type v,
1181 bignum_digit_type * q) /* return value */
1183 bignum_digit_type guess;
1184 bignum_digit_type comparand;
1185 bignum_digit_type v1 = (HD_HIGH (v));
1186 bignum_digit_type v2 = (HD_LOW (v));
1187 bignum_digit_type uj;
1188 bignum_digit_type uj_uj1;
1189 bignum_digit_type q1;
1190 bignum_digit_type q2;
1191 bignum_digit_type u [4];
1192 if (uh == 0)
1194 if (ul < v)
1196 (*q) = 0;
1197 return (ul);
1199 else if (ul == v)
1201 (*q) = 1;
1202 return (0);
1205 (u[0]) = (HD_HIGH (uh));
1206 (u[1]) = (HD_LOW (uh));
1207 (u[2]) = (HD_HIGH (ul));
1208 (u[3]) = (HD_LOW (ul));
1209 v1 = (HD_HIGH (v));
1210 v2 = (HD_LOW (v));
1211 BDD_STEP (q1, 0);
1212 BDD_STEP (q2, 1);
1213 (*q) = (HD_CONS (q1, q2));
1214 return (HD_CONS ((u[2]), (u[3])));
1217 #undef BDD_STEP
1219 #define BDDS_MULSUB(vn, un, carry_in) \
1221 product = ((vn * guess) + carry_in); \
1222 diff = (un - (HD_LOW (product))); \
1223 if (diff < 0) \
1225 un = (diff + BIGNUM_RADIX_ROOT); \
1226 carry = ((HD_HIGH (product)) + 1); \
1228 else \
1230 un = diff; \
1231 carry = (HD_HIGH (product)); \
1235 #define BDDS_ADD(vn, un, carry_in) \
1237 sum = (vn + un + carry_in); \
1238 if (sum < BIGNUM_RADIX_ROOT) \
1240 un = sum; \
1241 carry = 0; \
1243 else \
1245 un = (sum - BIGNUM_RADIX_ROOT); \
1246 carry = 1; \
1250 bignum_digit_type
1251 bignum_digit_divide_subtract(bignum_digit_type v1, bignum_digit_type v2,
1252 bignum_digit_type guess, bignum_digit_type * u)
1255 bignum_digit_type product;
1256 bignum_digit_type diff;
1257 bignum_digit_type carry;
1258 BDDS_MULSUB (v2, (u[2]), 0);
1259 BDDS_MULSUB (v1, (u[1]), carry);
1260 if (carry == 0)
1261 return (guess);
1262 diff = ((u[0]) - carry);
1263 if (diff < 0)
1264 (u[0]) = (diff + BIGNUM_RADIX);
1265 else
1267 (u[0]) = diff;
1268 return (guess);
1272 bignum_digit_type sum;
1273 bignum_digit_type carry;
1274 BDDS_ADD(v2, (u[2]), 0);
1275 BDDS_ADD(v1, (u[1]), carry);
1276 if (carry == 1)
1277 (u[0]) += 1;
1279 return (guess - 1);
1282 #undef BDDS_MULSUB
1283 #undef BDDS_ADD
1285 /* allocates memory */
1286 void
1287 bignum_divide_unsigned_small_denominator(bignum_type numerator,
1288 bignum_digit_type denominator,
1289 bignum_type * quotient,
1290 bignum_type * remainder,
1291 int q_negative_p,
1292 int r_negative_p)
1294 REGISTER_BIGNUM(numerator);
1295 bignum_type q = (bignum_new_sign (numerator, q_negative_p));
1296 UNREGISTER_BIGNUM(numerator);
1298 bignum_digit_type r = (bignum_destructive_scale_down (q, denominator));
1300 q = (bignum_trim (q));
1302 if (remainder != ((bignum_type *) 0))
1304 REGISTER_BIGNUM(q);
1305 (*remainder) = (bignum_digit_to_bignum (r, r_negative_p));
1306 UNREGISTER_BIGNUM(q);
1309 (*quotient) = q;
1311 return;
1314 /* Given (denominator > 1), it is fairly easy to show that
1315 (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1316 that all digits are < BIGNUM_RADIX. */
1318 bignum_digit_type
1319 bignum_destructive_scale_down(bignum_type bignum, bignum_digit_type denominator)
1321 bignum_digit_type numerator;
1322 bignum_digit_type remainder = 0;
1323 bignum_digit_type two_digits;
1324 #define quotient_high remainder
1325 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
1326 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum)));
1327 BIGNUM_ASSERT ((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1328 while (start < scan)
1330 two_digits = (*--scan);
1331 numerator = (HD_CONS (remainder, (HD_HIGH (two_digits))));
1332 quotient_high = (numerator / denominator);
1333 numerator = (HD_CONS ((numerator % denominator), (HD_LOW (two_digits))));
1334 (*scan) = (HD_CONS (quotient_high, (numerator / denominator)));
1335 remainder = (numerator % denominator);
1337 return (remainder);
1338 #undef quotient_high
1341 /* allocates memory */
1342 bignum_type
1343 bignum_remainder_unsigned_small_denominator(
1344 bignum_type n, bignum_digit_type d, int negative_p)
1346 bignum_digit_type two_digits;
1347 bignum_digit_type * start = (BIGNUM_START_PTR (n));
1348 bignum_digit_type * scan = (start + (BIGNUM_LENGTH (n)));
1349 bignum_digit_type r = 0;
1350 BIGNUM_ASSERT ((d > 1) && (d < BIGNUM_RADIX_ROOT));
1351 while (start < scan)
1353 two_digits = (*--scan);
1355 ((HD_CONS (((HD_CONS (r, (HD_HIGH (two_digits)))) % d),
1356 (HD_LOW (two_digits))))
1357 % d);
1359 return (bignum_digit_to_bignum (r, negative_p));
1362 /* allocates memory */
1363 bignum_type
1364 bignum_digit_to_bignum(bignum_digit_type digit, int negative_p)
1366 if (digit == 0)
1367 return (BIGNUM_ZERO ());
1368 else
1370 bignum_type result = (allot_bignum (1, negative_p));
1371 (BIGNUM_REF (result, 0)) = digit;
1372 return (result);
1376 /* allocates memory */
1377 bignum_type
1378 allot_bignum(bignum_length_type length, int negative_p)
1380 BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX));
1381 bignum_type result = allot_array_internal(BIGNUM_TYPE,length + 1);
1382 BIGNUM_SET_NEGATIVE_P (result, negative_p);
1383 return (result);
1386 /* allocates memory */
1387 bignum_type
1388 allot_bignum_zeroed(bignum_length_type length, int negative_p)
1390 bignum_type result = allot_bignum(length,negative_p);
1391 bignum_digit_type * scan = (BIGNUM_START_PTR (result));
1392 bignum_digit_type * end = (scan + length);
1393 while (scan < end)
1394 (*scan++) = 0;
1395 return (result);
1398 #define BIGNUM_REDUCE_LENGTH(source, length) \
1399 source = reallot_array(source,length + 1)
1401 /* allocates memory */
1402 bignum_type
1403 bignum_shorten_length(bignum_type bignum, bignum_length_type length)
1405 bignum_length_type current_length = (BIGNUM_LENGTH (bignum));
1406 BIGNUM_ASSERT ((length >= 0) || (length <= current_length));
1407 if (length < current_length)
1409 BIGNUM_REDUCE_LENGTH (bignum, length);
1410 BIGNUM_SET_NEGATIVE_P (bignum, (length != 0) && (BIGNUM_NEGATIVE_P (bignum)));
1412 return (bignum);
1415 /* allocates memory */
1416 bignum_type
1417 bignum_trim(bignum_type bignum)
1419 bignum_digit_type * start = (BIGNUM_START_PTR (bignum));
1420 bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum)));
1421 bignum_digit_type * scan = end;
1422 while ((start <= scan) && ((*--scan) == 0))
1424 scan += 1;
1425 if (scan < end)
1427 bignum_length_type length = (scan - start);
1428 BIGNUM_REDUCE_LENGTH (bignum, length);
1429 BIGNUM_SET_NEGATIVE_P (bignum, (length != 0) && (BIGNUM_NEGATIVE_P (bignum)));
1431 return (bignum);
1434 /* Copying */
1436 /* allocates memory */
1437 bignum_type
1438 bignum_new_sign(bignum_type bignum, int negative_p)
1440 REGISTER_BIGNUM(bignum);
1441 bignum_type result =
1442 (allot_bignum ((BIGNUM_LENGTH (bignum)), negative_p));
1443 UNREGISTER_BIGNUM(bignum);
1445 bignum_destructive_copy (bignum, result);
1446 return (result);
1449 /* allocates memory */
1450 bignum_type
1451 bignum_maybe_new_sign(bignum_type bignum, int negative_p)
1453 if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p))
1454 return (bignum);
1455 else
1457 bignum_type result =
1458 (allot_bignum ((BIGNUM_LENGTH (bignum)), negative_p));
1459 bignum_destructive_copy (bignum, result);
1460 return (result);
1464 void
1465 bignum_destructive_copy(bignum_type source, bignum_type target)
1467 bignum_digit_type * scan_source = (BIGNUM_START_PTR (source));
1468 bignum_digit_type * end_source =
1469 (scan_source + (BIGNUM_LENGTH (source)));
1470 bignum_digit_type * scan_target = (BIGNUM_START_PTR (target));
1471 while (scan_source < end_source)
1472 (*scan_target++) = (*scan_source++);
1473 return;
1477 * Added bitwise operations (and oddp).
1480 /* allocates memory */
1481 bignum_type
1482 bignum_bitwise_not(bignum_type x)
1484 return bignum_subtract(BIGNUM_ONE(1), x);
1487 /* allocates memory */
1488 bignum_type
1489 bignum_arithmetic_shift(bignum_type arg1, F_FIXNUM n)
1491 if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1492 return bignum_bitwise_not(bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1493 else
1494 return bignum_magnitude_ash(arg1, n);
1497 #define AND_OP 0
1498 #define IOR_OP 1
1499 #define XOR_OP 2
1501 /* allocates memory */
1502 bignum_type
1503 bignum_bitwise_and(bignum_type arg1, bignum_type arg2)
1505 return(
1506 (BIGNUM_NEGATIVE_P (arg1))
1507 ? (BIGNUM_NEGATIVE_P (arg2))
1508 ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1509 : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1510 : (BIGNUM_NEGATIVE_P (arg2))
1511 ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1512 : bignum_pospos_bitwise_op(AND_OP, arg1, arg2)
1516 /* allocates memory */
1517 bignum_type
1518 bignum_bitwise_ior(bignum_type arg1, bignum_type arg2)
1520 return(
1521 (BIGNUM_NEGATIVE_P (arg1))
1522 ? (BIGNUM_NEGATIVE_P (arg2))
1523 ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1524 : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1525 : (BIGNUM_NEGATIVE_P (arg2))
1526 ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1527 : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2)
1531 /* allocates memory */
1532 bignum_type
1533 bignum_bitwise_xor(bignum_type arg1, bignum_type arg2)
1535 return(
1536 (BIGNUM_NEGATIVE_P (arg1))
1537 ? (BIGNUM_NEGATIVE_P (arg2))
1538 ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1539 : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1540 : (BIGNUM_NEGATIVE_P (arg2))
1541 ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1542 : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2)
1546 /* allocates memory */
1547 /* ash for the magnitude */
1548 /* assume arg1 is a big number, n is a long */
1549 bignum_type
1550 bignum_magnitude_ash(bignum_type arg1, F_FIXNUM n)
1552 bignum_type result = NULL;
1553 bignum_digit_type *scan1;
1554 bignum_digit_type *scanr;
1555 bignum_digit_type *end;
1557 F_FIXNUM digit_offset,bit_offset;
1559 if (BIGNUM_ZERO_P (arg1)) return (arg1);
1561 if (n > 0) {
1562 digit_offset = n / BIGNUM_DIGIT_LENGTH;
1563 bit_offset = n % BIGNUM_DIGIT_LENGTH;
1565 REGISTER_BIGNUM(arg1);
1566 result = allot_bignum_zeroed (BIGNUM_LENGTH (arg1) + digit_offset + 1,
1567 BIGNUM_NEGATIVE_P(arg1));
1568 UNREGISTER_BIGNUM(arg1);
1570 scanr = BIGNUM_START_PTR (result) + digit_offset;
1571 scan1 = BIGNUM_START_PTR (arg1);
1572 end = scan1 + BIGNUM_LENGTH (arg1);
1574 while (scan1 < end) {
1575 *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1576 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1577 scanr++;
1578 *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1579 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1582 else if (n < 0
1583 && (-n >= (BIGNUM_LENGTH (arg1) * (bignum_length_type) BIGNUM_DIGIT_LENGTH)))
1584 result = BIGNUM_ZERO ();
1586 else if (n < 0) {
1587 digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1588 bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1590 REGISTER_BIGNUM(arg1);
1591 result = allot_bignum_zeroed (BIGNUM_LENGTH (arg1) - digit_offset,
1592 BIGNUM_NEGATIVE_P(arg1));
1593 UNREGISTER_BIGNUM(arg1);
1595 scanr = BIGNUM_START_PTR (result);
1596 scan1 = BIGNUM_START_PTR (arg1) + digit_offset;
1597 end = scanr + BIGNUM_LENGTH (result) - 1;
1599 while (scanr < end) {
1600 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset ;
1601 *scanr = (*scanr |
1602 *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) & BIGNUM_DIGIT_MASK;
1603 scanr++;
1605 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset ;
1607 else if (n == 0) result = arg1;
1609 return (bignum_trim (result));
1612 /* allocates memory */
1613 bignum_type
1614 bignum_pospos_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
1616 bignum_type result;
1617 bignum_length_type max_length;
1619 bignum_digit_type *scan1, *end1, digit1;
1620 bignum_digit_type *scan2, *end2, digit2;
1621 bignum_digit_type *scanr, *endr;
1623 max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2))
1624 ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2);
1626 REGISTER_BIGNUM(arg1);
1627 REGISTER_BIGNUM(arg2);
1628 result = allot_bignum(max_length, 0);
1629 UNREGISTER_BIGNUM(arg2);
1630 UNREGISTER_BIGNUM(arg1);
1632 scanr = BIGNUM_START_PTR(result);
1633 scan1 = BIGNUM_START_PTR(arg1);
1634 scan2 = BIGNUM_START_PTR(arg2);
1635 endr = scanr + max_length;
1636 end1 = scan1 + BIGNUM_LENGTH(arg1);
1637 end2 = scan2 + BIGNUM_LENGTH(arg2);
1639 while (scanr < endr) {
1640 digit1 = (scan1 < end1) ? *scan1++ : 0;
1641 digit2 = (scan2 < end2) ? *scan2++ : 0;
1642 *scanr++ = (op == AND_OP) ? digit1 & digit2 :
1643 (op == IOR_OP) ? digit1 | digit2 :
1644 digit1 ^ digit2;
1646 return bignum_trim(result);
1649 /* allocates memory */
1650 bignum_type
1651 bignum_posneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
1653 bignum_type result;
1654 bignum_length_type max_length;
1656 bignum_digit_type *scan1, *end1, digit1;
1657 bignum_digit_type *scan2, *end2, digit2, carry2;
1658 bignum_digit_type *scanr, *endr;
1660 char neg_p = op == IOR_OP || op == XOR_OP;
1662 max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1)
1663 ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2) + 1;
1665 REGISTER_BIGNUM(arg1);
1666 REGISTER_BIGNUM(arg2);
1667 result = allot_bignum(max_length, neg_p);
1668 UNREGISTER_BIGNUM(arg2);
1669 UNREGISTER_BIGNUM(arg1);
1671 scanr = BIGNUM_START_PTR(result);
1672 scan1 = BIGNUM_START_PTR(arg1);
1673 scan2 = BIGNUM_START_PTR(arg2);
1674 endr = scanr + max_length;
1675 end1 = scan1 + BIGNUM_LENGTH(arg1);
1676 end2 = scan2 + BIGNUM_LENGTH(arg2);
1678 carry2 = 1;
1680 while (scanr < endr) {
1681 digit1 = (scan1 < end1) ? *scan1++ : 0;
1682 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK)
1683 + carry2;
1685 if (digit2 < BIGNUM_RADIX)
1686 carry2 = 0;
1687 else
1689 digit2 = (digit2 - BIGNUM_RADIX);
1690 carry2 = 1;
1693 *scanr++ = (op == AND_OP) ? digit1 & digit2 :
1694 (op == IOR_OP) ? digit1 | digit2 :
1695 digit1 ^ digit2;
1698 if (neg_p)
1699 bignum_negate_magnitude(result);
1701 return bignum_trim(result);
1704 /* allocates memory */
1705 bignum_type
1706 bignum_negneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2)
1708 bignum_type result;
1709 bignum_length_type max_length;
1711 bignum_digit_type *scan1, *end1, digit1, carry1;
1712 bignum_digit_type *scan2, *end2, digit2, carry2;
1713 bignum_digit_type *scanr, *endr;
1715 char neg_p = op == AND_OP || op == IOR_OP;
1717 max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2))
1718 ? BIGNUM_LENGTH(arg1) + 1 : BIGNUM_LENGTH(arg2) + 1;
1720 REGISTER_BIGNUM(arg1);
1721 REGISTER_BIGNUM(arg2);
1722 result = allot_bignum(max_length, neg_p);
1723 UNREGISTER_BIGNUM(arg2);
1724 UNREGISTER_BIGNUM(arg1);
1726 scanr = BIGNUM_START_PTR(result);
1727 scan1 = BIGNUM_START_PTR(arg1);
1728 scan2 = BIGNUM_START_PTR(arg2);
1729 endr = scanr + max_length;
1730 end1 = scan1 + BIGNUM_LENGTH(arg1);
1731 end2 = scan2 + BIGNUM_LENGTH(arg2);
1733 carry1 = 1;
1734 carry2 = 1;
1736 while (scanr < endr) {
1737 digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1738 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1740 if (digit1 < BIGNUM_RADIX)
1741 carry1 = 0;
1742 else
1744 digit1 = (digit1 - BIGNUM_RADIX);
1745 carry1 = 1;
1748 if (digit2 < BIGNUM_RADIX)
1749 carry2 = 0;
1750 else
1752 digit2 = (digit2 - BIGNUM_RADIX);
1753 carry2 = 1;
1756 *scanr++ = (op == AND_OP) ? digit1 & digit2 :
1757 (op == IOR_OP) ? digit1 | digit2 :
1758 digit1 ^ digit2;
1761 if (neg_p)
1762 bignum_negate_magnitude(result);
1764 return bignum_trim(result);
1767 void
1768 bignum_negate_magnitude(bignum_type arg)
1770 bignum_digit_type *scan;
1771 bignum_digit_type *end;
1772 bignum_digit_type digit;
1773 bignum_digit_type carry;
1775 scan = BIGNUM_START_PTR(arg);
1776 end = scan + BIGNUM_LENGTH(arg);
1778 carry = 1;
1780 while (scan < end) {
1781 digit = (~*scan & BIGNUM_DIGIT_MASK) + carry;
1783 if (digit < BIGNUM_RADIX)
1784 carry = 0;
1785 else
1787 digit = (digit - BIGNUM_RADIX);
1788 carry = 1;
1791 *scan++ = digit;
1795 /* Allocates memory */
1796 bignum_type
1797 bignum_integer_length(bignum_type bignum)
1799 bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1);
1800 bignum_digit_type digit = (BIGNUM_REF (bignum, index));
1802 REGISTER_BIGNUM(bignum);
1803 bignum_type result = (allot_bignum (2, 0));
1804 UNREGISTER_BIGNUM(bignum);
1806 (BIGNUM_REF (result, 0)) = index;
1807 (BIGNUM_REF (result, 1)) = 0;
1808 bignum_destructive_scale_up (result, BIGNUM_DIGIT_LENGTH);
1809 while (digit > 1)
1811 bignum_destructive_add (result, ((bignum_digit_type) 1));
1812 digit >>= 1;
1814 return (bignum_trim (result));
1817 /* Allocates memory */
1819 bignum_logbitp(int shift, bignum_type arg)
1821 return((BIGNUM_NEGATIVE_P (arg))
1822 ? !bignum_unsigned_logbitp (shift, bignum_bitwise_not (arg))
1823 : bignum_unsigned_logbitp (shift,arg));
1827 bignum_unsigned_logbitp(int shift, bignum_type bignum)
1829 bignum_length_type len = (BIGNUM_LENGTH (bignum));
1830 bignum_digit_type digit;
1831 int index = shift / BIGNUM_DIGIT_LENGTH;
1832 int p;
1833 if (index >= len)
1834 return 0;
1835 digit = (BIGNUM_REF (bignum, index));
1836 p = shift % BIGNUM_DIGIT_LENGTH;
1837 return digit & (1 << p);
1840 /* Allocates memory */
1841 bignum_type
1842 digit_stream_to_bignum(unsigned int n_digits,
1843 unsigned int (*producer)(unsigned int),
1844 unsigned int radix,
1845 int negative_p)
1847 BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT));
1848 if (n_digits == 0)
1849 return (BIGNUM_ZERO ());
1850 if (n_digits == 1)
1852 F_FIXNUM digit = ((F_FIXNUM) ((*producer) (0)));
1853 return (fixnum_to_bignum (negative_p ? (- digit) : digit));
1856 bignum_length_type length;
1858 unsigned int radix_copy = radix;
1859 unsigned int log_radix = 0;
1860 while (radix_copy > 0)
1862 radix_copy >>= 1;
1863 log_radix += 1;
1865 /* This length will be at least as large as needed. */
1866 length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix));
1869 bignum_type result = (allot_bignum_zeroed (length, negative_p));
1870 while ((n_digits--) > 0)
1872 bignum_destructive_scale_up (result, ((bignum_digit_type) radix));
1873 bignum_destructive_add
1874 (result, ((bignum_digit_type) ((*producer) (n_digits))));
1876 return (bignum_trim (result));