2 Copyright (C) 1989-94 Massachusetts Institute of Technology
3 Portions copyright (C) 2004-2008 Slava Pestov
5 This material was developed by the Scheme project at the Massachusetts
6 Institute of Technology, Department of Electrical Engineering and
7 Computer Science. Permission to copy and modify this software, to
8 redistribute either the original software or a modified version, and
9 to use this software for any purpose is granted, subject to the
10 following restrictions and understandings.
12 1. Any copy made of this software must include this copyright notice
15 2. Users of this software agree to make their best efforts (a) to
16 return to the MIT Scheme project any improvements or extensions that
17 they make, so that these may be included in future releases; and (b)
18 to inform MIT of noteworthy uses of this software.
20 3. All materials developed as a consequence of the use of this
21 software shall duly acknowledge such use, in accordance with the usual
22 standards of acknowledging credit in academic research.
24 4. MIT has made no warrantee or representation that the operation of
25 this software will be error-free, and MIT is under no obligation to
26 provide any services, by way of maintenance, update, or otherwise.
28 5. In conjunction with products arising from the use of this material,
29 there shall be no use of the name of the Massachusetts Institute of
30 Technology nor of any adaptation thereof in any advertising,
31 promotional, or sales literature without prior written consent from
34 /* Changes for Scheme 48:
35 * - Converted to ANSI.
36 * - Added bitwise operations.
37 * - Added s48 to the beginning of all externally visible names.
38 * - Cached the bignum representations of -1, 0, and 1.
41 /* Changes for Factor:
42 * - Adapt bignumint.h for Factor memory manager
43 * - Add more bignum <-> C type conversions
44 * - Remove unused functions
45 * - Add local variable GC root recording
46 * - Remove s48 prefix from function names
47 * - Various fixes for Win64
49 * - Added bignum_gcd implementation
58 int factor_vm::bignum_equal_p(bignum* x, bignum* y) {
59 return ((BIGNUM_ZERO_P(x))
61 : ((!(BIGNUM_ZERO_P(y))) &&
62 ((BIGNUM_NEGATIVE_P(x)) ? (BIGNUM_NEGATIVE_P(y))
63 : (!(BIGNUM_NEGATIVE_P(y)))) &&
64 (bignum_equal_p_unsigned(x, y))));
67 enum bignum_comparison factor_vm::bignum_compare(bignum* x, bignum* y) {
68 return ((BIGNUM_ZERO_P(x)) ? ((BIGNUM_ZERO_P(y)) ? bignum_comparison_equal
69 : (BIGNUM_NEGATIVE_P(y))
70 ? bignum_comparison_greater
71 : bignum_comparison_less)
73 ? ((BIGNUM_NEGATIVE_P(x)) ? bignum_comparison_less
74 : bignum_comparison_greater)
75 : (BIGNUM_NEGATIVE_P(x))
76 ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_compare_unsigned(y, x))
77 : (bignum_comparison_less))
78 : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_comparison_greater)
79 : (bignum_compare_unsigned(x, y))));
82 /* Allocates memory */
83 bignum* factor_vm::bignum_add(bignum* x, bignum* y) {
85 (BIGNUM_ZERO_P(x)) ? (y) : (BIGNUM_ZERO_P(y))
87 : ((BIGNUM_NEGATIVE_P(x))
88 ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_add_unsigned(x, y, 1))
89 : (bignum_subtract_unsigned(y, x)))
90 : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_subtract_unsigned(x, y))
91 : (bignum_add_unsigned(x, y, 0)))));
94 /* Allocates memory */
95 bignum* factor_vm::bignum_subtract(bignum* x, bignum* y) {
96 return ((BIGNUM_ZERO_P(x))
97 ? ((BIGNUM_ZERO_P(y)) ? (y) : (bignum_new_sign(
98 y, (!(BIGNUM_NEGATIVE_P(y))))))
101 : ((BIGNUM_NEGATIVE_P(x))
102 ? ((BIGNUM_NEGATIVE_P(y))
103 ? (bignum_subtract_unsigned(y, x))
104 : (bignum_add_unsigned(x, y, 1)))
105 : ((BIGNUM_NEGATIVE_P(y))
106 ? (bignum_add_unsigned(x, y, 0))
107 : (bignum_subtract_unsigned(x, y))))));
111 bignum *factor_vm::bignum_square(bignum* x_)
113 return bignum_multiply(x_, x_);
116 /* Allocates memory */
117 bignum *factor_vm::bignum_square(bignum* x_)
119 data_root<bignum> x(x_, this);
121 bignum_length_type length = (BIGNUM_LENGTH (x));
122 bignum * z = (allot_bignum_zeroed ((length + length), 0));
124 bignum_digit_type * scan_z = BIGNUM_START_PTR (z);
125 bignum_digit_type * scan_x = BIGNUM_START_PTR (x);
126 bignum_digit_type * end_x = scan_x + length;
128 for (int i = 0; i < length; ++i) {
129 bignum_twodigit_type carry;
130 bignum_twodigit_type f = BIGNUM_REF (x, i);
131 bignum_digit_type *pz = scan_z + (i << 1);
132 bignum_digit_type *px = scan_x + i + 1;
135 *pz++ = carry & BIGNUM_DIGIT_MASK;
136 carry >>= BIGNUM_DIGIT_LENGTH;
137 BIGNUM_ASSERT (carry <= BIGNUM_DIGIT_MASK);
142 carry += *pz + *px++ * f;
143 *pz++ = carry & BIGNUM_DIGIT_MASK;
144 carry >>= BIGNUM_DIGIT_LENGTH;
145 BIGNUM_ASSERT (carry <= (BIGNUM_DIGIT_MASK << 1));
149 *pz++ = carry & BIGNUM_DIGIT_MASK;
150 carry >>= BIGNUM_DIGIT_LENGTH;
153 *pz += carry & BIGNUM_DIGIT_MASK;
154 BIGNUM_ASSERT ((carry >> BIGNUM_DIGIT_LENGTH) == 0);
156 return (bignum_trim (z));
160 /* Allocates memory */
161 bignum* factor_vm::bignum_multiply(bignum* x, bignum* y) {
165 return bignum_square(x);
169 bignum_length_type x_length = (BIGNUM_LENGTH(x));
170 bignum_length_type y_length = (BIGNUM_LENGTH(y));
171 int negative_p = ((BIGNUM_NEGATIVE_P(x)) ? (!(BIGNUM_NEGATIVE_P(y)))
172 : (BIGNUM_NEGATIVE_P(y)));
173 if (BIGNUM_ZERO_P(x))
175 if (BIGNUM_ZERO_P(y))
178 bignum_digit_type digit = (BIGNUM_REF(x, 0));
180 return (bignum_maybe_new_sign(y, negative_p));
181 if (digit < BIGNUM_RADIX_ROOT)
182 return (bignum_multiply_unsigned_small_factor(y, digit, negative_p));
185 bignum_digit_type digit = (BIGNUM_REF(y, 0));
187 return (bignum_maybe_new_sign(x, negative_p));
188 if (digit < BIGNUM_RADIX_ROOT)
189 return (bignum_multiply_unsigned_small_factor(x, digit, negative_p));
191 return (bignum_multiply_unsigned(x, y, negative_p));
194 /* Allocates memory */
195 void factor_vm::bignum_divide(bignum* numerator, bignum* denominator,
196 bignum** quotient, bignum** remainder) {
197 if (BIGNUM_ZERO_P(denominator)) {
198 divide_by_zero_error();
201 if (BIGNUM_ZERO_P(numerator)) {
202 (*quotient) = numerator;
203 (*remainder) = numerator;
205 int r_negative_p = (BIGNUM_NEGATIVE_P(numerator));
207 ((BIGNUM_NEGATIVE_P(denominator)) ? (!r_negative_p) : r_negative_p);
208 switch (bignum_compare_unsigned(numerator, denominator)) {
209 case bignum_comparison_equal: {
210 (*quotient) = (BIGNUM_ONE(q_negative_p));
211 (*remainder) = (BIGNUM_ZERO());
214 case bignum_comparison_less: {
215 (*quotient) = (BIGNUM_ZERO());
216 (*remainder) = numerator;
219 case bignum_comparison_greater: {
220 if ((BIGNUM_LENGTH(denominator)) == 1) {
221 bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
223 (*quotient) = (bignum_maybe_new_sign(numerator, q_negative_p));
224 (*remainder) = (BIGNUM_ZERO());
226 } else if (digit < BIGNUM_RADIX_ROOT) {
227 bignum_divide_unsigned_small_denominator(numerator, digit, quotient,
228 remainder, q_negative_p,
232 bignum_divide_unsigned_medium_denominator(
233 numerator, digit, quotient, remainder, q_negative_p,
238 bignum_divide_unsigned_large_denominator(
239 numerator, denominator, quotient, remainder, q_negative_p,
247 /* Allocates memory */
248 bignum* factor_vm::bignum_quotient(bignum* numerator, bignum* denominator) {
249 if (BIGNUM_ZERO_P(denominator)) {
250 divide_by_zero_error();
251 return (BIGNUM_OUT_OF_BAND);
253 if (BIGNUM_ZERO_P(numerator))
257 ((BIGNUM_NEGATIVE_P(denominator)) ? (!(BIGNUM_NEGATIVE_P(numerator)))
258 : (BIGNUM_NEGATIVE_P(numerator)));
259 switch (bignum_compare_unsigned(numerator, denominator)) {
260 case bignum_comparison_equal:
261 return (BIGNUM_ONE(q_negative_p));
262 case bignum_comparison_less:
263 return (BIGNUM_ZERO());
264 case bignum_comparison_greater:
265 default: /* to appease gcc -Wall */
268 if ((BIGNUM_LENGTH(denominator)) == 1) {
269 bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
271 return (bignum_maybe_new_sign(numerator, q_negative_p));
272 if (digit < BIGNUM_RADIX_ROOT)
273 bignum_divide_unsigned_small_denominator(
274 numerator, digit, ("ient), ((bignum**)0), q_negative_p, 0);
276 bignum_divide_unsigned_medium_denominator(
277 numerator, digit, ("ient), ((bignum**)0), q_negative_p, 0);
279 bignum_divide_unsigned_large_denominator(
280 numerator, denominator, ("ient), ((bignum**)0), q_negative_p,
288 /* Allocates memory */
289 bignum* factor_vm::bignum_remainder(bignum* numerator, bignum* denominator) {
290 if (BIGNUM_ZERO_P(denominator)) {
291 divide_by_zero_error();
292 return (BIGNUM_OUT_OF_BAND);
294 if (BIGNUM_ZERO_P(numerator))
296 switch (bignum_compare_unsigned(numerator, denominator)) {
297 case bignum_comparison_equal:
298 return (BIGNUM_ZERO());
299 case bignum_comparison_less:
301 case bignum_comparison_greater:
302 default: /* to appease gcc -Wall */
305 if ((BIGNUM_LENGTH(denominator)) == 1) {
306 bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
308 return (BIGNUM_ZERO());
309 if (digit < BIGNUM_RADIX_ROOT)
310 return (bignum_remainder_unsigned_small_denominator(
311 numerator, digit, (BIGNUM_NEGATIVE_P(numerator))));
312 bignum_divide_unsigned_medium_denominator(
313 numerator, digit, ((bignum**)0), (&remainder), 0,
314 (BIGNUM_NEGATIVE_P(numerator)));
316 bignum_divide_unsigned_large_denominator(
317 numerator, denominator, ((bignum**)0), (&remainder), 0,
318 (BIGNUM_NEGATIVE_P(numerator)));
324 /* cell_to_bignum, fixnum_to_bignum, long_long_to_bignum, ulong_long_to_bignum
326 /* Allocates memory */
327 #define FOO_TO_BIGNUM(name, type, stype, utype) \
328 bignum* factor_vm::name##_to_bignum(type n) { \
330 /* Special cases win when these small constants are cached. */ \
332 return (BIGNUM_ZERO()); \
334 return (BIGNUM_ONE(0)); \
335 if (n < (type) 0 && n == (type) - 1) \
336 return (BIGNUM_ONE(1)); \
338 utype accumulator = \
339 ((negative_p = (n < (type) 0)) ? ((type)(-(stype) n)) : n); \
340 if (accumulator < BIGNUM_RADIX) \
342 bignum* result = allot_bignum(1, negative_p); \
343 bignum_digit_type* scan = (BIGNUM_START_PTR(result)); \
344 *scan = (accumulator & BIGNUM_DIGIT_MASK); \
347 bignum_digit_type result_digits[BIGNUM_DIGITS_FOR(type)]; \
348 bignum_digit_type* end_digits = result_digits; \
350 (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); \
351 accumulator >>= BIGNUM_DIGIT_LENGTH; \
352 } while (accumulator != 0); \
354 (allot_bignum((end_digits - result_digits), negative_p)); \
355 bignum_digit_type* scan_digits = result_digits; \
356 bignum_digit_type* scan_result = (BIGNUM_START_PTR(result)); \
357 while (scan_digits < end_digits) \
358 (*scan_result++) = (*scan_digits++); \
364 FOO_TO_BIGNUM(cell, cell, fixnum, cell)
365 FOO_TO_BIGNUM(fixnum, fixnum, fixnum, cell)
366 FOO_TO_BIGNUM(long_long, int64_t, int64_t, uint64_t)
367 FOO_TO_BIGNUM(ulong_long, uint64_t, int64_t, uint64_t)
369 /* cannot allocate memory */
370 /* bignum_to_cell, fixnum_to_cell, long_long_to_cell, ulong_long_to_cell */
371 #define BIGNUM_TO_FOO(name, type, stype, utype) \
372 type factor_vm::bignum_to_##name(bignum* bn) { \
373 if (BIGNUM_ZERO_P(bn)) \
376 utype accumulator = 0; \
377 bignum_digit_type* start = (BIGNUM_START_PTR(bn)); \
378 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn))); \
379 while (start < scan) \
380 accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); \
381 return ((BIGNUM_NEGATIVE_P(bn)) ? ((type)(-(stype) accumulator)) \
386 BIGNUM_TO_FOO(cell, cell, fixnum, cell)
387 BIGNUM_TO_FOO(fixnum, fixnum, fixnum, cell)
388 BIGNUM_TO_FOO(long_long, int64_t, int64_t, uint64_t)
389 BIGNUM_TO_FOO(ulong_long, uint64_t, int64_t, uint64_t)
391 /* cannot allocate memory */
392 fixnum factor_vm::bignum_to_fixnum_strict(bignum* bn) {
393 fixnum len = BIGNUM_LENGTH(bn);
394 bignum_digit_type *digits = BIGNUM_START_PTR(bn);
395 if ((len == 1 && digits[0] > fixnum_max) || (len > 1)) {
396 general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag<bignum>(bn), false_object);
398 fixnum fix = bignum_to_fixnum(bn);
399 FACTOR_ASSERT(fix <= fixnum_max && fix >= fixnum_min);
403 #define DTB_WRITE_DIGIT(factor) \
405 significand *= (factor); \
406 digit = ((bignum_digit_type) significand); \
408 significand -= ((double)digit); \
411 #define inf std::numeric_limits<double>::infinity()
413 /* Allocates memory */
414 bignum* factor_vm::double_to_bignum(double x) {
415 if (x == inf || x == -inf || x != x)
416 return (BIGNUM_ZERO());
418 double significand = (frexp(x, (&exponent)));
420 return (BIGNUM_ZERO());
422 return (BIGNUM_ONE(x < 0));
424 significand = (-significand);
426 bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent));
427 bignum* result = (allot_bignum(length, (x < 0)));
428 bignum_digit_type* start = (BIGNUM_START_PTR(result));
429 bignum_digit_type* scan = (start + length);
430 bignum_digit_type digit;
431 int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
433 DTB_WRITE_DIGIT((fixnum)1 << odd_bits);
434 while (start < scan) {
435 if (significand == 0) {
440 DTB_WRITE_DIGIT(BIGNUM_RADIX);
446 #undef DTB_WRITE_DIGIT
450 int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
451 bignum_length_type length = (BIGNUM_LENGTH(x));
452 if (length != (BIGNUM_LENGTH(y)))
455 bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
456 bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
457 bignum_digit_type* end_x = (scan_x + length);
458 while (scan_x < end_x)
459 if ((*scan_x++) != (*scan_y++))
465 enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
467 bignum_length_type x_length = (BIGNUM_LENGTH(x));
468 bignum_length_type y_length = (BIGNUM_LENGTH(y));
469 if (x_length < y_length)
470 return (bignum_comparison_less);
471 if (x_length > y_length)
472 return (bignum_comparison_greater);
474 bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
475 bignum_digit_type* scan_x = (start_x + x_length);
476 bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
477 while (start_x < scan_x) {
478 bignum_digit_type digit_x = (*--scan_x);
479 bignum_digit_type digit_y = (*--scan_y);
480 if (digit_x < digit_y)
481 return (bignum_comparison_less);
482 if (digit_x > digit_y)
483 return (bignum_comparison_greater);
486 return (bignum_comparison_equal);
491 /* Allocates memory, fixed! */
492 bignum* factor_vm::bignum_add_unsigned(bignum* x_, bignum* y_, int negative_p) {
495 data_root<bignum> x(x_, this);
496 data_root<bignum> y(y_, this);
497 if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
501 bignum_length_type x_length = (BIGNUM_LENGTH(x));
503 bignum* r = (allot_bignum((x_length + 1), negative_p));
505 bignum_digit_type sum;
506 bignum_digit_type carry = 0;
507 bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
508 bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
510 bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
511 bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
512 while (scan_y < end_y) {
513 sum = ((*scan_x++) + (*scan_y++) + carry);
514 if (sum < BIGNUM_RADIX) {
518 (*scan_r++) = (sum - BIGNUM_RADIX);
524 bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
526 while (scan_x < end_x) {
527 sum = ((*scan_x++) + 1);
528 if (sum < BIGNUM_RADIX) {
533 (*scan_r++) = (sum - BIGNUM_RADIX);
535 while (scan_x < end_x)
536 (*scan_r++) = (*scan_x++);
542 return (bignum_shorten_length(r, x_length));
548 /* Allocates memory */
549 bignum* factor_vm::bignum_subtract_unsigned(bignum* x_, bignum* y_) {
551 data_root<bignum> x(x_, this);
552 data_root<bignum> y(y_, this);
555 switch (bignum_compare_unsigned(x.untagged(), y.untagged())) {
556 case bignum_comparison_equal:
557 return (BIGNUM_ZERO());
558 case bignum_comparison_less:
562 case bignum_comparison_greater:
567 bignum_length_type x_length = (BIGNUM_LENGTH(x));
569 bignum* r = (allot_bignum(x_length, negative_p));
571 bignum_digit_type difference;
572 bignum_digit_type borrow = 0;
573 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
574 bignum_digit_type* scan_r = BIGNUM_START_PTR(r);
576 bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
577 bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
578 while (scan_y < end_y) {
579 difference = (((*scan_x++) - (*scan_y++)) - borrow);
580 if (difference < 0) {
581 (*scan_r++) = (difference + BIGNUM_RADIX);
584 (*scan_r++) = difference;
590 bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
592 while (scan_x < end_x) {
593 difference = ((*scan_x++) - borrow);
595 (*scan_r++) = (difference + BIGNUM_RADIX);
597 (*scan_r++) = difference;
602 BIGNUM_ASSERT(borrow == 0);
603 while (scan_x < end_x)
604 (*scan_r++) = (*scan_x++);
606 return (bignum_trim(r));
611 Maximum value for product_low or product_high:
612 ((R * R) + (R * (R - 2)) + (R - 1))
613 Maximum value for carry: ((R * (R - 1)) + (R - 1))
614 where R == BIGNUM_RADIX_ROOT */
616 /* Allocates memory */
617 bignum* factor_vm::bignum_multiply_unsigned(bignum* x_, bignum* y_,
620 data_root<bignum> x(x_, this);
621 data_root<bignum> y(y_, this);
623 if (BIGNUM_LENGTH(x) > BIGNUM_LENGTH(y)) {
627 bignum_digit_type carry;
628 bignum_digit_type y_digit_low;
629 bignum_digit_type y_digit_high;
630 bignum_digit_type x_digit_low;
631 bignum_digit_type x_digit_high;
632 bignum_digit_type product_low;
633 bignum_digit_type* scan_r;
634 bignum_digit_type* scan_y;
635 bignum_length_type x_length = BIGNUM_LENGTH(x);
636 bignum_length_type y_length = BIGNUM_LENGTH(y);
638 bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p));
640 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
641 bignum_digit_type* end_x = (scan_x + x_length);
642 bignum_digit_type* start_y = BIGNUM_START_PTR(y);
643 bignum_digit_type* end_y = (start_y + y_length);
644 bignum_digit_type* start_r = (BIGNUM_START_PTR(r));
645 #define x_digit x_digit_high
646 #define y_digit y_digit_high
647 #define product_high carry
648 while (scan_x < end_x) {
649 x_digit = (*scan_x++);
650 x_digit_low = (HD_LOW(x_digit));
651 x_digit_high = (HD_HIGH(x_digit));
654 scan_r = (start_r++);
655 while (scan_y < end_y) {
656 y_digit = (*scan_y++);
657 y_digit_low = (HD_LOW(y_digit));
658 y_digit_high = (HD_HIGH(y_digit));
660 ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry)));
662 ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) +
663 (HD_HIGH(product_low)) + (HD_HIGH(carry)));
664 (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
665 carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high)));
669 return (bignum_trim(r));
676 /* Allocates memory */
677 bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x_,
680 data_root<bignum> x(x_, this);
682 bignum_length_type length_x = (BIGNUM_LENGTH(x));
684 bignum* p = (allot_bignum((length_x + 1), negative_p));
686 bignum_destructive_copy(x.untagged(), p);
687 (BIGNUM_REF(p, length_x)) = 0;
688 bignum_destructive_scale_up(p, y);
689 return (bignum_trim(p));
692 void factor_vm::bignum_destructive_add(bignum* bn, bignum_digit_type n) {
693 bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
694 bignum_digit_type digit;
695 digit = ((*scan) + n);
696 if (digit < BIGNUM_RADIX) {
700 (*scan++) = (digit - BIGNUM_RADIX);
702 digit = ((*scan) + 1);
703 if (digit < BIGNUM_RADIX) {
707 (*scan++) = (digit - BIGNUM_RADIX);
711 void factor_vm::bignum_destructive_scale_up(bignum* bn,
712 bignum_digit_type factor) {
713 bignum_digit_type carry = 0;
714 bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
715 bignum_digit_type two_digits;
716 bignum_digit_type product_low;
717 #define product_high carry
718 bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bn)));
719 BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
721 two_digits = (*scan);
722 product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry)));
723 product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) +
725 (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
726 carry = (HD_HIGH(product_high));
728 /* A carry here would be an overflow, i.e. it would not fit.
729 Hopefully the callers allocate enough space that this will
732 BIGNUM_ASSERT(carry == 0);
739 /* For help understanding this algorithm, see:
740 Knuth, Donald E., "The Art of Computer Programming",
741 volume 2, "Seminumerical Algorithms"
742 section 4.3.1, "Multiple-Precision Arithmetic". */
744 /* Allocates memory */
745 void factor_vm::bignum_divide_unsigned_large_denominator(
746 bignum* numerator_, bignum* denominator_, bignum** quotient,
747 bignum** remainder, int q_negative_p, int r_negative_p) {
749 data_root<bignum> numerator(numerator_, this);
750 data_root<bignum> denominator(denominator_, this);
752 bignum_length_type length_n = ((BIGNUM_LENGTH(numerator)) + 1);
753 bignum_length_type length_d = (BIGNUM_LENGTH(denominator));
756 if (quotient != ((bignum**)0)) {
757 q_ = allot_bignum(length_n - length_d, q_negative_p);
759 q_ = BIGNUM_OUT_OF_BAND;
762 data_root<bignum> q(q_, this);
764 data_root<bignum> u(allot_bignum(length_n, r_negative_p), this);
767 BIGNUM_ASSERT(length_d > 1);
769 bignum_digit_type v1 = (BIGNUM_REF((denominator), (length_d - 1)));
770 while (v1 < (BIGNUM_RADIX / 2)) {
776 bignum_destructive_copy(numerator.untagged(), u.untagged());
777 (BIGNUM_REF(u, (length_n - 1))) = 0;
778 bignum_divide_unsigned_normalized(u.untagged(),
779 denominator.untagged(),
782 bignum* v = (allot_bignum(length_d, 0));
784 bignum_destructive_normalization(numerator.untagged(),
787 bignum_destructive_normalization(denominator.untagged(),
790 bignum_divide_unsigned_normalized(u.untagged(), v, q.untagged());
791 if (remainder != ((bignum**)0))
792 bignum_destructive_unnormalization(u.untagged(), shift);
796 q = bignum_trim(q.untagged());
799 u = bignum_trim(u.untagged());
801 if (quotient != ((bignum**)0))
802 (*quotient) = q.untagged();
804 if (remainder != ((bignum**)0))
805 (*remainder) = u.untagged();
810 void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
812 bignum_length_type u_length = (BIGNUM_LENGTH(u));
813 bignum_length_type v_length = (BIGNUM_LENGTH(v));
814 bignum_digit_type* u_start = (BIGNUM_START_PTR(u));
815 bignum_digit_type* u_scan = (u_start + u_length);
816 bignum_digit_type* u_scan_limit = (u_start + v_length);
817 bignum_digit_type* u_scan_start = (u_scan - v_length);
818 bignum_digit_type* v_start = (BIGNUM_START_PTR(v));
819 bignum_digit_type* v_end = (v_start + v_length);
820 bignum_digit_type* q_scan = NULL;
821 bignum_digit_type v1 = (v_end[-1]);
822 bignum_digit_type v2 = (v_end[-2]);
823 bignum_digit_type ph; /* high half of double-digit product */
824 bignum_digit_type pl; /* low half of double-digit product */
825 bignum_digit_type guess;
826 bignum_digit_type gh; /* high half-digit of guess */
827 bignum_digit_type ch; /* high half of double-digit comparand */
828 bignum_digit_type v2l = (HD_LOW(v2));
829 bignum_digit_type v2h = (HD_HIGH(v2));
830 bignum_digit_type cl; /* low half of double-digit comparand */
831 #define gl ph /* low half-digit of guess */
834 bignum_digit_type gm; /* memory loc for reference parameter */
835 if (q != BIGNUM_OUT_OF_BAND)
836 q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q)));
837 while (u_scan_limit < u_scan) {
841 (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
842 guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
844 ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm)));
848 ch = ((u_scan[-1]) + v1);
849 guess = (BIGNUM_RADIX - 1);
852 /* product = (guess * v2); */
853 gl = (HD_LOW(guess));
854 gh = (HD_HIGH(guess));
856 ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl)));
857 pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))));
858 ph = ((v2h * gh) + (HD_HIGH(ph)));
859 /* if (comparand >= product) */
860 if ((ch > ph) || ((ch == ph) && (cl >= pl)))
863 /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
865 /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
866 if (ch >= BIGNUM_RADIX)
869 qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start)));
870 if (q != BIGNUM_OUT_OF_BAND)
879 bignum_digit_type factor_vm::bignum_divide_subtract(
880 bignum_digit_type* v_start, bignum_digit_type* v_end,
881 bignum_digit_type guess, bignum_digit_type* u_start) {
882 bignum_digit_type* v_scan = v_start;
883 bignum_digit_type* u_scan = u_start;
884 bignum_digit_type carry = 0;
888 bignum_digit_type gl = (HD_LOW(guess));
889 bignum_digit_type gh = (HD_HIGH(guess));
891 bignum_digit_type pl;
892 bignum_digit_type vl;
896 while (v_scan < v_end) {
900 pl = ((vl * gl) + (HD_LOW(carry)));
901 ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry)));
902 diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))));
904 (*u_scan++) = (diff + BIGNUM_RADIX);
905 carry = ((vh * gh) + (HD_HIGH(ph)) + 1);
908 carry = ((vh * gh) + (HD_HIGH(ph)));
913 diff = ((*u_scan) - carry);
915 (*u_scan) = (diff + BIGNUM_RADIX);
924 /* Subtraction generated carry, implying guess is one too large.
925 Add v back in to bring it back down. */
929 while (v_scan < v_end) {
930 bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
931 if (sum < BIGNUM_RADIX) {
935 (*u_scan++) = (sum - BIGNUM_RADIX);
940 bignum_digit_type sum = ((*u_scan) + carry);
941 (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
946 /* Allocates memory */
947 void factor_vm::bignum_divide_unsigned_medium_denominator(
948 bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
949 bignum** remainder, int q_negative_p, int r_negative_p) {
951 data_root<bignum> numerator(numerator_, this);
953 bignum_length_type length_n = (BIGNUM_LENGTH(numerator));
956 /* Because `bignum_digit_divide' requires a normalized denominator. */
957 while (denominator < (BIGNUM_RADIX / 2)) {
962 bignum_length_type length_q = (shift == 0) ? length_n : length_n + 1;
963 data_root<bignum> q(allot_bignum(length_q, q_negative_p), this);
965 bignum_destructive_copy(numerator.untagged(), q.untagged());
967 bignum_destructive_normalization(numerator.untagged(), q.untagged(), shift);
970 bignum_digit_type r = 0;
971 bignum_digit_type* start = (BIGNUM_START_PTR(q));
972 bignum_digit_type* scan = (start + length_q);
973 bignum_digit_type qj;
975 while (start < scan) {
976 r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
980 q = bignum_trim(q.untagged());
982 if (remainder != ((bignum**)0)) {
986 (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
989 if (quotient != ((bignum**)0))
990 (*quotient) = q.untagged();
995 void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
997 bignum_digit_type digit;
998 bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
999 bignum_digit_type carry = 0;
1000 bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1001 bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1002 bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target)));
1003 int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
1004 bignum_digit_type mask = (((cell)1 << shift_right) - 1);
1005 while (scan_source < end_source) {
1006 digit = (*scan_source++);
1007 (*scan_target++) = (((digit & mask) << shift_left) | carry);
1008 carry = (digit >> shift_right);
1010 if (scan_target < end_target)
1011 (*scan_target) = carry;
1013 BIGNUM_ASSERT(carry == 0);
1017 void factor_vm::bignum_destructive_unnormalization(bignum* bn,
1019 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1020 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1021 bignum_digit_type digit;
1022 bignum_digit_type carry = 0;
1023 int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1024 bignum_digit_type mask = (((fixnum)1 << shift_right) - 1);
1025 while (start < scan) {
1027 (*scan) = ((digit >> shift_right) | carry);
1028 carry = ((digit & mask) << shift_left);
1030 BIGNUM_ASSERT(carry == 0);
1034 /* This is a reduced version of the division algorithm, applied to the
1035 case of dividing two bignum digits by one bignum digit. It is
1036 assumed that the numerator, denominator are normalized. */
1038 #define BDD_STEP(qn, j) \
1042 uj_uj1 = (HD_CONS(uj, (u[j + 1]))); \
1043 guess = (uj_uj1 / v1); \
1044 comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2]))); \
1046 guess = (BIGNUM_RADIX_ROOT - 1); \
1047 comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2]))); \
1049 while ((guess * v2) > comparand) { \
1051 comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH); \
1052 if (comparand >= BIGNUM_RADIX) \
1055 qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \
1058 bignum_digit_type factor_vm::bignum_digit_divide(
1059 bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
1060 bignum_digit_type* q) /* return value */
1062 bignum_digit_type guess;
1063 bignum_digit_type comparand;
1064 bignum_digit_type v1 = (HD_HIGH(v));
1065 bignum_digit_type v2 = (HD_LOW(v));
1066 bignum_digit_type uj;
1067 bignum_digit_type uj_uj1;
1068 bignum_digit_type q1;
1069 bignum_digit_type q2;
1070 bignum_digit_type u[4];
1075 } else if (ul == v) {
1080 (u[0]) = (HD_HIGH(uh));
1081 (u[1]) = (HD_LOW(uh));
1082 (u[2]) = (HD_HIGH(ul));
1083 (u[3]) = (HD_LOW(ul));
1088 (*q) = (HD_CONS(q1, q2));
1089 return (HD_CONS((u[2]), (u[3])));
1094 #define BDDS_MULSUB(vn, un, carry_in) \
1096 product = ((vn * guess) + carry_in); \
1097 diff = (un - (HD_LOW(product))); \
1099 un = (diff + BIGNUM_RADIX_ROOT); \
1100 carry = ((HD_HIGH(product)) + 1); \
1103 carry = (HD_HIGH(product)); \
1107 #define BDDS_ADD(vn, un, carry_in) \
1109 sum = (vn + un + carry_in); \
1110 if (sum < BIGNUM_RADIX_ROOT) { \
1114 un = (sum - BIGNUM_RADIX_ROOT); \
1119 bignum_digit_type factor_vm::bignum_digit_divide_subtract(
1120 bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
1121 bignum_digit_type* u) {
1123 bignum_digit_type product;
1124 bignum_digit_type diff;
1125 bignum_digit_type carry;
1126 BDDS_MULSUB(v2, (u[2]), 0);
1127 BDDS_MULSUB(v1, (u[1]), carry);
1130 diff = ((u[0]) - carry);
1132 (u[0]) = (diff + BIGNUM_RADIX);
1139 bignum_digit_type sum;
1140 bignum_digit_type carry;
1141 BDDS_ADD(v2, (u[2]), 0);
1142 BDDS_ADD(v1, (u[1]), carry);
1152 /* Allocates memory */
1153 void factor_vm::bignum_divide_unsigned_small_denominator(
1154 bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
1155 bignum** remainder, int q_negative_p, int r_negative_p) {
1156 data_root<bignum> numerator(numerator_, this);
1158 bignum* q_ = bignum_new_sign(numerator.untagged(), q_negative_p);
1159 data_root<bignum> q(q_, this);
1161 bignum_digit_type r = bignum_destructive_scale_down(q.untagged(), denominator);
1163 q = bignum_trim(q.untagged());
1165 if (remainder != ((bignum**)0))
1166 (*remainder) = bignum_digit_to_bignum(r, r_negative_p);
1168 (*quotient) = q.untagged();
1173 /* Given (denominator > 1), it is fairly easy to show that
1174 (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1175 that all digits are < BIGNUM_RADIX. */
1177 bignum_digit_type factor_vm::bignum_destructive_scale_down(
1178 bignum* bn, bignum_digit_type denominator) {
1179 bignum_digit_type numerator;
1180 bignum_digit_type remainder = 0;
1181 bignum_digit_type two_digits;
1182 #define quotient_high remainder
1183 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1184 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1185 BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1186 while (start < scan) {
1187 two_digits = (*--scan);
1188 numerator = (HD_CONS(remainder, (HD_HIGH(two_digits))));
1189 quotient_high = (numerator / denominator);
1190 numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits))));
1191 (*scan) = (HD_CONS(quotient_high, (numerator / denominator)));
1192 remainder = (numerator % denominator);
1195 #undef quotient_high
1198 /* Allocates memory */
1199 bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
1200 bignum* n, bignum_digit_type d, int negative_p) {
1201 bignum_digit_type two_digits;
1202 bignum_digit_type* start = (BIGNUM_START_PTR(n));
1203 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
1204 bignum_digit_type r = 0;
1205 BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT));
1206 while (start < scan) {
1207 two_digits = (*--scan);
1208 r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d),
1209 (HD_LOW(two_digits)))) %
1212 return (bignum_digit_to_bignum(r, negative_p));
1215 /* Allocates memory */
1216 bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
1219 return (BIGNUM_ZERO());
1221 bignum* result = (allot_bignum(1, negative_p));
1222 (BIGNUM_REF(result, 0)) = digit;
1227 /* Allocates memory */
1228 bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
1229 BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
1230 bignum* result = allot_uninitialized_array<bignum>(length + 1);
1231 BIGNUM_SET_NEGATIVE_P(result, negative_p);
1235 /* Allocates memory */
1236 bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
1238 bignum* result = allot_bignum(length, negative_p);
1239 bignum_digit_type* scan = (BIGNUM_START_PTR(result));
1240 bignum_digit_type* end = (scan + length);
1246 /* Allocates memory */
1247 bignum* factor_vm::bignum_shorten_length(bignum* bn,
1248 bignum_length_type length) {
1249 bignum_length_type current_length = (BIGNUM_LENGTH(bn));
1250 BIGNUM_ASSERT((length >= 0) || (length <= current_length));
1251 if (length < current_length) {
1252 bn = reallot_array(bn, length + 1);
1253 BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1258 /* Allocates memory */
1259 bignum* factor_vm::bignum_trim(bignum* bn) {
1260 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1261 bignum_digit_type* end = (start + (BIGNUM_LENGTH(bn)));
1262 bignum_digit_type* scan = end;
1263 while ((start <= scan) && ((*--scan) == 0))
1267 bignum_length_type length = (scan - start);
1268 bn = reallot_array(bn, length + 1);
1269 BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1276 /* Allocates memory */
1277 bignum* factor_vm::bignum_new_sign(bignum* x_, int negative_p) {
1278 data_root<bignum> x(x_, this);
1279 bignum* result = allot_bignum(BIGNUM_LENGTH(x), negative_p);
1280 bignum_destructive_copy(x.untagged(), result);
1284 /* Allocates memory */
1285 bignum* factor_vm::bignum_maybe_new_sign(bignum* x_, int negative_p) {
1286 if ((BIGNUM_NEGATIVE_P(x_)) ? negative_p : (!negative_p))
1289 return bignum_new_sign(x_, negative_p);
1293 void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
1294 bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1295 bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1296 bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1297 while (scan_source < end_source)
1298 (*scan_target++) = (*scan_source++);
1303 * Added bitwise operations (and oddp).
1306 /* Allocates memory */
1307 bignum* factor_vm::bignum_bitwise_not(bignum* x_) {
1310 bignum_length_type size = BIGNUM_LENGTH(x_);
1311 int is_negative = BIGNUM_NEGATIVE_P(x_);
1312 data_root<bignum> x(x_, this);
1313 data_root<bignum> y(allot_bignum(size, is_negative ? 0 : 1), this);
1315 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
1316 bignum_digit_type* end_x = scan_x + size;
1317 bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
1320 while (scan_x < end_x) {
1322 *scan_y++ = BIGNUM_RADIX - 1;
1325 *scan_y++ = *scan_x++ - 1;
1331 while (scan_x < end_x) {
1332 if (*scan_x == (BIGNUM_RADIX - 1)) {
1336 *scan_y++ = *scan_x++ + 1;
1343 while (scan_x < end_x) {
1344 *scan_y++ = *scan_x++;
1348 bignum* ret = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
1349 bignum_destructive_copy(y.untagged(), ret);
1350 bignum_digit_type* ret_start = BIGNUM_START_PTR(ret);
1351 *(ret_start + size) = 1;
1354 return bignum_trim(y.untagged());
1358 /* Allocates memory */
1359 bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
1360 if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1361 return bignum_bitwise_not(
1362 bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1364 return bignum_magnitude_ash(arg1, n);
1371 /* Allocates memory */
1372 bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
1373 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1374 ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1375 : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1376 : (BIGNUM_NEGATIVE_P(arg2))
1377 ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1378 : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
1381 /* Allocates memory */
1382 bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
1383 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1384 ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1385 : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1386 : (BIGNUM_NEGATIVE_P(arg2))
1387 ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1388 : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
1391 /* Allocates memory */
1392 bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
1393 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1394 ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1395 : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1396 : (BIGNUM_NEGATIVE_P(arg2))
1397 ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1398 : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
1401 /* Allocates memory */
1402 /* ash for the magnitude */
1403 /* assume arg1 is a big number, n is a long */
1404 bignum* factor_vm::bignum_magnitude_ash(bignum* arg1_, fixnum n) {
1406 data_root<bignum> arg1(arg1_, this);
1408 bignum* result = NULL;
1409 bignum_digit_type* scan1;
1410 bignum_digit_type* scanr;
1411 bignum_digit_type* end;
1413 fixnum digit_offset, bit_offset;
1415 if (BIGNUM_ZERO_P(arg1))
1416 return arg1.untagged();
1419 digit_offset = n / BIGNUM_DIGIT_LENGTH;
1420 bit_offset = n % BIGNUM_DIGIT_LENGTH;
1422 result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
1423 BIGNUM_NEGATIVE_P(arg1));
1425 scanr = BIGNUM_START_PTR(result) + digit_offset;
1426 scan1 = BIGNUM_START_PTR(arg1);
1427 end = scan1 + BIGNUM_LENGTH(arg1);
1429 while (scan1 < end) {
1430 *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1431 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1433 *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1434 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1436 } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
1437 BIGNUM_DIGIT_LENGTH))) {
1438 result = BIGNUM_ZERO();
1440 digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1441 bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1443 result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
1444 BIGNUM_NEGATIVE_P(arg1));
1446 scanr = BIGNUM_START_PTR(result);
1447 scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
1448 end = scanr + BIGNUM_LENGTH(result) - 1;
1450 while (scanr < end) {
1451 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1452 *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
1456 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1457 } else if (n == 0) {
1458 result = arg1.untagged();
1461 return bignum_trim(result);
1464 /* Allocates memory */
1465 bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1_,
1467 data_root<bignum> arg1(arg1_, this);
1468 data_root<bignum> arg2(arg2_, this);
1470 bignum_length_type max_length;
1472 bignum_digit_type* scan1, *end1, digit1;
1473 bignum_digit_type* scan2, *end2, digit2;
1474 bignum_digit_type* scanr, *endr;
1477 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
1478 : BIGNUM_LENGTH(arg2);
1480 bignum* result = allot_bignum(max_length, 0);
1482 scanr = BIGNUM_START_PTR(result);
1483 scan1 = BIGNUM_START_PTR(arg1);
1484 scan2 = BIGNUM_START_PTR(arg2);
1485 endr = scanr + max_length;
1486 end1 = scan1 + BIGNUM_LENGTH(arg1);
1487 end2 = scan2 + BIGNUM_LENGTH(arg2);
1489 while (scanr < endr) {
1490 digit1 = (scan1 < end1) ? *scan1++ : 0;
1491 digit2 = (scan2 < end2) ? *scan2++ : 0;
1493 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1496 return bignum_trim(result);
1499 /* Allocates memory */
1500 bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1_,
1502 data_root<bignum> arg1(arg1_, this);
1503 data_root<bignum> arg2(arg2_, this);
1505 bignum_length_type max_length;
1507 bignum_digit_type* scan1, *end1, digit1;
1508 bignum_digit_type* scan2, *end2, digit2, carry2;
1509 bignum_digit_type* scanr, *endr;
1511 char neg_p = op == IOR_OP || op == XOR_OP;
1514 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
1515 : BIGNUM_LENGTH(arg2) + 1;
1517 bignum* result = allot_bignum(max_length, neg_p);
1519 scanr = BIGNUM_START_PTR(result);
1520 scan1 = BIGNUM_START_PTR(arg1);
1521 scan2 = BIGNUM_START_PTR(arg2);
1522 endr = scanr + max_length;
1523 end1 = scan1 + BIGNUM_LENGTH(arg1);
1524 end2 = scan2 + BIGNUM_LENGTH(arg2);
1528 while (scanr < endr) {
1529 digit1 = (scan1 < end1) ? *scan1++ : 0;
1530 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1532 if (digit2 < BIGNUM_RADIX)
1535 digit2 = (digit2 - BIGNUM_RADIX);
1540 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1545 bignum_negate_magnitude(result);
1547 return bignum_trim(result);
1550 /* Allocates memory */
1551 bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1_,
1553 data_root<bignum> arg1(arg1_, this);
1554 data_root<bignum> arg2(arg2_, this);
1556 bignum_length_type max_length;
1558 bignum_digit_type* scan1, *end1, digit1, carry1;
1559 bignum_digit_type* scan2, *end2, digit2, carry2;
1560 bignum_digit_type* scanr, *endr;
1562 char neg_p = op == AND_OP || op == IOR_OP;
1565 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
1566 : BIGNUM_LENGTH(arg2) + 1;
1568 bignum* result = allot_bignum(max_length, neg_p);
1570 scanr = BIGNUM_START_PTR(result);
1571 scan1 = BIGNUM_START_PTR(arg1);
1572 scan2 = BIGNUM_START_PTR(arg2);
1573 endr = scanr + max_length;
1574 end1 = scan1 + BIGNUM_LENGTH(arg1);
1575 end2 = scan2 + BIGNUM_LENGTH(arg2);
1580 while (scanr < endr) {
1581 digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1582 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1584 if (digit1 < BIGNUM_RADIX)
1587 digit1 = (digit1 - BIGNUM_RADIX);
1591 if (digit2 < BIGNUM_RADIX)
1594 digit2 = (digit2 - BIGNUM_RADIX);
1599 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1604 bignum_negate_magnitude(result);
1606 return bignum_trim(result);
1609 void factor_vm::bignum_negate_magnitude(bignum* arg) {
1610 bignum_digit_type* scan;
1611 bignum_digit_type* end;
1612 bignum_digit_type digit;
1613 bignum_digit_type carry;
1615 scan = BIGNUM_START_PTR(arg);
1616 end = scan + BIGNUM_LENGTH(arg);
1620 while (scan < end) {
1621 digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;
1623 if (digit < BIGNUM_RADIX)
1626 digit = (digit - BIGNUM_RADIX);
1634 /* Allocates memory */
1635 bignum* factor_vm::bignum_integer_length(bignum* x_) {
1636 data_root<bignum> x(x_, this);
1637 bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
1638 bignum_digit_type digit = (BIGNUM_REF(x, index));
1639 bignum_digit_type carry = 0;
1647 if (index < BIGNUM_RADIX_ROOT) {
1648 result = allot_bignum(1, 0);
1649 (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
1651 result = allot_bignum(2, 0);
1652 (BIGNUM_REF(result, 0)) = index;
1653 (BIGNUM_REF(result, 1)) = 0;
1654 bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
1655 bignum_destructive_add(result, carry);
1657 return (bignum_trim(result));
1660 /* Allocates memory */
1661 int factor_vm::bignum_logbitp(int shift, bignum* arg) {
1662 return ((BIGNUM_NEGATIVE_P(arg))
1663 ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
1664 : bignum_unsigned_logbitp(shift, arg));
1667 int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bn) {
1668 bignum_length_type len = (BIGNUM_LENGTH(bn));
1669 int index = shift / BIGNUM_DIGIT_LENGTH;
1672 bignum_digit_type digit = (BIGNUM_REF(bn, index));
1673 int p = shift % BIGNUM_DIGIT_LENGTH;
1674 bignum_digit_type mask = ((fixnum)1) << p;
1675 return (digit & mask) ? 1 : 0;
1679 /* Allocates memory. */
1680 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1682 data_root<bignum> a(a_, this);
1683 data_root<bignum> b(b_, this);
1685 /* Copies of a and b with that are both positive. */
1686 data_root<bignum> ac(bignum_maybe_new_sign(a.untagged(), 0), this);
1687 data_root<bignum> bc(bignum_maybe_new_sign(b.untagged(), 0), this);
1689 if (bignum_compare(ac.untagged(), bc.untagged()) == bignum_comparison_less) {
1693 while (BIGNUM_LENGTH(bc) != 0) {
1694 data_root<bignum> d(bignum_remainder(ac.untagged(), bc.untagged()), this);
1695 if (d.untagged() == BIGNUM_OUT_OF_BAND) {
1696 return d.untagged();
1701 return ac.untagged();
1704 /* Allocates memory */
1705 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1706 data_root<bignum> a(a_, this);
1707 data_root<bignum> b(b_, this);
1708 bignum_twodigit_type x, y, q, s, t, A, B, C, D;
1710 bignum_length_type size_a, size_b, size_c;
1711 bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
1712 bignum_digit_type* a_end, *b_end, *c_end;
1714 /* clone the bignums so we can modify them in-place */
1715 size_a = BIGNUM_LENGTH(a);
1716 data_root<bignum> c(allot_bignum(size_a, 0), this);
1717 // c = allot_bignum(size_a, 0);
1718 scan_a = BIGNUM_START_PTR(a);
1719 a_end = scan_a + size_a;
1720 scan_c = BIGNUM_START_PTR(c);
1721 while (scan_a < a_end)
1722 (*scan_c++) = (*scan_a++);
1724 size_b = BIGNUM_LENGTH(b);
1725 data_root<bignum> d(allot_bignum(size_b, 0), this);
1726 scan_b = BIGNUM_START_PTR(b);
1727 b_end = scan_b + size_b;
1728 scan_d = BIGNUM_START_PTR(d);
1729 while (scan_b < b_end)
1730 (*scan_d++) = (*scan_b++);
1733 /* Initial reduction: make sure that 0 <= b <= a. */
1734 if (bignum_compare(a.untagged(), b.untagged()) == bignum_comparison_less) {
1736 std::swap(size_a, size_b);
1739 while (size_a > 1) {
1740 nbits = log2(BIGNUM_REF(a, size_a - 1));
1741 x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
1742 (BIGNUM_REF(a, size_a - 2) >> nbits));
1743 y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
1745 ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
1748 /* inner loop of Lehmer's algorithm; */
1757 q = (x + (A - 1)) / (y - C);
1777 /* no progress; do a Euclidean step */
1779 return bignum_trim(a.untagged());
1781 data_root<bignum> e(bignum_trim(a.untagged()), this);
1782 data_root<bignum> f(bignum_trim(b.untagged()), this);
1783 c = bignum_remainder(e.untagged(), f.untagged());
1784 if (c.untagged() == BIGNUM_OUT_OF_BAND) {
1785 return c.untagged();
1789 scan_a = BIGNUM_START_PTR(a);
1790 scan_b = BIGNUM_START_PTR(b);
1791 a_end = scan_a + size_a;
1792 b_end = scan_b + size_b;
1793 while (scan_b < b_end)
1794 *(scan_a++) = *(scan_b++);
1795 while (scan_a < a_end)
1800 scan_b = BIGNUM_START_PTR(b);
1801 scan_c = BIGNUM_START_PTR(c);
1802 size_c = BIGNUM_LENGTH(c);
1803 c_end = scan_c + size_c;
1804 while (scan_c < c_end)
1805 *(scan_b++) = *(scan_c++);
1806 while (scan_b < b_end)
1814 a, b = A*b - B*a, D*a - C*b if k is odd
1815 a, b = A*a - B*b, D*b - C*a if k is even
1817 scan_a = BIGNUM_START_PTR(a);
1818 scan_b = BIGNUM_START_PTR(b);
1821 a_end = scan_a + size_a;
1822 b_end = scan_b + size_b;
1826 while (scan_b < b_end) {
1827 s += (A * *scan_b) - (B * *scan_a);
1828 t += (D * *scan_a++) - (C * *scan_b++);
1829 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1830 *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1831 s >>= BIGNUM_DIGIT_LENGTH;
1832 t >>= BIGNUM_DIGIT_LENGTH;
1834 while (scan_a < a_end) {
1836 t += (D * *scan_a++);
1837 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1838 //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1839 s >>= BIGNUM_DIGIT_LENGTH;
1840 t >>= BIGNUM_DIGIT_LENGTH;
1843 while (scan_b < b_end) {
1844 s += (A * *scan_a) - (B * *scan_b);
1845 t += (D * *scan_b++) - (C * *scan_a++);
1846 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1847 *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1848 s >>= BIGNUM_DIGIT_LENGTH;
1849 t >>= BIGNUM_DIGIT_LENGTH;
1851 while (scan_a < a_end) {
1853 t -= (C * *scan_a++);
1854 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1855 //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1856 s >>= BIGNUM_DIGIT_LENGTH;
1857 t >>= BIGNUM_DIGIT_LENGTH;
1860 BIGNUM_ASSERT(s == 0);
1861 BIGNUM_ASSERT(t == 0);
1863 // update size_a and size_b to remove any zeros at end
1864 while (size_a > 0 && *(--scan_a) == 0)
1866 while (size_b > 0 && *(--scan_b) == 0)
1869 BIGNUM_ASSERT(size_a >= size_b);
1872 /* a fits into a fixnum, so b must too */
1873 fixnum xx = bignum_to_fixnum(a.untagged());
1874 fixnum yy = bignum_to_fixnum(b.untagged());
1877 /* usual Euclidean algorithm for longs */
1884 return fixnum_to_bignum(xx);