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 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 bool bignum_fits_fixnum_p(bignum* bn) {
392 fixnum len = BIGNUM_LENGTH(bn);
397 bignum_digit_type dig = BIGNUM_START_PTR(bn)[0];
398 return (BIGNUM_NEGATIVE_P(bn) && dig <= -fixnum_min) ||
399 (!BIGNUM_NEGATIVE_P(bn) && dig <= fixnum_max);
402 cell bignum_maybe_to_fixnum(bignum* bn) {
403 if (bignum_fits_fixnum_p(bn))
404 return tag_fixnum(bignum_to_fixnum(bn));
405 return tag<bignum>(bn);
408 /* cannot allocate memory */
409 fixnum factor_vm::bignum_to_fixnum_strict(bignum* bn) {
411 if (!bignum_fits_fixnum_p(bn)) {
412 general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag<bignum>(bn), false_object);
414 fixnum fix = bignum_to_fixnum(bn);
415 FACTOR_ASSERT(fix <= fixnum_max && fix >= fixnum_min);
419 #define DTB_WRITE_DIGIT(factor) \
421 significand *= (factor); \
422 digit = ((bignum_digit_type) significand); \
424 significand -= ((double)digit); \
427 #define inf std::numeric_limits<double>::infinity()
429 /* Allocates memory */
430 bignum* factor_vm::double_to_bignum(double x) {
431 if (x == inf || x == -inf || x != x)
432 return (BIGNUM_ZERO());
434 double significand = (frexp(x, (&exponent)));
436 return (BIGNUM_ZERO());
438 return (BIGNUM_ONE(x < 0));
440 significand = (-significand);
442 bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent));
443 bignum* result = (allot_bignum(length, (x < 0)));
444 bignum_digit_type* start = (BIGNUM_START_PTR(result));
445 bignum_digit_type* scan = (start + length);
446 bignum_digit_type digit;
447 int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
449 DTB_WRITE_DIGIT((fixnum)1 << odd_bits);
450 while (start < scan) {
451 if (significand == 0) {
456 DTB_WRITE_DIGIT(BIGNUM_RADIX);
462 #undef DTB_WRITE_DIGIT
466 int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
467 bignum_length_type length = (BIGNUM_LENGTH(x));
468 if (length != (BIGNUM_LENGTH(y)))
471 bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
472 bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
473 bignum_digit_type* end_x = (scan_x + length);
474 while (scan_x < end_x)
475 if ((*scan_x++) != (*scan_y++))
481 enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
483 bignum_length_type x_length = (BIGNUM_LENGTH(x));
484 bignum_length_type y_length = (BIGNUM_LENGTH(y));
485 if (x_length < y_length)
486 return (bignum_comparison_less);
487 if (x_length > y_length)
488 return (bignum_comparison_greater);
490 bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
491 bignum_digit_type* scan_x = (start_x + x_length);
492 bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
493 while (start_x < scan_x) {
494 bignum_digit_type digit_x = (*--scan_x);
495 bignum_digit_type digit_y = (*--scan_y);
496 if (digit_x < digit_y)
497 return (bignum_comparison_less);
498 if (digit_x > digit_y)
499 return (bignum_comparison_greater);
502 return (bignum_comparison_equal);
507 /* Allocates memory */
508 bignum* factor_vm::bignum_add_unsigned(bignum* x_, bignum* y_, int negative_p) {
510 data_root<bignum> x(x_, this);
511 data_root<bignum> y(y_, this);
512 if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
516 bignum_length_type x_length = (BIGNUM_LENGTH(x));
518 bignum* r = (allot_bignum((x_length + 1), negative_p));
520 bignum_digit_type sum;
521 bignum_digit_type carry = 0;
522 bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
523 bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
525 bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
526 bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
527 while (scan_y < end_y) {
528 sum = ((*scan_x++) + (*scan_y++) + carry);
529 if (sum < BIGNUM_RADIX) {
533 (*scan_r++) = (sum - BIGNUM_RADIX);
539 bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
541 while (scan_x < end_x) {
542 sum = ((*scan_x++) + 1);
543 if (sum < BIGNUM_RADIX) {
548 (*scan_r++) = (sum - BIGNUM_RADIX);
550 while (scan_x < end_x)
551 (*scan_r++) = (*scan_x++);
557 return (bignum_shorten_length(r, x_length));
563 /* Allocates memory */
564 bignum* factor_vm::bignum_subtract_unsigned(bignum* x_, bignum* y_) {
566 data_root<bignum> x(x_, this);
567 data_root<bignum> y(y_, this);
570 switch (bignum_compare_unsigned(x.untagged(), y.untagged())) {
571 case bignum_comparison_equal:
572 return (BIGNUM_ZERO());
573 case bignum_comparison_less:
577 case bignum_comparison_greater:
582 bignum_length_type x_length = (BIGNUM_LENGTH(x));
584 bignum* r = (allot_bignum(x_length, negative_p));
586 bignum_digit_type difference;
587 bignum_digit_type borrow = 0;
588 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
589 bignum_digit_type* scan_r = BIGNUM_START_PTR(r);
591 bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
592 bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
593 while (scan_y < end_y) {
594 difference = (((*scan_x++) - (*scan_y++)) - borrow);
595 if (difference < 0) {
596 (*scan_r++) = (difference + BIGNUM_RADIX);
599 (*scan_r++) = difference;
605 bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
607 while (scan_x < end_x) {
608 difference = ((*scan_x++) - borrow);
610 (*scan_r++) = (difference + BIGNUM_RADIX);
612 (*scan_r++) = difference;
617 BIGNUM_ASSERT(borrow == 0);
618 while (scan_x < end_x)
619 (*scan_r++) = (*scan_x++);
621 return (bignum_trim(r));
626 Maximum value for product_low or product_high:
627 ((R * R) + (R * (R - 2)) + (R - 1))
628 Maximum value for carry: ((R * (R - 1)) + (R - 1))
629 where R == BIGNUM_RADIX_ROOT */
631 /* Allocates memory */
632 bignum* factor_vm::bignum_multiply_unsigned(bignum* x_, bignum* y_,
635 data_root<bignum> x(x_, this);
636 data_root<bignum> y(y_, this);
638 if (BIGNUM_LENGTH(y) > BIGNUM_LENGTH(x)) {
642 bignum_digit_type carry;
643 bignum_digit_type y_digit_low;
644 bignum_digit_type y_digit_high;
645 bignum_digit_type x_digit_low;
646 bignum_digit_type x_digit_high;
647 bignum_digit_type product_low;
648 bignum_digit_type* scan_r;
649 bignum_digit_type* scan_y;
650 bignum_length_type x_length = BIGNUM_LENGTH(x);
651 bignum_length_type y_length = BIGNUM_LENGTH(y);
653 bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p));
655 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
656 bignum_digit_type* end_x = (scan_x + x_length);
657 bignum_digit_type* start_y = BIGNUM_START_PTR(y);
658 bignum_digit_type* end_y = (start_y + y_length);
659 bignum_digit_type* start_r = (BIGNUM_START_PTR(r));
660 #define x_digit x_digit_high
661 #define y_digit y_digit_high
662 #define product_high carry
663 while (scan_x < end_x) {
664 x_digit = (*scan_x++);
665 x_digit_low = (HD_LOW(x_digit));
666 x_digit_high = (HD_HIGH(x_digit));
669 scan_r = (start_r++);
670 while (scan_y < end_y) {
671 y_digit = (*scan_y++);
672 y_digit_low = (HD_LOW(y_digit));
673 y_digit_high = (HD_HIGH(y_digit));
675 ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry)));
677 ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) +
678 (HD_HIGH(product_low)) + (HD_HIGH(carry)));
679 (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
680 carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high)));
684 return (bignum_trim(r));
691 /* Allocates memory */
692 bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x_,
695 data_root<bignum> x(x_, this);
697 bignum_length_type length_x = (BIGNUM_LENGTH(x));
699 bignum* p = (allot_bignum((length_x + 1), negative_p));
701 bignum_destructive_copy(x.untagged(), p);
702 (BIGNUM_REF(p, length_x)) = 0;
703 bignum_destructive_scale_up(p, y);
704 return (bignum_trim(p));
707 void factor_vm::bignum_destructive_add(bignum* bn, bignum_digit_type n) {
708 bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
709 bignum_digit_type digit;
710 digit = ((*scan) + n);
711 if (digit < BIGNUM_RADIX) {
715 (*scan++) = (digit - BIGNUM_RADIX);
717 digit = ((*scan) + 1);
718 if (digit < BIGNUM_RADIX) {
722 (*scan++) = (digit - BIGNUM_RADIX);
726 void factor_vm::bignum_destructive_scale_up(bignum* bn,
727 bignum_digit_type factor) {
728 bignum_digit_type carry = 0;
729 bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
730 bignum_digit_type two_digits;
731 bignum_digit_type product_low;
732 #define product_high carry
733 bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bn)));
734 BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
736 two_digits = (*scan);
737 product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry)));
738 product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) +
740 (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
741 carry = (HD_HIGH(product_high));
743 /* A carry here would be an overflow, i.e. it would not fit.
744 Hopefully the callers allocate enough space that this will
747 BIGNUM_ASSERT(carry == 0);
754 /* For help understanding this algorithm, see:
755 Knuth, Donald E., "The Art of Computer Programming",
756 volume 2, "Seminumerical Algorithms"
757 section 4.3.1, "Multiple-Precision Arithmetic". */
759 /* Allocates memory */
760 void factor_vm::bignum_divide_unsigned_large_denominator(
761 bignum* numerator_, bignum* denominator_,
762 bignum** quotient, bignum** remainder,
763 int q_negative_p, int r_negative_p) {
765 data_root<bignum> numerator(numerator_, this);
766 data_root<bignum> denominator(denominator_, this);
768 bignum_length_type length_n = BIGNUM_LENGTH(numerator) + 1;
769 bignum_length_type length_d = BIGNUM_LENGTH(denominator);
771 data_root<bignum> u(allot_bignum(length_n, r_negative_p), this);
774 BIGNUM_ASSERT(length_d > 1);
776 bignum_digit_type v1 = BIGNUM_REF(denominator.untagged(), length_d - 1);
777 while (v1 < (BIGNUM_RADIX / 2)) {
783 if (quotient != NULL) {
784 bignum *q_ = allot_bignum(length_n - length_d, q_negative_p);
785 data_root<bignum> q(q_, this);
788 bignum_destructive_copy(numerator.untagged(), u.untagged());
789 (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
790 bignum_divide_unsigned_normalized(u.untagged(),
791 denominator.untagged(),
794 bignum* v = allot_bignum(length_d, 0);
795 bignum_destructive_normalization(numerator.untagged(),
798 bignum_destructive_normalization(denominator.untagged(), v, shift);
799 bignum_divide_unsigned_normalized(u.untagged(), v, q.untagged());
800 if (remainder != NULL)
801 bignum_destructive_unnormalization(u.untagged(), shift);
804 q = bignum_trim(q.untagged());
805 *quotient = q.untagged();
809 bignum_destructive_copy(numerator.untagged(), u.untagged());
810 (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
811 bignum_divide_unsigned_normalized(u.untagged(),
812 denominator.untagged(),
815 bignum* v = allot_bignum(length_d, 0);
816 bignum_destructive_normalization(numerator.untagged(),
819 bignum_destructive_normalization(denominator.untagged(),
822 bignum_divide_unsigned_normalized(u.untagged(), v, NULL);
823 if (remainder != NULL)
824 bignum_destructive_unnormalization(u.untagged(), shift);
828 u = bignum_trim(u.untagged());
829 if (remainder != NULL)
830 *remainder = u.untagged();
833 void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
835 bignum_length_type u_length = (BIGNUM_LENGTH(u));
836 bignum_length_type v_length = (BIGNUM_LENGTH(v));
837 bignum_digit_type* u_start = (BIGNUM_START_PTR(u));
838 bignum_digit_type* u_scan = (u_start + u_length);
839 bignum_digit_type* u_scan_limit = (u_start + v_length);
840 bignum_digit_type* u_scan_start = (u_scan - v_length);
841 bignum_digit_type* v_start = (BIGNUM_START_PTR(v));
842 bignum_digit_type* v_end = (v_start + v_length);
843 bignum_digit_type* q_scan = NULL;
844 bignum_digit_type v1 = (v_end[-1]);
845 bignum_digit_type v2 = (v_end[-2]);
846 bignum_digit_type ph; /* high half of double-digit product */
847 bignum_digit_type pl; /* low half of double-digit product */
848 bignum_digit_type guess;
849 bignum_digit_type gh; /* high half-digit of guess */
850 bignum_digit_type ch; /* high half of double-digit comparand */
851 bignum_digit_type v2l = (HD_LOW(v2));
852 bignum_digit_type v2h = (HD_HIGH(v2));
853 bignum_digit_type cl; /* low half of double-digit comparand */
854 #define gl ph /* low half-digit of guess */
857 bignum_digit_type gm; /* memory loc for reference parameter */
858 if (q != BIGNUM_OUT_OF_BAND)
859 q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q)));
860 while (u_scan_limit < u_scan) {
864 (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
865 guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
867 ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm)));
871 ch = ((u_scan[-1]) + v1);
872 guess = (BIGNUM_RADIX - 1);
875 /* product = (guess * v2); */
876 gl = (HD_LOW(guess));
877 gh = (HD_HIGH(guess));
879 ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl)));
880 pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))));
881 ph = ((v2h * gh) + (HD_HIGH(ph)));
882 /* if (comparand >= product) */
883 if ((ch > ph) || ((ch == ph) && (cl >= pl)))
886 /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
888 /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
889 if (ch >= BIGNUM_RADIX)
892 qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start)));
893 if (q != BIGNUM_OUT_OF_BAND)
902 bignum_digit_type factor_vm::bignum_divide_subtract(
903 bignum_digit_type* v_start, bignum_digit_type* v_end,
904 bignum_digit_type guess, bignum_digit_type* u_start) {
905 bignum_digit_type* v_scan = v_start;
906 bignum_digit_type* u_scan = u_start;
907 bignum_digit_type carry = 0;
911 bignum_digit_type gl = (HD_LOW(guess));
912 bignum_digit_type gh = (HD_HIGH(guess));
914 bignum_digit_type pl;
915 bignum_digit_type vl;
919 while (v_scan < v_end) {
923 pl = ((vl * gl) + (HD_LOW(carry)));
924 ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry)));
925 diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))));
927 (*u_scan++) = (diff + BIGNUM_RADIX);
928 carry = ((vh * gh) + (HD_HIGH(ph)) + 1);
931 carry = ((vh * gh) + (HD_HIGH(ph)));
936 diff = ((*u_scan) - carry);
938 (*u_scan) = (diff + BIGNUM_RADIX);
947 /* Subtraction generated carry, implying guess is one too large.
948 Add v back in to bring it back down. */
952 while (v_scan < v_end) {
953 bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
954 if (sum < BIGNUM_RADIX) {
958 (*u_scan++) = (sum - BIGNUM_RADIX);
963 bignum_digit_type sum = ((*u_scan) + carry);
964 (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
969 /* Allocates memory */
970 void factor_vm::bignum_divide_unsigned_medium_denominator(
971 bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
972 bignum** remainder, int q_negative_p, int r_negative_p) {
974 data_root<bignum> numerator(numerator_, this);
976 bignum_length_type length_n = (BIGNUM_LENGTH(numerator));
979 /* Because `bignum_digit_divide' requires a normalized denominator. */
980 while (denominator < (BIGNUM_RADIX / 2)) {
985 bignum_length_type length_q = (shift == 0) ? length_n : length_n + 1;
986 data_root<bignum> q(allot_bignum(length_q, q_negative_p), this);
988 bignum_destructive_copy(numerator.untagged(), q.untagged());
990 bignum_destructive_normalization(numerator.untagged(), q.untagged(), shift);
993 bignum_digit_type r = 0;
994 bignum_digit_type* start = (BIGNUM_START_PTR(q));
995 bignum_digit_type* scan = (start + length_q);
996 bignum_digit_type qj;
998 while (start < scan) {
999 r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
1003 q = bignum_trim(q.untagged());
1005 if (remainder != ((bignum**)0)) {
1009 (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
1012 if (quotient != ((bignum**)0))
1013 (*quotient) = q.untagged();
1018 void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
1020 bignum_digit_type digit;
1021 bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1022 bignum_digit_type carry = 0;
1023 bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1024 bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1025 bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target)));
1026 int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
1027 bignum_digit_type mask = (((cell)1 << shift_right) - 1);
1028 while (scan_source < end_source) {
1029 digit = (*scan_source++);
1030 (*scan_target++) = (((digit & mask) << shift_left) | carry);
1031 carry = (digit >> shift_right);
1033 if (scan_target < end_target)
1034 (*scan_target) = carry;
1036 BIGNUM_ASSERT(carry == 0);
1040 void factor_vm::bignum_destructive_unnormalization(bignum* bn,
1042 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1043 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1044 bignum_digit_type digit;
1045 bignum_digit_type carry = 0;
1046 int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1047 bignum_digit_type mask = (((fixnum)1 << shift_right) - 1);
1048 while (start < scan) {
1050 (*scan) = ((digit >> shift_right) | carry);
1051 carry = ((digit & mask) << shift_left);
1053 BIGNUM_ASSERT(carry == 0);
1057 /* This is a reduced version of the division algorithm, applied to the
1058 case of dividing two bignum digits by one bignum digit. It is
1059 assumed that the numerator, denominator are normalized. */
1061 #define BDD_STEP(qn, j) \
1065 uj_uj1 = (HD_CONS(uj, (u[j + 1]))); \
1066 guess = (uj_uj1 / v1); \
1067 comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2]))); \
1069 guess = (BIGNUM_RADIX_ROOT - 1); \
1070 comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2]))); \
1072 while ((guess * v2) > comparand) { \
1074 comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH); \
1075 if (comparand >= BIGNUM_RADIX) \
1078 qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \
1081 bignum_digit_type factor_vm::bignum_digit_divide(
1082 bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
1083 bignum_digit_type* q) /* return value */
1085 bignum_digit_type guess;
1086 bignum_digit_type comparand;
1087 bignum_digit_type v1 = (HD_HIGH(v));
1088 bignum_digit_type v2 = (HD_LOW(v));
1089 bignum_digit_type uj;
1090 bignum_digit_type uj_uj1;
1091 bignum_digit_type q1;
1092 bignum_digit_type q2;
1093 bignum_digit_type u[4];
1098 } else if (ul == v) {
1103 (u[0]) = (HD_HIGH(uh));
1104 (u[1]) = (HD_LOW(uh));
1105 (u[2]) = (HD_HIGH(ul));
1106 (u[3]) = (HD_LOW(ul));
1111 (*q) = (HD_CONS(q1, q2));
1112 return (HD_CONS((u[2]), (u[3])));
1117 #define BDDS_MULSUB(vn, un, carry_in) \
1119 product = ((vn * guess) + carry_in); \
1120 diff = (un - (HD_LOW(product))); \
1122 un = (diff + BIGNUM_RADIX_ROOT); \
1123 carry = ((HD_HIGH(product)) + 1); \
1126 carry = (HD_HIGH(product)); \
1130 #define BDDS_ADD(vn, un, carry_in) \
1132 sum = (vn + un + carry_in); \
1133 if (sum < BIGNUM_RADIX_ROOT) { \
1137 un = (sum - BIGNUM_RADIX_ROOT); \
1142 bignum_digit_type factor_vm::bignum_digit_divide_subtract(
1143 bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
1144 bignum_digit_type* u) {
1146 bignum_digit_type product;
1147 bignum_digit_type diff;
1148 bignum_digit_type carry;
1149 BDDS_MULSUB(v2, (u[2]), 0);
1150 BDDS_MULSUB(v1, (u[1]), carry);
1153 diff = ((u[0]) - carry);
1155 (u[0]) = (diff + BIGNUM_RADIX);
1162 bignum_digit_type sum;
1163 bignum_digit_type carry;
1164 BDDS_ADD(v2, (u[2]), 0);
1165 BDDS_ADD(v1, (u[1]), carry);
1175 /* Allocates memory */
1176 void factor_vm::bignum_divide_unsigned_small_denominator(
1177 bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
1178 bignum** remainder, int q_negative_p, int r_negative_p) {
1179 data_root<bignum> numerator(numerator_, this);
1181 bignum* q_ = bignum_new_sign(numerator.untagged(), q_negative_p);
1182 data_root<bignum> q(q_, this);
1184 bignum_digit_type r = bignum_destructive_scale_down(q.untagged(), denominator);
1186 q = bignum_trim(q.untagged());
1188 if (remainder != ((bignum**)0))
1189 (*remainder) = bignum_digit_to_bignum(r, r_negative_p);
1191 (*quotient) = q.untagged();
1196 /* Given (denominator > 1), it is fairly easy to show that
1197 (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1198 that all digits are < BIGNUM_RADIX. */
1200 bignum_digit_type factor_vm::bignum_destructive_scale_down(
1201 bignum* bn, bignum_digit_type denominator) {
1202 bignum_digit_type numerator;
1203 bignum_digit_type remainder = 0;
1204 bignum_digit_type two_digits;
1205 #define quotient_high remainder
1206 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1207 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1208 BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1209 while (start < scan) {
1210 two_digits = (*--scan);
1211 numerator = (HD_CONS(remainder, (HD_HIGH(two_digits))));
1212 quotient_high = (numerator / denominator);
1213 numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits))));
1214 (*scan) = (HD_CONS(quotient_high, (numerator / denominator)));
1215 remainder = (numerator % denominator);
1218 #undef quotient_high
1221 /* Allocates memory */
1222 bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
1223 bignum* n, bignum_digit_type d, int negative_p) {
1224 bignum_digit_type two_digits;
1225 bignum_digit_type* start = (BIGNUM_START_PTR(n));
1226 bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
1227 bignum_digit_type r = 0;
1228 BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT));
1229 while (start < scan) {
1230 two_digits = (*--scan);
1231 r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d),
1232 (HD_LOW(two_digits)))) %
1235 return (bignum_digit_to_bignum(r, negative_p));
1238 /* Allocates memory */
1239 bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
1242 return (BIGNUM_ZERO());
1244 bignum* result = (allot_bignum(1, negative_p));
1245 (BIGNUM_REF(result, 0)) = digit;
1250 /* Allocates memory */
1251 bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
1252 BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
1253 bignum* result = allot_uninitialized_array<bignum>(length + 1);
1254 BIGNUM_SET_NEGATIVE_P(result, negative_p);
1258 /* Allocates memory */
1259 bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
1261 bignum* result = allot_bignum(length, negative_p);
1262 bignum_digit_type* scan = (BIGNUM_START_PTR(result));
1263 bignum_digit_type* end = (scan + length);
1269 /* Allocates memory */
1270 bignum* factor_vm::bignum_shorten_length(bignum* bn,
1271 bignum_length_type length) {
1272 bignum_length_type current_length = (BIGNUM_LENGTH(bn));
1273 BIGNUM_ASSERT((length >= 0) || (length <= current_length));
1274 if (length < current_length) {
1275 bn = reallot_array(bn, length + 1);
1276 BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1281 /* Allocates memory */
1282 bignum* factor_vm::bignum_trim(bignum* bn) {
1283 bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1284 bignum_digit_type* end = (start + (BIGNUM_LENGTH(bn)));
1285 bignum_digit_type* scan = end;
1286 while ((start <= scan) && ((*--scan) == 0))
1290 bignum_length_type length = (scan - start);
1291 bn = reallot_array(bn, length + 1);
1292 BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1299 /* Allocates memory */
1300 bignum* factor_vm::bignum_new_sign(bignum* x_, int negative_p) {
1301 data_root<bignum> x(x_, this);
1302 bignum* result = allot_bignum(BIGNUM_LENGTH(x), negative_p);
1303 bignum_destructive_copy(x.untagged(), result);
1307 /* Allocates memory */
1308 bignum* factor_vm::bignum_maybe_new_sign(bignum* x_, int negative_p) {
1309 if ((BIGNUM_NEGATIVE_P(x_)) ? negative_p : (!negative_p))
1312 return bignum_new_sign(x_, negative_p);
1316 void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
1317 bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1318 bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1319 bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1320 while (scan_source < end_source)
1321 (*scan_target++) = (*scan_source++);
1326 * Added bitwise operations (and oddp).
1329 /* Allocates memory */
1330 bignum* factor_vm::bignum_bitwise_not(bignum* x_) {
1333 bignum_length_type size = BIGNUM_LENGTH(x_);
1334 int is_negative = BIGNUM_NEGATIVE_P(x_);
1335 data_root<bignum> x(x_, this);
1336 data_root<bignum> y(allot_bignum(size, is_negative ? 0 : 1), this);
1338 bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
1339 bignum_digit_type* end_x = scan_x + size;
1340 bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
1343 while (scan_x < end_x) {
1345 *scan_y++ = BIGNUM_RADIX - 1;
1348 *scan_y++ = *scan_x++ - 1;
1354 while (scan_x < end_x) {
1355 if (*scan_x == (BIGNUM_RADIX - 1)) {
1359 *scan_y++ = *scan_x++ + 1;
1366 while (scan_x < end_x) {
1367 *scan_y++ = *scan_x++;
1371 bignum* ret = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
1372 bignum_destructive_copy(y.untagged(), ret);
1373 bignum_digit_type* ret_start = BIGNUM_START_PTR(ret);
1374 *(ret_start + size) = 1;
1377 return bignum_trim(y.untagged());
1381 /* Allocates memory */
1382 bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
1383 if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1384 return bignum_bitwise_not(
1385 bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1387 return bignum_magnitude_ash(arg1, n);
1394 /* Allocates memory */
1395 bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
1396 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1397 ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1398 : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1399 : (BIGNUM_NEGATIVE_P(arg2))
1400 ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1401 : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
1404 /* Allocates memory */
1405 bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
1406 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1407 ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1408 : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1409 : (BIGNUM_NEGATIVE_P(arg2))
1410 ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1411 : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
1414 /* Allocates memory */
1415 bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
1416 return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1417 ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1418 : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1419 : (BIGNUM_NEGATIVE_P(arg2))
1420 ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1421 : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
1424 /* Allocates memory */
1425 /* ash for the magnitude */
1426 /* assume arg1 is a big number, n is a long */
1427 bignum* factor_vm::bignum_magnitude_ash(bignum* arg1_, fixnum n) {
1429 data_root<bignum> arg1(arg1_, this);
1431 bignum* result = NULL;
1432 bignum_digit_type* scan1;
1433 bignum_digit_type* scanr;
1434 bignum_digit_type* end;
1436 fixnum digit_offset, bit_offset;
1438 if (BIGNUM_ZERO_P(arg1))
1439 return arg1.untagged();
1442 digit_offset = n / BIGNUM_DIGIT_LENGTH;
1443 bit_offset = n % BIGNUM_DIGIT_LENGTH;
1445 result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
1446 BIGNUM_NEGATIVE_P(arg1));
1448 scanr = BIGNUM_START_PTR(result) + digit_offset;
1449 scan1 = BIGNUM_START_PTR(arg1);
1450 end = scan1 + BIGNUM_LENGTH(arg1);
1452 while (scan1 < end) {
1453 *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1454 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1456 *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1457 *scanr = *scanr & BIGNUM_DIGIT_MASK;
1459 } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
1460 BIGNUM_DIGIT_LENGTH))) {
1461 result = BIGNUM_ZERO();
1463 digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1464 bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1466 result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
1467 BIGNUM_NEGATIVE_P(arg1));
1469 scanr = BIGNUM_START_PTR(result);
1470 scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
1471 end = scanr + BIGNUM_LENGTH(result) - 1;
1473 while (scanr < end) {
1474 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1475 *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
1479 *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1480 } else if (n == 0) {
1481 result = arg1.untagged();
1484 return bignum_trim(result);
1487 /* Allocates memory */
1488 bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1_,
1490 data_root<bignum> arg1(arg1_, this);
1491 data_root<bignum> arg2(arg2_, this);
1493 bignum_length_type max_length;
1495 bignum_digit_type* scan1, *end1, digit1;
1496 bignum_digit_type* scan2, *end2, digit2;
1497 bignum_digit_type* scanr, *endr;
1500 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
1501 : BIGNUM_LENGTH(arg2);
1503 bignum* result = allot_bignum(max_length, 0);
1505 scanr = BIGNUM_START_PTR(result);
1506 scan1 = BIGNUM_START_PTR(arg1);
1507 scan2 = BIGNUM_START_PTR(arg2);
1508 endr = scanr + max_length;
1509 end1 = scan1 + BIGNUM_LENGTH(arg1);
1510 end2 = scan2 + BIGNUM_LENGTH(arg2);
1512 while (scanr < endr) {
1513 digit1 = (scan1 < end1) ? *scan1++ : 0;
1514 digit2 = (scan2 < end2) ? *scan2++ : 0;
1516 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1519 return bignum_trim(result);
1522 /* Allocates memory */
1523 bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1_,
1525 data_root<bignum> arg1(arg1_, this);
1526 data_root<bignum> arg2(arg2_, this);
1528 bignum_length_type max_length;
1530 bignum_digit_type* scan1, *end1, digit1;
1531 bignum_digit_type* scan2, *end2, digit2, carry2;
1532 bignum_digit_type* scanr, *endr;
1534 char neg_p = op == IOR_OP || op == XOR_OP;
1537 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
1538 : BIGNUM_LENGTH(arg2) + 1;
1540 bignum* result = allot_bignum(max_length, neg_p);
1542 scanr = BIGNUM_START_PTR(result);
1543 scan1 = BIGNUM_START_PTR(arg1);
1544 scan2 = BIGNUM_START_PTR(arg2);
1545 endr = scanr + max_length;
1546 end1 = scan1 + BIGNUM_LENGTH(arg1);
1547 end2 = scan2 + BIGNUM_LENGTH(arg2);
1551 while (scanr < endr) {
1552 digit1 = (scan1 < end1) ? *scan1++ : 0;
1553 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1555 if (digit2 < BIGNUM_RADIX)
1558 digit2 = (digit2 - BIGNUM_RADIX);
1563 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1568 bignum_negate_magnitude(result);
1570 return bignum_trim(result);
1573 /* Allocates memory */
1574 bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1_,
1576 data_root<bignum> arg1(arg1_, this);
1577 data_root<bignum> arg2(arg2_, this);
1579 bignum_length_type max_length;
1581 bignum_digit_type* scan1, *end1, digit1, carry1;
1582 bignum_digit_type* scan2, *end2, digit2, carry2;
1583 bignum_digit_type* scanr, *endr;
1585 char neg_p = op == AND_OP || op == IOR_OP;
1588 (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
1589 : BIGNUM_LENGTH(arg2) + 1;
1591 bignum* result = allot_bignum(max_length, neg_p);
1593 scanr = BIGNUM_START_PTR(result);
1594 scan1 = BIGNUM_START_PTR(arg1);
1595 scan2 = BIGNUM_START_PTR(arg2);
1596 endr = scanr + max_length;
1597 end1 = scan1 + BIGNUM_LENGTH(arg1);
1598 end2 = scan2 + BIGNUM_LENGTH(arg2);
1603 while (scanr < endr) {
1604 digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1605 digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1607 if (digit1 < BIGNUM_RADIX)
1610 digit1 = (digit1 - BIGNUM_RADIX);
1614 if (digit2 < BIGNUM_RADIX)
1617 digit2 = (digit2 - BIGNUM_RADIX);
1622 (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1627 bignum_negate_magnitude(result);
1629 return bignum_trim(result);
1632 void factor_vm::bignum_negate_magnitude(bignum* arg) {
1633 bignum_digit_type* scan;
1634 bignum_digit_type* end;
1635 bignum_digit_type digit;
1636 bignum_digit_type carry;
1638 scan = BIGNUM_START_PTR(arg);
1639 end = scan + BIGNUM_LENGTH(arg);
1643 while (scan < end) {
1644 digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;
1646 if (digit < BIGNUM_RADIX)
1649 digit = (digit - BIGNUM_RADIX);
1657 /* Allocates memory */
1658 bignum* factor_vm::bignum_integer_length(bignum* x_) {
1659 data_root<bignum> x(x_, this);
1660 bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
1661 bignum_digit_type digit = (BIGNUM_REF(x, index));
1662 bignum_digit_type carry = 0;
1670 if (index < BIGNUM_RADIX_ROOT) {
1671 result = allot_bignum(1, 0);
1672 (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
1674 result = allot_bignum(2, 0);
1675 (BIGNUM_REF(result, 0)) = index;
1676 (BIGNUM_REF(result, 1)) = 0;
1677 bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
1678 bignum_destructive_add(result, carry);
1680 return (bignum_trim(result));
1683 /* Allocates memory */
1684 int factor_vm::bignum_logbitp(int shift, bignum* arg) {
1685 return ((BIGNUM_NEGATIVE_P(arg))
1686 ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
1687 : bignum_unsigned_logbitp(shift, arg));
1690 int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bn) {
1691 bignum_length_type len = (BIGNUM_LENGTH(bn));
1692 int index = shift / BIGNUM_DIGIT_LENGTH;
1695 bignum_digit_type digit = (BIGNUM_REF(bn, index));
1696 int p = shift % BIGNUM_DIGIT_LENGTH;
1697 bignum_digit_type mask = ((fixnum)1) << p;
1698 return (digit & mask) ? 1 : 0;
1702 /* Allocates memory. */
1703 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1705 data_root<bignum> a(a_, this);
1706 data_root<bignum> b(b_, this);
1708 /* Copies of a and b with that are both positive. */
1709 data_root<bignum> ac(bignum_maybe_new_sign(a.untagged(), 0), this);
1710 data_root<bignum> bc(bignum_maybe_new_sign(b.untagged(), 0), this);
1712 if (bignum_compare(ac.untagged(), bc.untagged()) == bignum_comparison_less) {
1716 while (BIGNUM_LENGTH(bc) != 0) {
1717 data_root<bignum> d(bignum_remainder(ac.untagged(), bc.untagged()), this);
1718 if (d.untagged() == BIGNUM_OUT_OF_BAND) {
1719 return d.untagged();
1724 return ac.untagged();
1727 /* Allocates memory */
1728 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1729 data_root<bignum> a(a_, this);
1730 data_root<bignum> b(b_, this);
1731 bignum_twodigit_type x, y, q, s, t, A, B, C, D;
1733 bignum_length_type size_a, size_b, size_c;
1734 bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
1735 bignum_digit_type* a_end, *b_end, *c_end;
1737 /* clone the bignums so we can modify them in-place */
1738 size_a = BIGNUM_LENGTH(a);
1739 data_root<bignum> c(allot_bignum(size_a, 0), this);
1740 // c = allot_bignum(size_a, 0);
1741 scan_a = BIGNUM_START_PTR(a);
1742 a_end = scan_a + size_a;
1743 scan_c = BIGNUM_START_PTR(c);
1744 while (scan_a < a_end)
1745 (*scan_c++) = (*scan_a++);
1747 size_b = BIGNUM_LENGTH(b);
1748 data_root<bignum> d(allot_bignum(size_b, 0), this);
1749 scan_b = BIGNUM_START_PTR(b);
1750 b_end = scan_b + size_b;
1751 scan_d = BIGNUM_START_PTR(d);
1752 while (scan_b < b_end)
1753 (*scan_d++) = (*scan_b++);
1756 /* Initial reduction: make sure that 0 <= b <= a. */
1757 if (bignum_compare(a.untagged(), b.untagged()) == bignum_comparison_less) {
1759 std::swap(size_a, size_b);
1762 while (size_a > 1) {
1763 nbits = log2(BIGNUM_REF(a, size_a - 1));
1764 x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
1765 (BIGNUM_REF(a, size_a - 2) >> nbits));
1766 y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
1768 ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
1771 /* inner loop of Lehmer's algorithm; */
1780 q = (x + (A - 1)) / (y - C);
1800 /* no progress; do a Euclidean step */
1802 return bignum_trim(a.untagged());
1804 data_root<bignum> e(bignum_trim(a.untagged()), this);
1805 data_root<bignum> f(bignum_trim(b.untagged()), this);
1806 c = bignum_remainder(e.untagged(), f.untagged());
1807 if (c.untagged() == BIGNUM_OUT_OF_BAND) {
1808 return c.untagged();
1812 scan_a = BIGNUM_START_PTR(a);
1813 scan_b = BIGNUM_START_PTR(b);
1814 a_end = scan_a + size_a;
1815 b_end = scan_b + size_b;
1816 while (scan_b < b_end)
1817 *(scan_a++) = *(scan_b++);
1818 while (scan_a < a_end)
1823 scan_b = BIGNUM_START_PTR(b);
1824 scan_c = BIGNUM_START_PTR(c);
1825 size_c = BIGNUM_LENGTH(c);
1826 c_end = scan_c + size_c;
1827 while (scan_c < c_end)
1828 *(scan_b++) = *(scan_c++);
1829 while (scan_b < b_end)
1837 a, b = A*b - B*a, D*a - C*b if k is odd
1838 a, b = A*a - B*b, D*b - C*a if k is even
1840 scan_a = BIGNUM_START_PTR(a);
1841 scan_b = BIGNUM_START_PTR(b);
1844 a_end = scan_a + size_a;
1845 b_end = scan_b + size_b;
1849 while (scan_b < b_end) {
1850 s += (A * *scan_b) - (B * *scan_a);
1851 t += (D * *scan_a++) - (C * *scan_b++);
1852 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1853 *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1854 s >>= BIGNUM_DIGIT_LENGTH;
1855 t >>= BIGNUM_DIGIT_LENGTH;
1857 while (scan_a < a_end) {
1859 t += (D * *scan_a++);
1860 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1861 //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1862 s >>= BIGNUM_DIGIT_LENGTH;
1863 t >>= BIGNUM_DIGIT_LENGTH;
1866 while (scan_b < b_end) {
1867 s += (A * *scan_a) - (B * *scan_b);
1868 t += (D * *scan_b++) - (C * *scan_a++);
1869 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1870 *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1871 s >>= BIGNUM_DIGIT_LENGTH;
1872 t >>= BIGNUM_DIGIT_LENGTH;
1874 while (scan_a < a_end) {
1876 t -= (C * *scan_a++);
1877 *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1878 //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1879 s >>= BIGNUM_DIGIT_LENGTH;
1880 t >>= BIGNUM_DIGIT_LENGTH;
1883 BIGNUM_ASSERT(s == 0);
1884 BIGNUM_ASSERT(t == 0);
1886 // update size_a and size_b to remove any zeros at end
1887 while (size_a > 0 && *(--scan_a) == 0)
1889 while (size_b > 0 && *(--scan_b) == 0)
1892 BIGNUM_ASSERT(size_a >= size_b);
1895 /* a fits into a fixnum, so b must too */
1896 fixnum xx = bignum_to_fixnum(a.untagged());
1897 fixnum yy = bignum_to_fixnum(b.untagged());
1900 /* usual Euclidean algorithm for longs */
1907 return fixnum_to_bignum(xx);