// Copyright (C) 1989-94 Massachusetts Institute of Technology // Portions copyright (C) 2004-2008 Slava Pestov // This material was developed by the Scheme project at the Massachusetts // Institute of Technology, Department of Electrical Engineering and // Computer Science. Permission to copy and modify this software, to // redistribute either the original software or a modified version, and // to use this software for any purpose is granted, subject to the // following restrictions and understandings. // 1. Any copy made of this software must include this copyright notice // in full. // 2. Users of this software agree to make their best efforts (a) to // return to the MIT Scheme project any improvements or extensions that // they make, so that these may be included in future releases; and (b) // to inform MIT of noteworthy uses of this software. // 3. All materials developed as a consequence of the use of this // software shall duly acknowledge such use, in accordance with the usual // standards of acknowledging credit in academic research. // 4. MIT has made no warrantee or representation that the operation of // this software will be error-free, and MIT is under no obligation to // provide any services, by way of maintenance, update, or otherwise. // 5. In conjunction with products arising from the use of this material, // there shall be no use of the name of the Massachusetts Institute of // Technology nor of any adaptation thereof in any advertising, // promotional, or sales literature without prior written consent from // MIT in each case. // Changes for Scheme 48: // * - Converted to ANSI. // * - Added bitwise operations. // * - Added s48 to the beginning of all externally visible names. // * - Cached the bignum representations of -1, 0, and 1. // Changes for Factor: // * - Adapt bignumint.h for Factor memory manager // * - Add more bignum <-> C type conversions // * - Remove unused functions // * - Add local variable GC root recording // * - Remove s48 prefix from function names // * - Various fixes for Win64 // * - Port to C++ // * - Added bignum_gcd implementation #include "master.hpp" namespace factor { // Exports int factor_vm::bignum_equal_p(bignum* x, bignum* y) { return ((BIGNUM_ZERO_P(x)) ? (BIGNUM_ZERO_P(y)) : ((!(BIGNUM_ZERO_P(y))) && ((BIGNUM_NEGATIVE_P(x)) ? (BIGNUM_NEGATIVE_P(y)) : (!(BIGNUM_NEGATIVE_P(y)))) && (bignum_equal_p_unsigned(x, y)))); } enum bignum_comparison factor_vm::bignum_compare(bignum* x, bignum* y) { return ((BIGNUM_ZERO_P(x)) ? ((BIGNUM_ZERO_P(y)) ? BIGNUM_COMPARISON_EQUAL : (BIGNUM_NEGATIVE_P(y)) ? BIGNUM_COMPARISON_GREATER : BIGNUM_COMPARISON_LESS) : (BIGNUM_ZERO_P(y)) ? ((BIGNUM_NEGATIVE_P(x)) ? BIGNUM_COMPARISON_LESS : BIGNUM_COMPARISON_GREATER) : (BIGNUM_NEGATIVE_P(x)) ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_compare_unsigned(y, x)) : (BIGNUM_COMPARISON_LESS)) : ((BIGNUM_NEGATIVE_P(y)) ? (BIGNUM_COMPARISON_GREATER) : (bignum_compare_unsigned(x, y)))); } // Allocates memory bignum* factor_vm::bignum_add(bignum* x, bignum* y) { return ( (BIGNUM_ZERO_P(x)) ? (y) : (BIGNUM_ZERO_P(y)) ? (x) : ((BIGNUM_NEGATIVE_P(x)) ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_add_unsigned(x, y, 1)) : (bignum_subtract_unsigned(y, x))) : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_subtract_unsigned(x, y)) : (bignum_add_unsigned(x, y, 0))))); } // Allocates memory bignum* factor_vm::bignum_subtract(bignum* x, bignum* y) { return ((BIGNUM_ZERO_P(x)) ? ((BIGNUM_ZERO_P(y)) ? (y) : (bignum_new_sign( y, (!(BIGNUM_NEGATIVE_P(y)))))) : ((BIGNUM_ZERO_P(y)) ? (x) : ((BIGNUM_NEGATIVE_P(x)) ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_subtract_unsigned(y, x)) : (bignum_add_unsigned(x, y, 1))) : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_add_unsigned(x, y, 0)) : (bignum_subtract_unsigned(x, y)))))); } #ifdef _WIN64 bignum *factor_vm::bignum_square(bignum* x_) { return bignum_multiply(x_, x_); } #else // Allocates memory bignum *factor_vm::bignum_square(bignum* x_) { data_root x(x_, this); bignum_length_type length = (BIGNUM_LENGTH (x)); bignum * z = (allot_bignum_zeroed ((length + length), 0)); bignum_digit_type * scan_z = BIGNUM_START_PTR (z); bignum_digit_type * scan_x = BIGNUM_START_PTR (x); bignum_digit_type * end_x = scan_x + length; for (int i = 0; i < length; ++i) { bignum_twodigit_type carry; bignum_twodigit_type f = BIGNUM_REF (x, i); bignum_digit_type *pz = scan_z + (i << 1); bignum_digit_type *px = scan_x + i + 1; carry = *pz + f * f; *pz++ = carry & BIGNUM_DIGIT_MASK; carry >>= BIGNUM_DIGIT_LENGTH; BIGNUM_ASSERT (carry <= BIGNUM_DIGIT_MASK); f <<= 1; while (px < end_x) { carry += *pz + *px++ * f; *pz++ = carry & BIGNUM_DIGIT_MASK; carry >>= BIGNUM_DIGIT_LENGTH; BIGNUM_ASSERT (carry <= (BIGNUM_DIGIT_MASK << 1)); } if (carry) { carry += *pz; *pz++ = carry & BIGNUM_DIGIT_MASK; carry >>= BIGNUM_DIGIT_LENGTH; } if (carry) *pz += carry & BIGNUM_DIGIT_MASK; BIGNUM_ASSERT ((carry >> BIGNUM_DIGIT_LENGTH) == 0); } return (bignum_trim (z)); } #endif // Allocates memory bignum* factor_vm::bignum_multiply(bignum* x, bignum* y) { #ifndef _WIN64 if (x == y) { return bignum_square(x); } #endif bignum_length_type x_length = (BIGNUM_LENGTH(x)); bignum_length_type y_length = (BIGNUM_LENGTH(y)); int negative_p = ((BIGNUM_NEGATIVE_P(x)) ? (!(BIGNUM_NEGATIVE_P(y))) : (BIGNUM_NEGATIVE_P(y))); if (BIGNUM_ZERO_P(x)) return (x); if (BIGNUM_ZERO_P(y)) return (y); if (x_length == 1) { bignum_digit_type digit = (BIGNUM_REF(x, 0)); if (digit == 1) return (bignum_maybe_new_sign(y, negative_p)); if (digit < BIGNUM_RADIX_ROOT) return (bignum_multiply_unsigned_small_factor(y, digit, negative_p)); } if (y_length == 1) { bignum_digit_type digit = (BIGNUM_REF(y, 0)); if (digit == 1) return (bignum_maybe_new_sign(x, negative_p)); if (digit < BIGNUM_RADIX_ROOT) return (bignum_multiply_unsigned_small_factor(x, digit, negative_p)); } return (bignum_multiply_unsigned(x, y, negative_p)); } // Allocates memory void factor_vm::bignum_divide(bignum* numerator, bignum* denominator, bignum** quotient, bignum** remainder) { if (BIGNUM_ZERO_P(denominator)) { divide_by_zero_error(); return; } if (BIGNUM_ZERO_P(numerator)) { (*quotient) = numerator; (*remainder) = numerator; } else { int r_negative_p = (BIGNUM_NEGATIVE_P(numerator)); int q_negative_p = ((BIGNUM_NEGATIVE_P(denominator)) ? (!r_negative_p) : r_negative_p); switch (bignum_compare_unsigned(numerator, denominator)) { case BIGNUM_COMPARISON_EQUAL: { (*quotient) = (BIGNUM_ONE(q_negative_p)); (*remainder) = (BIGNUM_ZERO()); break; } case BIGNUM_COMPARISON_LESS: { (*quotient) = (BIGNUM_ZERO()); (*remainder) = numerator; break; } case BIGNUM_COMPARISON_GREATER: { if ((BIGNUM_LENGTH(denominator)) == 1) { bignum_digit_type digit = (BIGNUM_REF(denominator, 0)); if (digit == 1) { (*quotient) = (bignum_maybe_new_sign(numerator, q_negative_p)); (*remainder) = (BIGNUM_ZERO()); break; } else if (digit < BIGNUM_RADIX_ROOT) { bignum_divide_unsigned_small_denominator(numerator, digit, quotient, remainder, q_negative_p, r_negative_p); break; } else { bignum_divide_unsigned_medium_denominator( numerator, digit, quotient, remainder, q_negative_p, r_negative_p); break; } } bignum_divide_unsigned_large_denominator( numerator, denominator, quotient, remainder, q_negative_p, r_negative_p); break; } } } } // Allocates memory bignum* factor_vm::bignum_quotient(bignum* numerator, bignum* denominator) { if (BIGNUM_ZERO_P(denominator)) { divide_by_zero_error(); return (BIGNUM_OUT_OF_BAND); } if (BIGNUM_ZERO_P(numerator)) return numerator; { int q_negative_p = ((BIGNUM_NEGATIVE_P(denominator)) ? (!(BIGNUM_NEGATIVE_P(numerator))) : (BIGNUM_NEGATIVE_P(numerator))); switch (bignum_compare_unsigned(numerator, denominator)) { case BIGNUM_COMPARISON_EQUAL: return (BIGNUM_ONE(q_negative_p)); case BIGNUM_COMPARISON_LESS: return (BIGNUM_ZERO()); case BIGNUM_COMPARISON_GREATER: default: // to appease gcc -Wall { bignum* quotient; if ((BIGNUM_LENGTH(denominator)) == 1) { bignum_digit_type digit = (BIGNUM_REF(denominator, 0)); if (digit == 1) return (bignum_maybe_new_sign(numerator, q_negative_p)); if (digit < BIGNUM_RADIX_ROOT) bignum_divide_unsigned_small_denominator( numerator, digit, ("ient), ((bignum**)0), q_negative_p, 0); else bignum_divide_unsigned_medium_denominator( numerator, digit, ("ient), ((bignum**)0), q_negative_p, 0); } else bignum_divide_unsigned_large_denominator( numerator, denominator, ("ient), ((bignum**)0), q_negative_p, 0); return (quotient); } } } } // Allocates memory bignum* factor_vm::bignum_remainder(bignum* numerator, bignum* denominator) { if (BIGNUM_ZERO_P(denominator)) { divide_by_zero_error(); return (BIGNUM_OUT_OF_BAND); } if (BIGNUM_ZERO_P(numerator)) return numerator; switch (bignum_compare_unsigned(numerator, denominator)) { case BIGNUM_COMPARISON_EQUAL: return (BIGNUM_ZERO()); case BIGNUM_COMPARISON_LESS: return numerator; case BIGNUM_COMPARISON_GREATER: default: // to appease gcc -Wall { bignum* remainder; if ((BIGNUM_LENGTH(denominator)) == 1) { bignum_digit_type digit = (BIGNUM_REF(denominator, 0)); if (digit == 1) return (BIGNUM_ZERO()); if (digit < BIGNUM_RADIX_ROOT) return (bignum_remainder_unsigned_small_denominator( numerator, digit, (BIGNUM_NEGATIVE_P(numerator)))); bignum_divide_unsigned_medium_denominator( numerator, digit, ((bignum**)0), (&remainder), 0, (BIGNUM_NEGATIVE_P(numerator))); } else bignum_divide_unsigned_large_denominator( numerator, denominator, ((bignum**)0), (&remainder), 0, (BIGNUM_NEGATIVE_P(numerator))); return (remainder); } } } // cell_to_bignum, fixnum_to_bignum, long_long_to_bignum, ulong_long_to_bignum // Allocates memory #define FOO_TO_BIGNUM(name, type, stype, utype) \ bignum* factor_vm::name##_to_bignum(type n) { \ int negative_p; \ /* Special cases win when these small constants are cached. */ \ if (n == 0) \ return (BIGNUM_ZERO()); \ if (n == 1) \ return (BIGNUM_ONE(0)); \ if (n < (type) 0 && n == (type) - 1) \ return (BIGNUM_ONE(1)); \ { \ utype accumulator = \ ((negative_p = (n < (type) 0)) ? ((type)(-(stype) n)) : n); \ if (accumulator < BIGNUM_RADIX) \ { \ bignum* result = allot_bignum(1, negative_p); \ bignum_digit_type* scan = (BIGNUM_START_PTR(result)); \ *scan = (accumulator & BIGNUM_DIGIT_MASK); \ return (result); \ } else { \ bignum_digit_type result_digits[BIGNUM_DIGITS_FOR(type)]; \ bignum_digit_type* end_digits = result_digits; \ do { \ (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); \ accumulator >>= BIGNUM_DIGIT_LENGTH; \ } while (accumulator != 0); \ bignum* result = \ (allot_bignum((end_digits - result_digits), negative_p)); \ bignum_digit_type* scan_digits = result_digits; \ bignum_digit_type* scan_result = (BIGNUM_START_PTR(result)); \ while (scan_digits < end_digits) \ (*scan_result++) = (*scan_digits++); \ return (result); \ } \ } \ } FOO_TO_BIGNUM(cell, cell, fixnum, cell) FOO_TO_BIGNUM(fixnum, fixnum, fixnum, cell) FOO_TO_BIGNUM(long_long, int64_t, int64_t, uint64_t) FOO_TO_BIGNUM(ulong_long, uint64_t, int64_t, uint64_t) // cannot allocate memory // bignum_to_cell, fixnum_to_cell, long_long_to_cell, ulong_long_to_cell #define BIGNUM_TO_FOO(name, type, stype, utype) \ type bignum_to_##name(bignum* bn) { \ if (BIGNUM_ZERO_P(bn)) \ return (0); \ { \ utype accumulator = 0; \ bignum_digit_type* start = (BIGNUM_START_PTR(bn)); \ bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn))); \ while (start < scan) \ accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); \ return ((BIGNUM_NEGATIVE_P(bn)) ? ((type)(-(stype) accumulator)) \ : accumulator); \ } \ } BIGNUM_TO_FOO(cell, cell, fixnum, cell) BIGNUM_TO_FOO(fixnum, fixnum, fixnum, cell) BIGNUM_TO_FOO(long_long, int64_t, int64_t, uint64_t) BIGNUM_TO_FOO(ulong_long, uint64_t, int64_t, uint64_t) bool bignum_fits_fixnum_p(bignum* bn) { fixnum len = BIGNUM_LENGTH(bn); if (len == 0) return true; if (len > 1) return false; bignum_digit_type dig = BIGNUM_START_PTR(bn)[0]; return (BIGNUM_NEGATIVE_P(bn) && dig <= -fixnum_min) || (!BIGNUM_NEGATIVE_P(bn) && dig <= fixnum_max); } cell bignum_maybe_to_fixnum(bignum* bn) { if (bignum_fits_fixnum_p(bn)) return tag_fixnum(bignum_to_fixnum(bn)); return tag(bn); } // cannot allocate memory fixnum factor_vm::bignum_to_fixnum_strict(bignum* bn) { if (!bignum_fits_fixnum_p(bn)) { general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag(bn), false_object); } fixnum fix = bignum_to_fixnum(bn); FACTOR_ASSERT(fix <= fixnum_max && fix >= fixnum_min); return fix; } #define DTB_WRITE_DIGIT(factor) \ { \ significand *= (factor); \ digit = ((bignum_digit_type) significand); \ (*--scan) = digit; \ significand -= ((double)digit); \ } #define inf std::numeric_limits::infinity() // Allocates memory bignum* factor_vm::double_to_bignum(double x) { if (x == inf || x == -inf || x != x) return (BIGNUM_ZERO()); int exponent; double significand = (frexp(x, (&exponent))); if (exponent <= 0) return (BIGNUM_ZERO()); if (exponent == 1) return (BIGNUM_ONE(x < 0)); if (significand < 0) significand = (-significand); { bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent)); bignum* result = (allot_bignum(length, (x < 0))); bignum_digit_type* start = (BIGNUM_START_PTR(result)); bignum_digit_type* scan = (start + length); bignum_digit_type digit; int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH); if (odd_bits > 0) DTB_WRITE_DIGIT((fixnum)1 << odd_bits); while (start < scan) { if (significand == 0) { while (start < scan) (*--scan) = 0; break; } DTB_WRITE_DIGIT(BIGNUM_RADIX); } return (result); } } #undef DTB_WRITE_DIGIT // Comparisons int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) { bignum_length_type length = (BIGNUM_LENGTH(x)); if (length != (BIGNUM_LENGTH(y))) return (0); else { bignum_digit_type* scan_x = (BIGNUM_START_PTR(x)); bignum_digit_type* scan_y = (BIGNUM_START_PTR(y)); bignum_digit_type* end_x = (scan_x + length); while (scan_x < end_x) if ((*scan_x++) != (*scan_y++)) return (0); return (1); } } enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x, bignum* y) { bignum_length_type x_length = (BIGNUM_LENGTH(x)); bignum_length_type y_length = (BIGNUM_LENGTH(y)); if (x_length < y_length) return BIGNUM_COMPARISON_LESS; if (x_length > y_length) return BIGNUM_COMPARISON_GREATER; { bignum_digit_type* start_x = (BIGNUM_START_PTR(x)); bignum_digit_type* scan_x = (start_x + x_length); bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length); while (start_x < scan_x) { bignum_digit_type digit_x = (*--scan_x); bignum_digit_type digit_y = (*--scan_y); if (digit_x < digit_y) return BIGNUM_COMPARISON_LESS; if (digit_x > digit_y) return BIGNUM_COMPARISON_GREATER; } } return BIGNUM_COMPARISON_EQUAL; } // Addition // Allocates memory bignum* factor_vm::bignum_add_unsigned(bignum* x_, bignum* y_, int negative_p) { data_root x(x_, this); data_root y(y_, this); if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) { swap(x, y); } { bignum_length_type x_length = (BIGNUM_LENGTH(x)); bignum* r = (allot_bignum((x_length + 1), negative_p)); bignum_digit_type sum; bignum_digit_type carry = 0; bignum_digit_type* scan_x = (BIGNUM_START_PTR(x)); bignum_digit_type* scan_r = (BIGNUM_START_PTR(r)); { bignum_digit_type* scan_y = (BIGNUM_START_PTR(y)); bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y))); while (scan_y < end_y) { sum = ((*scan_x++) + (*scan_y++) + carry); if (sum < BIGNUM_RADIX) { (*scan_r++) = sum; carry = 0; } else { (*scan_r++) = (sum - BIGNUM_RADIX); carry = 1; } } } { bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length; if (carry != 0) while (scan_x < end_x) { sum = ((*scan_x++) + 1); if (sum < BIGNUM_RADIX) { (*scan_r++) = sum; carry = 0; break; } else (*scan_r++) = (sum - BIGNUM_RADIX); } while (scan_x < end_x) (*scan_r++) = (*scan_x++); } if (carry != 0) { (*scan_r) = 1; return (r); } return (bignum_shorten_length(r, x_length)); } } // Subtraction // Allocates memory bignum* factor_vm::bignum_subtract_unsigned(bignum* x_, bignum* y_) { data_root x(x_, this); data_root y(y_, this); int negative_p = 0; switch (bignum_compare_unsigned(x.untagged(), y.untagged())) { case BIGNUM_COMPARISON_EQUAL: return (BIGNUM_ZERO()); case BIGNUM_COMPARISON_LESS: swap(x, y); negative_p = 1; break; case BIGNUM_COMPARISON_GREATER: negative_p = 0; break; } { bignum_length_type x_length = (BIGNUM_LENGTH(x)); bignum* r = (allot_bignum(x_length, negative_p)); bignum_digit_type difference; bignum_digit_type borrow = 0; bignum_digit_type* scan_x = BIGNUM_START_PTR(x); bignum_digit_type* scan_r = BIGNUM_START_PTR(r); { bignum_digit_type* scan_y = BIGNUM_START_PTR(y); bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y))); while (scan_y < end_y) { difference = (((*scan_x++) - (*scan_y++)) - borrow); if (difference < 0) { (*scan_r++) = (difference + BIGNUM_RADIX); borrow = 1; } else { (*scan_r++) = difference; borrow = 0; } } } { bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length; if (borrow != 0) while (scan_x < end_x) { difference = ((*scan_x++) - borrow); if (difference < 0) (*scan_r++) = (difference + BIGNUM_RADIX); else { (*scan_r++) = difference; borrow = 0; break; } } BIGNUM_ASSERT(borrow == 0); while (scan_x < end_x) (*scan_r++) = (*scan_x++); } return (bignum_trim(r)); } } // Multiplication // Maximum value for product_low or product_high: // ((R * R) + (R * (R - 2)) + (R - 1)) // Maximum value for carry: ((R * (R - 1)) + (R - 1)) // where R == BIGNUM_RADIX_ROOT // Allocates memory bignum* factor_vm::bignum_multiply_unsigned(bignum* x_, bignum* y_, int negative_p) { data_root x(x_, this); data_root y(y_, this); if (BIGNUM_LENGTH(y) > BIGNUM_LENGTH(x)) { swap(x, y); } { bignum_digit_type carry; bignum_digit_type y_digit_low; bignum_digit_type y_digit_high; bignum_digit_type x_digit_low; bignum_digit_type x_digit_high; bignum_digit_type product_low; bignum_digit_type* scan_r; bignum_digit_type* scan_y; bignum_length_type x_length = BIGNUM_LENGTH(x); bignum_length_type y_length = BIGNUM_LENGTH(y); bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p)); bignum_digit_type* scan_x = BIGNUM_START_PTR(x); bignum_digit_type* end_x = (scan_x + x_length); bignum_digit_type* start_y = BIGNUM_START_PTR(y); bignum_digit_type* end_y = (start_y + y_length); bignum_digit_type* start_r = (BIGNUM_START_PTR(r)); #define x_digit x_digit_high #define y_digit y_digit_high #define product_high carry while (scan_x < end_x) { x_digit = (*scan_x++); x_digit_low = (HD_LOW(x_digit)); x_digit_high = (HD_HIGH(x_digit)); carry = 0; scan_y = start_y; scan_r = (start_r++); while (scan_y < end_y) { y_digit = (*scan_y++); y_digit_low = (HD_LOW(y_digit)); y_digit_high = (HD_HIGH(y_digit)); product_low = ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry))); product_high = ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) + (HD_HIGH(product_low)) + (HD_HIGH(carry))); (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low)))); carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high))); } (*scan_r) += carry; } return (bignum_trim(r)); #undef x_digit #undef y_digit #undef product_high } } // Allocates memory bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x_, bignum_digit_type y, int negative_p) { data_root x(x_, this); bignum_length_type length_x = (BIGNUM_LENGTH(x)); bignum* p = (allot_bignum((length_x + 1), negative_p)); bignum_destructive_copy(x.untagged(), p); (BIGNUM_REF(p, length_x)) = 0; bignum_destructive_scale_up(p, y); return (bignum_trim(p)); } void factor_vm::bignum_destructive_add(bignum* bn, bignum_digit_type n) { bignum_digit_type* scan = (BIGNUM_START_PTR(bn)); bignum_digit_type digit; digit = ((*scan) + n); if (digit < BIGNUM_RADIX) { (*scan) = digit; return; } (*scan++) = (digit - BIGNUM_RADIX); while (1) { digit = ((*scan) + 1); if (digit < BIGNUM_RADIX) { (*scan) = digit; return; } (*scan++) = (digit - BIGNUM_RADIX); } } void factor_vm::bignum_destructive_scale_up(bignum* bn, bignum_digit_type factor) { bignum_digit_type carry = 0; bignum_digit_type* scan = (BIGNUM_START_PTR(bn)); bignum_digit_type two_digits; bignum_digit_type product_low; #define product_high carry bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bn))); BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT)); while (scan < end) { two_digits = (*scan); product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry))); product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) + (HD_HIGH(carry))); (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low)))); carry = (HD_HIGH(product_high)); } // A carry here would be an overflow, i.e. it would not fit. // Hopefully the callers allocate enough space that this will // never happen. BIGNUM_ASSERT(carry == 0); return; #undef product_high } // Division // For help understanding this algorithm, see: // Knuth, Donald E., "The Art of Computer Programming", // volume 2, "Seminumerical Algorithms" // section 4.3.1, "Multiple-Precision Arithmetic". // Allocates memory void factor_vm::bignum_divide_unsigned_large_denominator( bignum* numerator_, bignum* denominator_, bignum** quotient, bignum** remainder, int q_negative_p, int r_negative_p) { data_root numerator(numerator_, this); data_root denominator(denominator_, this); bignum_length_type length_n = BIGNUM_LENGTH(numerator) + 1; bignum_length_type length_d = BIGNUM_LENGTH(denominator); data_root u(allot_bignum(length_n, r_negative_p), this); int shift = 0; BIGNUM_ASSERT(length_d > 1); { bignum_digit_type v1 = BIGNUM_REF(denominator.untagged(), length_d - 1); while (v1 < (BIGNUM_RADIX / 2)) { v1 <<= 1; shift += 1; } } if (quotient != NULL) { bignum *q_ = allot_bignum(length_n - length_d, q_negative_p); data_root q(q_, this); if (shift == 0) { bignum_destructive_copy(numerator.untagged(), u.untagged()); (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0; bignum_divide_unsigned_normalized(u.untagged(), denominator.untagged(), q.untagged()); } else { bignum* v = allot_bignum(length_d, 0); bignum_destructive_normalization(numerator.untagged(), u.untagged(), shift); bignum_destructive_normalization(denominator.untagged(), v, shift); bignum_divide_unsigned_normalized(u.untagged(), v, q.untagged()); if (remainder != NULL) bignum_destructive_unnormalization(u.untagged(), shift); } q.set_untagged(bignum_trim(q.untagged())); *quotient = q.untagged(); } else { if (shift == 0) { bignum_destructive_copy(numerator.untagged(), u.untagged()); (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0; bignum_divide_unsigned_normalized(u.untagged(), denominator.untagged(), NULL); } else { bignum* v = allot_bignum(length_d, 0); bignum_destructive_normalization(numerator.untagged(), u.untagged(), shift); bignum_destructive_normalization(denominator.untagged(), v, shift); bignum_divide_unsigned_normalized(u.untagged(), v, NULL); if (remainder != NULL) bignum_destructive_unnormalization(u.untagged(), shift); } } u.set_untagged(bignum_trim(u.untagged())); if (remainder != NULL) *remainder = u.untagged(); } void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v, bignum* q) { bignum_length_type u_length = (BIGNUM_LENGTH(u)); bignum_length_type v_length = (BIGNUM_LENGTH(v)); bignum_digit_type* u_start = (BIGNUM_START_PTR(u)); bignum_digit_type* u_scan = (u_start + u_length); bignum_digit_type* u_scan_limit = (u_start + v_length); bignum_digit_type* u_scan_start = (u_scan - v_length); bignum_digit_type* v_start = (BIGNUM_START_PTR(v)); bignum_digit_type* v_end = (v_start + v_length); bignum_digit_type* q_scan = NULL; bignum_digit_type v1 = (v_end[-1]); bignum_digit_type v2 = (v_end[-2]); bignum_digit_type ph; // high half of double-digit product bignum_digit_type pl; // low half of double-digit product bignum_digit_type guess; bignum_digit_type gh; // high half-digit of guess bignum_digit_type ch; // high half of double-digit comparand bignum_digit_type v2l = (HD_LOW(v2)); bignum_digit_type v2h = (HD_HIGH(v2)); bignum_digit_type cl; // low half of double-digit comparand #define gl ph // low half-digit of guess #define uj pl #define qj ph bignum_digit_type gm; // memory loc for reference parameter if (q != BIGNUM_OUT_OF_BAND) q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q))); while (u_scan_limit < u_scan) { uj = (*--u_scan); if (uj != v1) { // comparand = // (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2); // guess = (((uj * BIGNUM_RADIX) + uj1) / v1); cl = (u_scan[-2]); ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm))); guess = gm; } else { cl = (u_scan[-2]); ch = ((u_scan[-1]) + v1); guess = (BIGNUM_RADIX - 1); } while (1) { // product = (guess * v2); gl = (HD_LOW(guess)); gh = (HD_HIGH(guess)); pl = (v2l * gl); ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl))); pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))); ph = ((v2h * gh) + (HD_HIGH(ph))); // if (comparand >= product) if ((ch > ph) || ((ch == ph) && (cl >= pl))) break; guess -= 1; // comparand += (v1 << BIGNUM_DIGIT_LENGTH) ch += v1; // if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) if (ch >= BIGNUM_RADIX) break; } qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start))); if (q != BIGNUM_OUT_OF_BAND) (*--q_scan) = qj; } return; #undef gl #undef uj #undef qj } bignum_digit_type factor_vm::bignum_divide_subtract( bignum_digit_type* v_start, bignum_digit_type* v_end, bignum_digit_type guess, bignum_digit_type* u_start) { bignum_digit_type* v_scan = v_start; bignum_digit_type* u_scan = u_start; bignum_digit_type carry = 0; if (guess == 0) return (0); { bignum_digit_type gl = (HD_LOW(guess)); bignum_digit_type gh = (HD_HIGH(guess)); bignum_digit_type v; bignum_digit_type pl; bignum_digit_type vl; #define vh v #define ph carry #define diff pl while (v_scan < v_end) { v = (*v_scan++); vl = (HD_LOW(v)); vh = (HD_HIGH(v)); pl = ((vl * gl) + (HD_LOW(carry))); ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry))); diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))))); if (diff < 0) { (*u_scan++) = (diff + BIGNUM_RADIX); carry = ((vh * gh) + (HD_HIGH(ph)) + 1); } else { (*u_scan++) = diff; carry = ((vh * gh) + (HD_HIGH(ph))); } } if (carry == 0) return (guess); diff = ((*u_scan) - carry); if (diff < 0) (*u_scan) = (diff + BIGNUM_RADIX); else { (*u_scan) = diff; return (guess); } #undef vh #undef ph #undef diff } // Subtraction generated carry, implying guess is one too large. // Add v back in to bring it back down. v_scan = v_start; u_scan = u_start; carry = 0; while (v_scan < v_end) { bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry); if (sum < BIGNUM_RADIX) { (*u_scan++) = sum; carry = 0; } else { (*u_scan++) = (sum - BIGNUM_RADIX); carry = 1; } } if (carry == 1) { bignum_digit_type sum = ((*u_scan) + carry); (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX)); } return (guess - 1); } // Allocates memory void factor_vm::bignum_divide_unsigned_medium_denominator( bignum* numerator_, bignum_digit_type denominator, bignum** quotient, bignum** remainder, int q_negative_p, int r_negative_p) { data_root numerator(numerator_, this); bignum_length_type length_n = (BIGNUM_LENGTH(numerator)); int shift = 0; // Because `bignum_digit_divide' requires a normalized denominator. while (denominator < (BIGNUM_RADIX / 2)) { denominator <<= 1; shift += 1; } bignum_length_type length_q = (shift == 0) ? length_n : length_n + 1; data_root q(allot_bignum(length_q, q_negative_p), this); if (shift == 0) { bignum_destructive_copy(numerator.untagged(), q.untagged()); } else { bignum_destructive_normalization(numerator.untagged(), q.untagged(), shift); } { bignum_digit_type r = 0; bignum_digit_type* start = (BIGNUM_START_PTR(q)); bignum_digit_type* scan = (start + length_q); bignum_digit_type qj; while (start < scan) { r = (bignum_digit_divide(r, (*--scan), denominator, (&qj))); (*scan) = qj; } q.set_untagged(bignum_trim(q.untagged())); if (remainder != ((bignum**)0)) { if (shift != 0) r >>= shift; (*remainder) = (bignum_digit_to_bignum(r, r_negative_p)); } if (quotient != ((bignum**)0)) (*quotient) = q.untagged(); } return; } void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target, int shift_left) { bignum_digit_type digit; bignum_digit_type* scan_source = (BIGNUM_START_PTR(source)); bignum_digit_type carry = 0; bignum_digit_type* scan_target = (BIGNUM_START_PTR(target)); bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source))); bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target))); int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left); bignum_digit_type mask = (((cell)1 << shift_right) - 1); while (scan_source < end_source) { digit = (*scan_source++); (*scan_target++) = (((digit & mask) << shift_left) | carry); carry = (digit >> shift_right); } if (scan_target < end_target) (*scan_target) = carry; else BIGNUM_ASSERT(carry == 0); return; } void factor_vm::bignum_destructive_unnormalization(bignum* bn, int shift_right) { bignum_digit_type* start = (BIGNUM_START_PTR(bn)); bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn))); bignum_digit_type digit; bignum_digit_type carry = 0; int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right); bignum_digit_type mask = (((fixnum)1 << shift_right) - 1); while (start < scan) { digit = (*--scan); (*scan) = ((digit >> shift_right) | carry); carry = ((digit & mask) << shift_left); } BIGNUM_ASSERT(carry == 0); return; } // This is a reduced version of the division algorithm, applied to the // case of dividing two bignum digits by one bignum digit. It is // assumed that the numerator, denominator are normalized. #define BDD_STEP(qn, j) \ { \ uj = (u[j]); \ if (uj != v1) { \ uj_uj1 = (HD_CONS(uj, (u[j + 1]))); \ guess = (uj_uj1 / v1); \ comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2]))); \ } else { \ guess = (BIGNUM_RADIX_ROOT - 1); \ comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2]))); \ } \ while ((guess * v2) > comparand) { \ guess -= 1; \ comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH); \ if (comparand >= BIGNUM_RADIX) \ break; \ } \ qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \ } bignum_digit_type factor_vm::bignum_digit_divide( bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v, bignum_digit_type* q) // return value { bignum_digit_type guess; bignum_digit_type comparand; bignum_digit_type v1 = (HD_HIGH(v)); bignum_digit_type v2 = (HD_LOW(v)); bignum_digit_type uj; bignum_digit_type uj_uj1; bignum_digit_type q1; bignum_digit_type q2; bignum_digit_type u[4]; if (uh == 0) { if (ul < v) { (*q) = 0; return (ul); } else if (ul == v) { (*q) = 1; return (0); } } (u[0]) = (HD_HIGH(uh)); (u[1]) = (HD_LOW(uh)); (u[2]) = (HD_HIGH(ul)); (u[3]) = (HD_LOW(ul)); v1 = (HD_HIGH(v)); v2 = (HD_LOW(v)); BDD_STEP(q1, 0); BDD_STEP(q2, 1); (*q) = (HD_CONS(q1, q2)); return (HD_CONS((u[2]), (u[3]))); } #undef BDD_STEP #define BDDS_MULSUB(vn, un, carry_in) \ { \ product = ((vn * guess) + carry_in); \ diff = (un - (HD_LOW(product))); \ if (diff < 0) { \ un = (diff + BIGNUM_RADIX_ROOT); \ carry = ((HD_HIGH(product)) + 1); \ } else { \ un = diff; \ carry = (HD_HIGH(product)); \ } \ } #define BDDS_ADD(vn, un, carry_in) \ { \ sum = (vn + un + carry_in); \ if (sum < BIGNUM_RADIX_ROOT) { \ un = sum; \ carry = 0; \ } else { \ un = (sum - BIGNUM_RADIX_ROOT); \ carry = 1; \ } \ } bignum_digit_type factor_vm::bignum_digit_divide_subtract( bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess, bignum_digit_type* u) { { bignum_digit_type product; bignum_digit_type diff; bignum_digit_type carry; BDDS_MULSUB(v2, (u[2]), 0); BDDS_MULSUB(v1, (u[1]), carry); if (carry == 0) return (guess); diff = ((u[0]) - carry); if (diff < 0) (u[0]) = (diff + BIGNUM_RADIX); else { (u[0]) = diff; return (guess); } } { bignum_digit_type sum; bignum_digit_type carry; BDDS_ADD(v2, (u[2]), 0); BDDS_ADD(v1, (u[1]), carry); if (carry == 1) (u[0]) += 1; } return (guess - 1); } #undef BDDS_MULSUB #undef BDDS_ADD // Allocates memory void factor_vm::bignum_divide_unsigned_small_denominator( bignum* numerator_, bignum_digit_type denominator, bignum** quotient, bignum** remainder, int q_negative_p, int r_negative_p) { data_root numerator(numerator_, this); bignum* q_ = bignum_new_sign(numerator.untagged(), q_negative_p); data_root q(q_, this); bignum_digit_type r = bignum_destructive_scale_down(q.untagged(), denominator); q.set_untagged(bignum_trim(q.untagged())); if (remainder != ((bignum**)0)) (*remainder) = bignum_digit_to_bignum(r, r_negative_p); (*quotient) = q.untagged(); return; } // Given (denominator > 1), it is fairly easy to show that // (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see // that all digits are < BIGNUM_RADIX. bignum_digit_type factor_vm::bignum_destructive_scale_down( bignum* bn, bignum_digit_type denominator) { bignum_digit_type numerator; bignum_digit_type remainder = 0; bignum_digit_type two_digits; #define quotient_high remainder bignum_digit_type* start = (BIGNUM_START_PTR(bn)); bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn))); BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT)); while (start < scan) { two_digits = (*--scan); numerator = (HD_CONS(remainder, (HD_HIGH(two_digits)))); quotient_high = (numerator / denominator); numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits)))); (*scan) = (HD_CONS(quotient_high, (numerator / denominator))); remainder = (numerator % denominator); } return (remainder); #undef quotient_high } // Allocates memory bignum* factor_vm::bignum_remainder_unsigned_small_denominator( bignum* n, bignum_digit_type d, int negative_p) { bignum_digit_type two_digits; bignum_digit_type* start = (BIGNUM_START_PTR(n)); bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n))); bignum_digit_type r = 0; BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT)); while (start < scan) { two_digits = (*--scan); r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d), (HD_LOW(two_digits)))) % d); } return (bignum_digit_to_bignum(r, negative_p)); } // Allocates memory bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit, int negative_p) { if (digit == 0) return (BIGNUM_ZERO()); else { bignum* result = (allot_bignum(1, negative_p)); (BIGNUM_REF(result, 0)) = digit; return (result); } } // Allocates memory bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) { BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX)); bignum* result = allot_uninitialized_array(length + 1); BIGNUM_SET_NEGATIVE_P(result, negative_p); return (result); } // Allocates memory bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length, int negative_p) { bignum* result = allot_bignum(length, negative_p); bignum_digit_type* scan = (BIGNUM_START_PTR(result)); bignum_digit_type* end = (scan + length); while (scan < end) (*scan++) = 0; return (result); } // Allocates memory bignum* factor_vm::bignum_shorten_length(bignum* bn, bignum_length_type length) { bignum_length_type current_length = (BIGNUM_LENGTH(bn)); BIGNUM_ASSERT((length >= 0) || (length <= current_length)); if (length < current_length) { bn = reallot_array(bn, length + 1); BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn))); } return (bn); } // Allocates memory bignum* factor_vm::bignum_trim(bignum* bn) { bignum_digit_type* start = (BIGNUM_START_PTR(bn)); bignum_digit_type* end = (start + (BIGNUM_LENGTH(bn))); bignum_digit_type* scan = end; while ((start <= scan) && ((*--scan) == 0)) ; scan += 1; if (scan < end) { bignum_length_type length = (scan - start); bn = reallot_array(bn, length + 1); BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn))); } return (bn); } // Copying // Allocates memory bignum* factor_vm::bignum_new_sign(bignum* x_, int negative_p) { data_root x(x_, this); bignum* result = allot_bignum(BIGNUM_LENGTH(x), negative_p); bignum_destructive_copy(x.untagged(), result); return result; } // Allocates memory bignum* factor_vm::bignum_maybe_new_sign(bignum* x_, int negative_p) { if ((BIGNUM_NEGATIVE_P(x_)) ? negative_p : (!negative_p)) return x_; else { return bignum_new_sign(x_, negative_p); } } void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) { bignum_digit_type* scan_source = (BIGNUM_START_PTR(source)); bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source))); bignum_digit_type* scan_target = (BIGNUM_START_PTR(target)); while (scan_source < end_source) (*scan_target++) = (*scan_source++); return; } // * Added bitwise operations (and oddp). // Allocates memory bignum* factor_vm::bignum_bitwise_not(bignum* x_) { int carry = 1; bignum_length_type size = BIGNUM_LENGTH(x_); int is_negative = BIGNUM_NEGATIVE_P(x_); data_root x(x_, this); data_root y(allot_bignum(size, is_negative ? 0 : 1), this); bignum_digit_type* scan_x = BIGNUM_START_PTR(x); bignum_digit_type* end_x = scan_x + size; bignum_digit_type* scan_y = BIGNUM_START_PTR(y); if (is_negative) { while (scan_x < end_x) { if (*scan_x == 0) { *scan_y++ = BIGNUM_RADIX - 1; scan_x++; } else { *scan_y++ = *scan_x++ - 1; carry = 0; break; } } } else { while (scan_x < end_x) { if (*scan_x == (BIGNUM_RADIX - 1)) { *scan_y++ = 0; scan_x++; } else { *scan_y++ = *scan_x++ + 1; carry = 0; break; } } } while (scan_x < end_x) { *scan_y++ = *scan_x++; } if (carry) { bignum* ret = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y)); bignum_destructive_copy(y.untagged(), ret); bignum_digit_type* ret_start = BIGNUM_START_PTR(ret); *(ret_start + size) = 1; return ret; } else { return bignum_trim(y.untagged()); } } // Allocates memory bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) { if (BIGNUM_NEGATIVE_P(arg1) && n < 0) return bignum_bitwise_not( bignum_magnitude_ash(bignum_bitwise_not(arg1), n)); else return bignum_magnitude_ash(arg1, n); } #define AND_OP 0 #define IOR_OP 1 #define XOR_OP 2 // Allocates memory bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) { return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2)) ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2) : bignum_posneg_bitwise_op(AND_OP, arg2, arg1) : (BIGNUM_NEGATIVE_P(arg2)) ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2) : bignum_pospos_bitwise_op(AND_OP, arg1, arg2)); } // Allocates memory bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) { return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2)) ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2) : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1) : (BIGNUM_NEGATIVE_P(arg2)) ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2) : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2)); } // Allocates memory bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) { return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2)) ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2) : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1) : (BIGNUM_NEGATIVE_P(arg2)) ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2) : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2)); } // Allocates memory // ash for the magnitude // assume arg1 is a big number, n is a long bignum* factor_vm::bignum_magnitude_ash(bignum* arg1_, fixnum n) { data_root arg1(arg1_, this); bignum* result = NULL; bignum_digit_type* scan1; bignum_digit_type* scanr; bignum_digit_type* end; fixnum digit_offset, bit_offset; if (BIGNUM_ZERO_P(arg1)) return arg1.untagged(); if (n > 0) { digit_offset = n / BIGNUM_DIGIT_LENGTH; bit_offset = n % BIGNUM_DIGIT_LENGTH; result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1, BIGNUM_NEGATIVE_P(arg1)); scanr = BIGNUM_START_PTR(result) + digit_offset; scan1 = BIGNUM_START_PTR(arg1); end = scan1 + BIGNUM_LENGTH(arg1); while (scan1 < end) { *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset; *scanr = *scanr & BIGNUM_DIGIT_MASK; scanr++; *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset); *scanr = *scanr & BIGNUM_DIGIT_MASK; } } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type) BIGNUM_DIGIT_LENGTH))) { result = BIGNUM_ZERO(); } else if (n < 0) { digit_offset = -n / BIGNUM_DIGIT_LENGTH; bit_offset = -n % BIGNUM_DIGIT_LENGTH; result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset, BIGNUM_NEGATIVE_P(arg1)); scanr = BIGNUM_START_PTR(result); scan1 = BIGNUM_START_PTR(arg1) + digit_offset; end = scanr + BIGNUM_LENGTH(result) - 1; while (scanr < end) { *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset; *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) & BIGNUM_DIGIT_MASK; scanr++; } *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset; } else if (n == 0) { result = arg1.untagged(); } return bignum_trim(result); } // Allocates memory bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1_, bignum* arg2_) { data_root arg1(arg1_, this); data_root arg2(arg2_, this); bignum_length_type max_length; bignum_digit_type* scan1, *end1, digit1; bignum_digit_type* scan2, *end2, digit2; bignum_digit_type* scanr, *endr; max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2); bignum* result = allot_bignum(max_length, 0); scanr = BIGNUM_START_PTR(result); scan1 = BIGNUM_START_PTR(arg1); scan2 = BIGNUM_START_PTR(arg2); endr = scanr + max_length; end1 = scan1 + BIGNUM_LENGTH(arg1); end2 = scan2 + BIGNUM_LENGTH(arg2); while (scanr < endr) { digit1 = (scan1 < end1) ? *scan1++ : 0; digit2 = (scan2 < end2) ? *scan2++ : 0; *scanr++ = (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2 : digit1 ^ digit2; } return bignum_trim(result); } // Allocates memory bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1_, bignum* arg2_) { data_root arg1(arg1_, this); data_root arg2(arg2_, this); bignum_length_type max_length; bignum_digit_type* scan1, *end1, digit1; bignum_digit_type* scan2, *end2, digit2, carry2; bignum_digit_type* scanr, *endr; char neg_p = op == IOR_OP || op == XOR_OP; max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1) : BIGNUM_LENGTH(arg2) + 1; bignum* result = allot_bignum(max_length, neg_p); scanr = BIGNUM_START_PTR(result); scan1 = BIGNUM_START_PTR(arg1); scan2 = BIGNUM_START_PTR(arg2); endr = scanr + max_length; end1 = scan1 + BIGNUM_LENGTH(arg1); end2 = scan2 + BIGNUM_LENGTH(arg2); carry2 = 1; while (scanr < endr) { digit1 = (scan1 < end1) ? *scan1++ : 0; digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2; if (digit2 < BIGNUM_RADIX) carry2 = 0; else { digit2 = (digit2 - BIGNUM_RADIX); carry2 = 1; } *scanr++ = (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2 : digit1 ^ digit2; } if (neg_p) bignum_negate_magnitude(result); return bignum_trim(result); } // Allocates memory bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1_, bignum* arg2_) { data_root arg1(arg1_, this); data_root arg2(arg2_, this); bignum_length_type max_length; bignum_digit_type* scan1, *end1, digit1, carry1; bignum_digit_type* scan2, *end2, digit2, carry2; bignum_digit_type* scanr, *endr; char neg_p = op == AND_OP || op == IOR_OP; max_length = (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1 : BIGNUM_LENGTH(arg2) + 1; bignum* result = allot_bignum(max_length, neg_p); scanr = BIGNUM_START_PTR(result); scan1 = BIGNUM_START_PTR(arg1); scan2 = BIGNUM_START_PTR(arg2); endr = scanr + max_length; end1 = scan1 + BIGNUM_LENGTH(arg1); end2 = scan2 + BIGNUM_LENGTH(arg2); carry1 = 1; carry2 = 1; while (scanr < endr) { digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1; digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2; if (digit1 < BIGNUM_RADIX) carry1 = 0; else { digit1 = (digit1 - BIGNUM_RADIX); carry1 = 1; } if (digit2 < BIGNUM_RADIX) carry2 = 0; else { digit2 = (digit2 - BIGNUM_RADIX); carry2 = 1; } *scanr++ = (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2 : digit1 ^ digit2; } if (neg_p) bignum_negate_magnitude(result); return bignum_trim(result); } void factor_vm::bignum_negate_magnitude(bignum* arg) { bignum_digit_type* scan; bignum_digit_type* end; bignum_digit_type digit; bignum_digit_type carry; scan = BIGNUM_START_PTR(arg); end = scan + BIGNUM_LENGTH(arg); carry = 1; while (scan < end) { digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry; if (digit < BIGNUM_RADIX) carry = 0; else { digit = (digit - BIGNUM_RADIX); carry = 1; } *scan++ = digit; } } // Allocates memory bignum* factor_vm::bignum_integer_length(bignum* x_) { data_root x(x_, this); bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1); bignum_digit_type digit = (BIGNUM_REF(x, index)); bignum_digit_type carry = 0; bignum* result; while (digit > 1) { carry += 1; digit >>= 1; } if (index < BIGNUM_RADIX_ROOT) { result = allot_bignum(1, 0); (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry; } else { result = allot_bignum(2, 0); (BIGNUM_REF(result, 0)) = index; (BIGNUM_REF(result, 1)) = 0; bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH); bignum_destructive_add(result, carry); } return (bignum_trim(result)); } // Allocates memory int factor_vm::bignum_logbitp(int shift, bignum* arg) { return ((BIGNUM_NEGATIVE_P(arg)) ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg)) : bignum_unsigned_logbitp(shift, arg)); } int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bn) { bignum_length_type len = (BIGNUM_LENGTH(bn)); int index = shift / BIGNUM_DIGIT_LENGTH; if (index >= len) return 0; bignum_digit_type digit = (BIGNUM_REF(bn, index)); int p = shift % BIGNUM_DIGIT_LENGTH; bignum_digit_type mask = ((fixnum)1) << p; return (digit & mask) ? 1 : 0; } #ifdef _WIN64 // Allocates memory. bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) { data_root a(a_, this); data_root b(b_, this); // Copies of a and b with that are both positive. data_root ac(bignum_maybe_new_sign(a.untagged(), 0), this); data_root bc(bignum_maybe_new_sign(b.untagged(), 0), this); if (bignum_compare(ac.untagged(), bc.untagged()) == BIGNUM_COMPARISON_LESS) { swap(ac, bc); } while (BIGNUM_LENGTH(bc) != 0) { data_root d(bignum_remainder(ac.untagged(), bc.untagged()), this); if (d.untagged() == BIGNUM_OUT_OF_BAND) { return d.untagged(); } ac = bc; bc = d; } return ac.untagged(); } #else // Allocates memory bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) { data_root a(a_, this); data_root b(b_, this); bignum_twodigit_type x, y, q, s, t, A, B, C, D; int nbits, k; bignum_length_type size_a, size_b, size_c; bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d; bignum_digit_type* a_end, *b_end, *c_end; // clone the bignums so we can modify them in-place size_a = BIGNUM_LENGTH(a); data_root c(allot_bignum(size_a, 0), this); // c = allot_bignum(size_a, 0); scan_a = BIGNUM_START_PTR(a); a_end = scan_a + size_a; scan_c = BIGNUM_START_PTR(c); while (scan_a < a_end) (*scan_c++) = (*scan_a++); a = c; size_b = BIGNUM_LENGTH(b); data_root d(allot_bignum(size_b, 0), this); scan_b = BIGNUM_START_PTR(b); b_end = scan_b + size_b; scan_d = BIGNUM_START_PTR(d); while (scan_b < b_end) (*scan_d++) = (*scan_b++); b = d; // Initial reduction: make sure that 0 <= b <= a. if (bignum_compare(a.untagged(), b.untagged()) == BIGNUM_COMPARISON_LESS) { swap(a, b); std::swap(size_a, size_b); } while (size_a > 1) { nbits = log2(BIGNUM_REF(a, size_a - 1)); x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) | (BIGNUM_REF(a, size_a - 2) >> nbits)); y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) | (size_b >= size_a ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits) : 0)); // inner loop of Lehmer's algorithm; A = 1; B = 0; C = 0; D = 1; for (k = 0;; k++) { if (y - C == 0) break; q = (x + (A - 1)) / (y - C); s = B + (q * D); t = x - (q * y); if (s > t) break; x = y; y = t; t = A + (q * C); A = D; B = C; C = s; D = t; } if (k == 0) { // no progress; do a Euclidean step if (size_b == 0) { return bignum_trim(a.untagged()); } data_root e(bignum_trim(a.untagged()), this); data_root f(bignum_trim(b.untagged()), this); c.set_untagged(bignum_remainder(e.untagged(), f.untagged())); if (c.untagged() == BIGNUM_OUT_OF_BAND) { return c.untagged(); } // copy 'b' to 'a' scan_a = BIGNUM_START_PTR(a); scan_b = BIGNUM_START_PTR(b); a_end = scan_a + size_a; b_end = scan_b + size_b; while (scan_b < b_end) *(scan_a++) = *(scan_b++); while (scan_a < a_end) *(scan_a++) = 0; size_a = size_b; // copy 'c' to 'b' scan_b = BIGNUM_START_PTR(b); scan_c = BIGNUM_START_PTR(c); size_c = BIGNUM_LENGTH(c); c_end = scan_c + size_c; while (scan_c < c_end) *(scan_b++) = *(scan_c++); while (scan_b < b_end) *(scan_b++) = 0; size_b = size_c; continue; } // a, b = A*b - B*a, D*a - C*b if k is odd // a, b = A*a - B*b, D*b - C*a if k is even scan_a = BIGNUM_START_PTR(a); scan_b = BIGNUM_START_PTR(b); scan_c = scan_a; scan_d = scan_b; a_end = scan_a + size_a; b_end = scan_b + size_b; s = 0; t = 0; if (k & 1) { while (scan_b < b_end) { s += (A * *scan_b) - (B * *scan_a); t += (D * *scan_a++) - (C * *scan_b++); *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK); *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK); s >>= BIGNUM_DIGIT_LENGTH; t >>= BIGNUM_DIGIT_LENGTH; } while (scan_a < a_end) { s -= (B * *scan_a); t += (D * *scan_a++); *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK); //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); s >>= BIGNUM_DIGIT_LENGTH; t >>= BIGNUM_DIGIT_LENGTH; } } else { while (scan_b < b_end) { s += (A * *scan_a) - (B * *scan_b); t += (D * *scan_b++) - (C * *scan_a++); *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK); *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK); s >>= BIGNUM_DIGIT_LENGTH; t >>= BIGNUM_DIGIT_LENGTH; } while (scan_a < a_end) { s += (A * *scan_a); t -= (C * *scan_a++); *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK); //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); s >>= BIGNUM_DIGIT_LENGTH; t >>= BIGNUM_DIGIT_LENGTH; } } BIGNUM_ASSERT(s == 0); BIGNUM_ASSERT(t == 0); // update size_a and size_b to remove any zeros at end while (size_a > 0 && *(--scan_a) == 0) size_a--; while (size_b > 0 && *(--scan_b) == 0) size_b--; BIGNUM_ASSERT(size_a >= size_b); } // a fits into a fixnum, so b must too fixnum xx = bignum_to_fixnum(a.untagged()); fixnum yy = bignum_to_fixnum(b.untagged()); fixnum tt; // usual Euclidean algorithm for longs while (yy != 0) { tt = yy; yy = xx % yy; xx = tt; } return fixnum_to_bignum(xx); } #endif }