From 63f1365820c318d8631f99032b3d4e79606cc19d Mon Sep 17 00:00:00 2001 From: Slava Pestov Date: Wed, 25 Aug 2004 03:46:55 +0000 Subject: [PATCH] Porting Scheme48 bignums to Factor. --- native/array.c | 16 +- native/array.h | 3 +- native/bignum.h | 10 +- native/complex.c | 2 +- native/factor.h | 2 + native/misc.c | 10 + native/misc.h | 1 + native/primitives.c | 3 +- native/primitives.h | 2 +- native/read.c | 6 +- native/relocate.c | 1 + native/s48_bignum.c | 2011 ++++++++++++++++++++++++++++++++++++++++ native/s48_bignum.h | 94 ++ native/s48_bignumint.h | 127 +++ native/types.c | 4 +- 15 files changed, 2272 insertions(+), 20 deletions(-) create mode 100644 native/s48_bignum.c create mode 100644 native/s48_bignum.h create mode 100644 native/s48_bignumint.h diff --git a/native/array.c b/native/array.c index 66ddf20d1b..52491d656f 100644 --- a/native/array.c +++ b/native/array.c @@ -1,10 +1,9 @@ #include "factor.h" /* untagged */ -ARRAY* allot_array(CELL capacity) +ARRAY* allot_array(CELL type, CELL capacity) { - ARRAY* array = allot_object(ARRAY_TYPE, - sizeof(ARRAY) + capacity * CELLS); + ARRAY* array = allot_object(type,sizeof(ARRAY) + capacity * CELLS); array->capacity = capacity; return array; } @@ -14,7 +13,7 @@ ARRAY* array(CELL capacity, CELL fill) { int i; - ARRAY* array = allot_array(capacity); + ARRAY* array = allot_array(ARRAY_TYPE, capacity); for(i = 0; i < capacity; i++) put(AREF(array,i),fill); @@ -27,7 +26,7 @@ ARRAY* grow_array(ARRAY* array, CELL capacity, CELL fill) /* later on, do an optimization: if end of array is here, just grow */ int i; - ARRAY* new_array = allot_array(capacity); + ARRAY* new_array = allot_array(untag_header(array->header),capacity); memcpy(new_array + 1,array + 1,array->capacity * CELLS); @@ -37,6 +36,13 @@ ARRAY* grow_array(ARRAY* array, CELL capacity, CELL fill) return new_array; } +ARRAY* shrink_array(ARRAY* array, CELL capacity) +{ + ARRAY* new_array = allot_array(untag_header(array->header),capacity); + memcpy(new_array + 1,array + 1,capacity * CELLS); + return new_array; +} + void fixup_array(ARRAY* array) { int i = 0; diff --git a/native/array.h b/native/array.h index fa904aa451..0ad56f06f5 100644 --- a/native/array.h +++ b/native/array.h @@ -10,9 +10,10 @@ INLINE ARRAY* untag_array(CELL tagged) return (ARRAY*)UNTAG(tagged); } +ARRAY* allot_array(CELL type, CELL capacity); ARRAY* array(CELL capacity, CELL fill); - ARRAY* grow_array(ARRAY* array, CELL capacity, CELL fill); +ARRAY* shrink_array(ARRAY* array, CELL capacity); #define AREF(array,index) ((CELL)array + sizeof(ARRAY) + index * CELLS) diff --git a/native/bignum.h b/native/bignum.h index fb2cd7bc9f..b8208f325a 100644 --- a/native/bignum.h +++ b/native/bignum.h @@ -2,17 +2,17 @@ typedef long long BIGNUM_2; typedef struct { CELL header; -/* FIXME */ -#ifndef FACTOR_64 - CELL alignment; -#endif + CELL capacity; + CELL sign; + CELL fill; /* bad */ BIGNUM_2 n; } BIGNUM; /* untagged */ INLINE BIGNUM* allot_bignum() { - return allot_object(BIGNUM_TYPE,sizeof(BIGNUM)); + /* Bignums are really retrofitted arrays */ + return (BIGNUM*)allot_array(BIGNUM_TYPE,4); } /* untagged */ diff --git a/native/complex.c b/native/complex.c index b838e136d2..0949e8c9fc 100644 --- a/native/complex.c +++ b/native/complex.c @@ -163,7 +163,7 @@ CELL divfloat_complex(CELL x, CELL y) } #define INCOMPARABLE(x,y) general_error(ERROR_INCOMPARABLE, \ - tag_cons(cons(tag_complex(x),tag_complex(y)))); + tag_cons(cons(RETAG(x,COMPLEX_TYPE),RETAG(y,COMPLEX_TYPE)))); CELL less_complex(CELL x, CELL y) { diff --git a/native/factor.h b/native/factor.h index 3596e125b8..7d9099021f 100644 --- a/native/factor.h +++ b/native/factor.h @@ -48,6 +48,8 @@ and allows profiling. */ #include "word.h" #include "run.h" #include "fixnum.h" +#include "s48_bignumint.h" +#include "s48_bignum.h" #include "bignum.h" #include "ratio.h" #include "float.h" diff --git a/native/misc.c b/native/misc.c index 268aa16d71..5581aa5eb9 100644 --- a/native/misc.c +++ b/native/misc.c @@ -42,3 +42,13 @@ void primitive_random_int(void) { dpush(tag_object(bignum(random()))); } + +void primitive_dump(void) +{ + /* Take an object, and print its memory. Later, return a vector */ + CELL obj = dpop(); + CELL size = object_size(obj); + int i; + for(i = 0; i < size; i += CELLS) + fprintf(stderr,"%x\n",get(UNTAG(obj) + i)); +} diff --git a/native/misc.h b/native/misc.h index 8118ef2959..5c8e8e64ac 100644 --- a/native/misc.h +++ b/native/misc.h @@ -4,3 +4,4 @@ void primitive_eq(void); void primitive_millis(void); void primitive_init_random(void); void primitive_random_int(void); +void primitive_dump(void); diff --git a/native/primitives.c b/native/primitives.c index 17893c7648..57b553afef 100644 --- a/native/primitives.c +++ b/native/primitives.c @@ -140,7 +140,8 @@ XT primitives[] = { primitive_size_of, primitive_profiling, primitive_word_call_count, - primitive_set_word_call_count + primitive_set_word_call_count, + primitive_dump }; CELL primitive_to_xt(CELL primitive) diff --git a/native/primitives.h b/native/primitives.h index baeaa8997b..9990a9d11e 100644 --- a/native/primitives.h +++ b/native/primitives.h @@ -1,4 +1,4 @@ extern XT primitives[]; -#define PRIMITIVE_COUNT 140 +#define PRIMITIVE_COUNT 141 CELL primitive_to_xt(CELL primitive); diff --git a/native/read.c b/native/read.c index 2277dde10a..6981fb29b8 100644 --- a/native/read.c +++ b/native/read.c @@ -87,7 +87,7 @@ bool can_read_line(PORT* port) pending_io_error(port); if(port->type != PORT_READ && port->type != PORT_RECV) - general_error(ERROR_INCOMPATIBLE_PORT,port); + general_error(ERROR_INCOMPATIBLE_PORT,tag_object(port)); if(port->line_ready) return true; @@ -184,7 +184,7 @@ bool can_read_count(PORT* port, FIXNUM count) pending_io_error(port); if(port->type != PORT_READ && port->type != PORT_RECV) - general_error(ERROR_INCOMPATIBLE_PORT,port); + general_error(ERROR_INCOMPATIBLE_PORT,tag_object(port)); if(port->line != F && CAN_READ_COUNT(port,count)) return true; @@ -237,7 +237,7 @@ void primitive_read_count_8(void) PORT* port = untag_port(dpop()); FIXNUM len = to_fixnum(dpop()); if(port->count != len) - critical_error("read# counts don't match",port); + critical_error("read# counts don't match",tag_object(port)); pending_io_error(port); diff --git a/native/relocate.c b/native/relocate.c index cda6281d90..3b5276c758 100644 --- a/native/relocate.c +++ b/native/relocate.c @@ -26,6 +26,7 @@ void relocate_object() break; case PORT_TYPE: fixup_port((PORT*)relocating); + break; } relocating += size; diff --git a/native/s48_bignum.c b/native/s48_bignum.c new file mode 100644 index 0000000000..33520052f2 --- /dev/null +++ b/native/s48_bignum.c @@ -0,0 +1,2011 @@ +/* -*-C-*- + +$Id$ + +Copyright (c) 1989-94 Massachusetts Institute of Technology + +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: + * - Add s48_ prefix to file names + * - Adapt s48_bignumint.h for Factor memory manager + */ + +#include "factor.h" +#include +#include +#include /* abort */ +#include + +/* Forward references */ +static int bignum_equal_p_unsigned(bignum_type, bignum_type); +static enum bignum_comparison bignum_compare_unsigned(bignum_type, bignum_type); +static bignum_type bignum_add_unsigned(bignum_type, bignum_type, int); +static bignum_type bignum_subtract_unsigned(bignum_type, bignum_type); +static bignum_type bignum_multiply_unsigned(bignum_type, bignum_type, int); +static bignum_type bignum_multiply_unsigned_small_factor + (bignum_type, bignum_digit_type, int); +static void bignum_destructive_scale_up(bignum_type, bignum_digit_type); +static void bignum_destructive_add(bignum_type, bignum_digit_type); +static void bignum_divide_unsigned_large_denominator + (bignum_type, bignum_type, bignum_type *, bignum_type *, int, int); +static void bignum_destructive_normalization(bignum_type, bignum_type, int); +static void bignum_destructive_unnormalization(bignum_type, int); +static void bignum_divide_unsigned_normalized(bignum_type, bignum_type, bignum_type); +static bignum_digit_type bignum_divide_subtract + (bignum_digit_type *, bignum_digit_type *, bignum_digit_type, + bignum_digit_type *); +static void bignum_divide_unsigned_medium_denominator + (bignum_type, bignum_digit_type, bignum_type *, bignum_type *, int, int); +static bignum_digit_type bignum_digit_divide + (bignum_digit_type, bignum_digit_type, bignum_digit_type, bignum_digit_type *); +static bignum_digit_type bignum_digit_divide_subtract + (bignum_digit_type, bignum_digit_type, bignum_digit_type, bignum_digit_type *); +static void bignum_divide_unsigned_small_denominator + (bignum_type, bignum_digit_type, bignum_type *, bignum_type *, int, int); +static bignum_digit_type bignum_destructive_scale_down + (bignum_type, bignum_digit_type); +static bignum_type bignum_remainder_unsigned_small_denominator + (bignum_type, bignum_digit_type, int); +static bignum_type bignum_digit_to_bignum(bignum_digit_type, int); +static bignum_type bignum_allocate(bignum_length_type, int); +static bignum_type bignum_allocate_zeroed(bignum_length_type, int); +static bignum_type bignum_shorten_length(bignum_type, bignum_length_type); +static bignum_type bignum_trim(bignum_type); +static bignum_type bignum_copy(bignum_type); +static bignum_type bignum_new_sign(bignum_type, int); +static bignum_type bignum_maybe_new_sign(bignum_type, int); +static void bignum_destructive_copy(bignum_type, bignum_type); +/* Unused +static void bignum_destructive_zero(bignum_type); +*/ + +/* Added for bitwise operations. */ +static bignum_type bignum_magnitude_ash(bignum_type arg1, long n); +static bignum_type bignum_pospos_bitwise_op(int op, bignum_type, bignum_type); +static bignum_type bignum_posneg_bitwise_op(int op, bignum_type, bignum_type); +static bignum_type bignum_negneg_bitwise_op(int op, bignum_type, bignum_type); +static void bignum_negate_magnitude(bignum_type); +static long bignum_unsigned_logcount(bignum_type arg); +static int bignum_unsigned_logbitp(int shift, bignum_type bignum); + +static ARRAY* s48_bignum_zero; +static ARRAY* s48_bignum_pos_one; +static ARRAY* s48_bignum_neg_one; + +/* Exports */ + +/* + * We have to allocate the cached constants slightly differently because + * they have to be registered with the GC, which requires that we have + * tagged pointers to them. + */ + +void +s48_bignum_make_cached_constants(void) +{ + bignum_type temp; + + s48_bignum_zero = bignum_allocate(0,0); + + s48_bignum_pos_one = bignum_allocate(1,0); + (BIGNUM_REF (s48_bignum_pos_one, 0)) = 1; + + s48_bignum_neg_one = bignum_allocate(1,0); + (BIGNUM_REF (s48_bignum_neg_one, 0)) = 1; +} + +int +s48_bignum_equal_p(bignum_type x, bignum_type 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 +s48_bignum_test(bignum_type bignum) +{ + return + ((BIGNUM_ZERO_P (bignum)) + ? bignum_comparison_equal + : (BIGNUM_NEGATIVE_P (bignum)) + ? bignum_comparison_less + : bignum_comparison_greater); +} + +enum bignum_comparison +s48_bignum_compare(bignum_type x, bignum_type 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)))); +} + +bignum_type +s48_bignum_add(bignum_type x, bignum_type y) +{ + return + ((BIGNUM_ZERO_P (x)) + ? (BIGNUM_MAYBE_COPY (y)) + : (BIGNUM_ZERO_P (y)) + ? (BIGNUM_MAYBE_COPY (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))))); +} + +bignum_type +s48_bignum_subtract(bignum_type x, bignum_type y) +{ + return + ((BIGNUM_ZERO_P (x)) + ? ((BIGNUM_ZERO_P (y)) + ? (BIGNUM_MAYBE_COPY (y)) + : (bignum_new_sign (y, (! (BIGNUM_NEGATIVE_P (y)))))) + : ((BIGNUM_ZERO_P (y)) + ? (BIGNUM_MAYBE_COPY (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)))))); +} + +bignum_type +s48_bignum_negate(bignum_type x) +{ + return + ((BIGNUM_ZERO_P (x)) + ? (BIGNUM_MAYBE_COPY (x)) + : (bignum_new_sign (x, (! (BIGNUM_NEGATIVE_P (x)))))); +} + +bignum_type +s48_bignum_multiply(bignum_type x, bignum_type y) +{ + 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 (BIGNUM_MAYBE_COPY (x)); + if (BIGNUM_ZERO_P (y)) + return (BIGNUM_MAYBE_COPY (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)); +} + +static int +bignum_divide(bignum_type numerator, bignum_type denominator, + bignum_type * quotient, bignum_type * remainder) +{ + if (BIGNUM_ZERO_P (denominator)) + return (1); + if (BIGNUM_ZERO_P (numerator)) + { + (*quotient) = (BIGNUM_MAYBE_COPY (numerator)); + (*remainder) = (BIGNUM_MAYBE_COPY (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) = (BIGNUM_MAYBE_COPY (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; + } + } + } + return (0); +} + +int +s48_bignum_divide(bignum_type numerator, bignum_type denominator, + void* quotient, void * remainder) +{ + return bignum_divide(numerator, denominator, + (bignum_type *)quotient, (bignum_type *)remainder); +} + +bignum_type +s48_bignum_quotient(bignum_type numerator, bignum_type denominator) +{ + if (BIGNUM_ZERO_P (denominator)) + return (BIGNUM_OUT_OF_BAND); + if (BIGNUM_ZERO_P (numerator)) + return (BIGNUM_MAYBE_COPY (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_type 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_type *) 0), + q_negative_p, 0); + else + bignum_divide_unsigned_medium_denominator + (numerator, digit, + ("ient), ((bignum_type *) 0), + q_negative_p, 0); + } + else + bignum_divide_unsigned_large_denominator + (numerator, denominator, + ("ient), ((bignum_type *) 0), + q_negative_p, 0); + return (quotient); + } + } + } +} + +bignum_type +s48_bignum_remainder(bignum_type numerator, bignum_type denominator) +{ + if (BIGNUM_ZERO_P (denominator)) + return (BIGNUM_OUT_OF_BAND); + if (BIGNUM_ZERO_P (numerator)) + return (BIGNUM_MAYBE_COPY (numerator)); + switch (bignum_compare_unsigned (numerator, denominator)) + { + case bignum_comparison_equal: + return (BIGNUM_ZERO ()); + case bignum_comparison_less: + return (BIGNUM_MAYBE_COPY (numerator)); + case bignum_comparison_greater: + default: /* to appease gcc -Wall */ + { + bignum_type 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_type *) 0), (&remainder), + 0, (BIGNUM_NEGATIVE_P (numerator))); + } + else + bignum_divide_unsigned_large_denominator + (numerator, denominator, + ((bignum_type *) 0), (&remainder), + 0, (BIGNUM_NEGATIVE_P (numerator))); + return (remainder); + } + } +} + +/* These procedures depend on the non-portable type `unsigned long'. + If your compiler doesn't support this type, either define the + switch `BIGNUM_NO_ULONG' to disable them (in "bignum.h"), or write + new versions that don't use this type. */ + +#ifndef BIGNUM_NO_ULONG + +bignum_type +s48_long_to_bignum(long n) +{ + int negative_p; + bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG]; + bignum_digit_type * end_digits = result_digits; + /* Special cases win when these small constants are cached. */ + if (n == 0) return (BIGNUM_ZERO ()); + if (n == 1) return (BIGNUM_ONE (0)); + if (n == -1) return (BIGNUM_ONE (1)); + { + unsigned long accumulator = ((negative_p = (n < 0)) ? (-n) : n); + do + { + (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); + accumulator >>= BIGNUM_DIGIT_LENGTH; + } + while (accumulator != 0); + } + { + bignum_type result = + (bignum_allocate ((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); + } +} + +long +s48_bignum_to_long(bignum_type bignum) +{ + if (BIGNUM_ZERO_P (bignum)) + return (0); + { + unsigned long accumulator = 0; + bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); + while (start < scan) + accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); + return ((BIGNUM_NEGATIVE_P (bignum)) ? (-((long)accumulator)) : accumulator); + } +} + +bignum_type +ulong_to_bignum(unsigned long n) +{ + bignum_digit_type result_digits [BIGNUM_DIGITS_FOR_LONG]; + bignum_digit_type * end_digits = result_digits; + /* Special cases win when these small constants are cached. */ + if (n == 0) return (BIGNUM_ZERO ()); + if (n == 1) return (BIGNUM_ONE (0)); + { + unsigned long accumulator = n; + do + { + (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK); + accumulator >>= BIGNUM_DIGIT_LENGTH; + } + while (accumulator != 0); + } + { + bignum_type result = + (bignum_allocate ((end_digits - result_digits), 0)); + 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); + } +} + +unsigned long +s48_bignum_to_ulong(bignum_type bignum) +{ + if (BIGNUM_ZERO_P (bignum)) + return (0); + { + unsigned long accumulator = 0; + bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); + while (start < scan) + accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan)); + return (accumulator); + } +} + +#endif /* not BIGNUM_NO_ULONG */ + +#define DTB_WRITE_DIGIT(factor) \ +{ \ + significand *= (factor); \ + digit = ((bignum_digit_type) significand); \ + (*--scan) = digit; \ + significand -= ((double) digit); \ +} + +bignum_type +s48_double_to_bignum(double x) +{ + 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_type result = (bignum_allocate (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 (1L << 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 + +double +s48_bignum_to_double(bignum_type bignum) +{ + if (BIGNUM_ZERO_P (bignum)) + return (0); + { + double accumulator = 0; + bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); + while (start < scan) + accumulator = ((accumulator * BIGNUM_RADIX) + (*--scan)); + return ((BIGNUM_NEGATIVE_P (bignum)) ? (-accumulator) : accumulator); + } +} + +int +s48_bignum_fits_in_word_p(bignum_type bignum, long word_length, + int twos_complement_p) +{ + unsigned int n_bits = (twos_complement_p ? (word_length - 1) : word_length); + BIGNUM_ASSERT (n_bits > 0); + { + bignum_length_type length = (BIGNUM_LENGTH (bignum)); + bignum_length_type max_digits = (BIGNUM_BITS_TO_DIGITS (n_bits)); + bignum_digit_type msd, max; + return + ((length < max_digits) || + ((length == max_digits) && + ((((msd = (BIGNUM_REF (bignum, (length - 1)))) < + (max = (1L << (n_bits - ((length - 1) * BIGNUM_DIGIT_LENGTH))))) || + (twos_complement_p && + (msd == max) && + (BIGNUM_NEGATIVE_P (bignum))))))); + } +} + +bignum_type +s48_bignum_length_in_bits(bignum_type bignum) +{ + if (BIGNUM_ZERO_P (bignum)) + return (BIGNUM_ZERO ()); + { + bignum_length_type index = ((BIGNUM_LENGTH (bignum)) - 1); + bignum_digit_type digit = (BIGNUM_REF (bignum, index)); + bignum_type result = (bignum_allocate (2, 0)); + (BIGNUM_REF (result, 0)) = index; + (BIGNUM_REF (result, 1)) = 0; + bignum_destructive_scale_up (result, BIGNUM_DIGIT_LENGTH); + while (digit > 0) + { + bignum_destructive_add (result, ((bignum_digit_type) 1)); + digit >>= 1; + } + return (bignum_trim (result)); + } +} + +bignum_type +s48_bignum_length_upper_limit(void) +{ + bignum_type result = (bignum_allocate (2, 0)); + (BIGNUM_REF (result, 0)) = 0; + (BIGNUM_REF (result, 1)) = BIGNUM_DIGIT_LENGTH; + return (result); +} + +bignum_type +s48_digit_stream_to_bignum(unsigned int n_digits, + unsigned int *producer(bignum_procedure_context), + bignum_procedure_context context, + unsigned int radix, + int negative_p) +{ + BIGNUM_ASSERT ((radix > 1) && (radix <= BIGNUM_RADIX_ROOT)); + if (n_digits == 0) + return (BIGNUM_ZERO ()); + if (n_digits == 1) + { + long digit = ((long) ((*producer) (context))); + return (s48_long_to_bignum (negative_p ? (- digit) : digit)); + } + { + bignum_length_type length; + { + unsigned int radix_copy = radix; + unsigned int log_radix = 0; + while (radix_copy > 0) + { + radix_copy >>= 1; + log_radix += 1; + } + /* This length will be at least as large as needed. */ + length = (BIGNUM_BITS_TO_DIGITS (n_digits * log_radix)); + } + { + bignum_type result = (bignum_allocate_zeroed (length, negative_p)); + while ((n_digits--) > 0) + { + bignum_destructive_scale_up (result, ((bignum_digit_type) radix)); + bignum_destructive_add + (result, ((bignum_digit_type) ((*producer) (context)))); + } + return (bignum_trim (result)); + } + } +} + +long +s48_bignum_max_digit_stream_radix(void) +{ + return (BIGNUM_RADIX_ROOT); +} + +/* Comparisons */ + +static int +bignum_equal_p_unsigned(bignum_type x, bignum_type 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); + } +} + +static enum bignum_comparison +bignum_compare_unsigned(bignum_type x, bignum_type 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 */ + +static bignum_type +bignum_add_unsigned(bignum_type x, bignum_type y, int negative_p) +{ + if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x))) + { + bignum_type z = x; + x = y; + y = z; + } + { + bignum_length_type x_length = (BIGNUM_LENGTH (x)); + bignum_type r = (bignum_allocate ((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 */ + +static bignum_type +bignum_subtract_unsigned(bignum_type x, bignum_type y) +{ + int negative_p; + switch (bignum_compare_unsigned (x, y)) + { + case bignum_comparison_equal: + return (BIGNUM_ZERO ()); + case bignum_comparison_less: + { + bignum_type z = x; + x = y; + y = z; + } + negative_p = 1; + break; + case bignum_comparison_greater: + negative_p = 0; + break; + } + { + bignum_length_type x_length = (BIGNUM_LENGTH (x)); + bignum_type r = (bignum_allocate (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 */ + +static bignum_type +bignum_multiply_unsigned(bignum_type x, bignum_type y, int negative_p) +{ + if ((BIGNUM_LENGTH (y)) > (BIGNUM_LENGTH (x))) + { + bignum_type z = x; + x = y; + y = z; + } + { + 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_type r = + (bignum_allocate_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 + } +} + +static bignum_type +bignum_multiply_unsigned_small_factor(bignum_type x, bignum_digit_type y, + int negative_p) +{ + bignum_length_type length_x = (BIGNUM_LENGTH (x)); + bignum_type p = (bignum_allocate ((length_x + 1), negative_p)); + bignum_destructive_copy (x, p); + (BIGNUM_REF (p, length_x)) = 0; + bignum_destructive_scale_up (p, y); + return (bignum_trim (p)); +} + +static void +bignum_destructive_scale_up(bignum_type bignum, bignum_digit_type factor) +{ + bignum_digit_type carry = 0; + bignum_digit_type * scan = (BIGNUM_START_PTR (bignum)); + bignum_digit_type two_digits; + bignum_digit_type product_low; +#define product_high carry + bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum))); + 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 +} + +static void +bignum_destructive_add(bignum_type bignum, bignum_digit_type n) +{ + bignum_digit_type * scan = (BIGNUM_START_PTR (bignum)); + 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); + } +} + +/* 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". */ + +static void +bignum_divide_unsigned_large_denominator(bignum_type numerator, + bignum_type denominator, + bignum_type * quotient, + bignum_type * remainder, + int q_negative_p, + int r_negative_p) +{ + bignum_length_type length_n = ((BIGNUM_LENGTH (numerator)) + 1); + bignum_length_type length_d = (BIGNUM_LENGTH (denominator)); + bignum_type q = + ((quotient != ((bignum_type *) 0)) + ? (bignum_allocate ((length_n - length_d), q_negative_p)) + : BIGNUM_OUT_OF_BAND); + bignum_type u = (bignum_allocate (length_n, r_negative_p)); + int shift = 0; + BIGNUM_ASSERT (length_d > 1); + { + bignum_digit_type v1 = (BIGNUM_REF ((denominator), (length_d - 1))); + while (v1 < (BIGNUM_RADIX / 2)) + { + v1 <<= 1; + shift += 1; + } + } + if (shift == 0) + { + bignum_destructive_copy (numerator, u); + (BIGNUM_REF (u, (length_n - 1))) = 0; + bignum_divide_unsigned_normalized (u, denominator, q); + } + else + { + bignum_type v = (bignum_allocate (length_d, 0)); + bignum_destructive_normalization (numerator, u, shift); + bignum_destructive_normalization (denominator, v, shift); + bignum_divide_unsigned_normalized (u, v, q); + BIGNUM_DEALLOCATE (v); + if (remainder != ((bignum_type *) 0)) + bignum_destructive_unnormalization (u, shift); + } + if (quotient != ((bignum_type *) 0)) + (*quotient) = (bignum_trim (q)); + if (remainder != ((bignum_type *) 0)) + (*remainder) = (bignum_trim (u)); + else + BIGNUM_DEALLOCATE (u); + return; +} + +static void +bignum_divide_unsigned_normalized(bignum_type u, bignum_type v, bignum_type 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; + 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 +} + +static bignum_digit_type +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); +} + +static void +bignum_divide_unsigned_medium_denominator(bignum_type numerator, + bignum_digit_type denominator, + bignum_type * quotient, + bignum_type * remainder, + int q_negative_p, + int r_negative_p) +{ + bignum_length_type length_n = (BIGNUM_LENGTH (numerator)); + bignum_length_type length_q; + bignum_type q; + int shift = 0; + /* Because `bignum_digit_divide' requires a normalized denominator. */ + while (denominator < (BIGNUM_RADIX / 2)) + { + denominator <<= 1; + shift += 1; + } + if (shift == 0) + { + length_q = length_n; + q = (bignum_allocate (length_q, q_negative_p)); + bignum_destructive_copy (numerator, q); + } + else + { + length_q = (length_n + 1); + q = (bignum_allocate (length_q, q_negative_p)); + bignum_destructive_normalization (numerator, q, 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; + if (quotient != ((bignum_type *) 0)) + { + while (start < scan) + { + r = (bignum_digit_divide (r, (*--scan), denominator, (&qj))); + (*scan) = qj; + } + (*quotient) = (bignum_trim (q)); + } + else + { + while (start < scan) + r = (bignum_digit_divide (r, (*--scan), denominator, (&qj))); + BIGNUM_DEALLOCATE (q); + } + if (remainder != ((bignum_type *) 0)) + { + if (shift != 0) + r >>= shift; + (*remainder) = (bignum_digit_to_bignum (r, r_negative_p)); + } + } + return; +} + +static void +bignum_destructive_normalization(bignum_type source, bignum_type 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 = ((1L << 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; +} + +static void +bignum_destructive_unnormalization(bignum_type bignum, int shift_right) +{ + bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); + bignum_digit_type digit; + bignum_digit_type carry = 0; + int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right); + bignum_digit_type mask = ((1L << 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]))); \ +} + +static bignum_digit_type +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; \ + } \ +} + +static bignum_digit_type +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 + +static void +bignum_divide_unsigned_small_denominator(bignum_type numerator, + bignum_digit_type denominator, + bignum_type * quotient, + bignum_type * remainder, + int q_negative_p, + int r_negative_p) +{ + bignum_type q = (bignum_new_sign (numerator, q_negative_p)); + bignum_digit_type r = (bignum_destructive_scale_down (q, denominator)); + (*quotient) = (bignum_trim (q)); + if (remainder != ((bignum_type *) 0)) + (*remainder) = (bignum_digit_to_bignum (r, r_negative_p)); + 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. */ + +static bignum_digit_type +bignum_destructive_scale_down(bignum_type bignum, 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 (bignum)); + bignum_digit_type * scan = (start + (BIGNUM_LENGTH (bignum))); + 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 +} + +static bignum_type +bignum_remainder_unsigned_small_denominator( + bignum_type 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)); +} + +static bignum_type +bignum_digit_to_bignum(bignum_digit_type digit, int negative_p) +{ + if (digit == 0) + return (BIGNUM_ZERO ()); + else + { + bignum_type result = (bignum_allocate (1, negative_p)); + (BIGNUM_REF (result, 0)) = digit; + return (result); + } +} + +/* Allocation */ + +static bignum_type +bignum_allocate(bignum_length_type length, int negative_p) +{ + BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX)); + { + bignum_type result = (BIGNUM_ALLOCATE (length)); + BIGNUM_SET_NEGATIVE_P (result, negative_p); + return (result); + } +} + +static bignum_type +bignum_allocate_zeroed(bignum_length_type length, int negative_p) +{ + BIGNUM_ASSERT ((length >= 0) || (length < BIGNUM_RADIX)); + { + bignum_type result = (BIGNUM_ALLOCATE (length)); + bignum_digit_type * scan = (BIGNUM_START_PTR (result)); + bignum_digit_type * end = (scan + length); + BIGNUM_SET_NEGATIVE_P (result, negative_p); + while (scan < end) + (*scan++) = 0; + return (result); + } +} + +static bignum_type +bignum_shorten_length(bignum_type bignum, bignum_length_type length) +{ + bignum_length_type current_length = (BIGNUM_LENGTH (bignum)); + BIGNUM_ASSERT ((length >= 0) || (length <= current_length)); + if (length < current_length) + { + BIGNUM_REDUCE_LENGTH (bignum, bignum, length); + BIGNUM_SET_NEGATIVE_P (bignum, (length != 0) && (BIGNUM_NEGATIVE_P (bignum))); + } + return (bignum); +} + +static bignum_type +bignum_trim(bignum_type bignum) +{ + bignum_digit_type * start = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * end = (start + (BIGNUM_LENGTH (bignum))); + bignum_digit_type * scan = end; + while ((start <= scan) && ((*--scan) == 0)) + ; + scan += 1; + if (scan < end) + { + bignum_length_type length = (scan - start); + BIGNUM_REDUCE_LENGTH (bignum, bignum, length); + BIGNUM_SET_NEGATIVE_P (bignum, (length != 0) && (BIGNUM_NEGATIVE_P (bignum))); + } + return (bignum); +} + +/* Copying */ + +static bignum_type +bignum_copy(bignum_type source) +{ + bignum_type target = + (bignum_allocate ((BIGNUM_LENGTH (source)), (BIGNUM_NEGATIVE_P (source)))); + bignum_destructive_copy (source, target); + return (target); +} + +static bignum_type +bignum_new_sign(bignum_type bignum, int negative_p) +{ + bignum_type result = + (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p)); + bignum_destructive_copy (bignum, result); + return (result); +} + +static bignum_type +bignum_maybe_new_sign(bignum_type bignum, int negative_p) +{ +#ifndef BIGNUM_FORCE_NEW_RESULTS + if ((BIGNUM_NEGATIVE_P (bignum)) ? negative_p : (! negative_p)) + return (bignum); + else +#endif /* not BIGNUM_FORCE_NEW_RESULTS */ + { + bignum_type result = + (bignum_allocate ((BIGNUM_LENGTH (bignum)), negative_p)); + bignum_destructive_copy (bignum, result); + return (result); + } +} + +static void +bignum_destructive_copy(bignum_type source, bignum_type 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; +} + +/* Unused +static void +bignum_destructive_zero(bignum_type bignum) +{ + bignum_digit_type * scan = (BIGNUM_START_PTR (bignum)); + bignum_digit_type * end = (scan + (BIGNUM_LENGTH (bignum))); + while (scan < end) + (*scan++) = 0; + return; +} +*/ + +/* + * Added bitwise operations (and oddp). + */ + +int +s48_bignum_oddp (bignum_type bignum) +{ + return (BIGNUM_LENGTH (bignum) > 0) && (BIGNUM_REF (bignum, 0) & 1); +} + +bignum_type +s48_bignum_bitwise_not(bignum_type x) +{ + return s48_bignum_subtract(BIGNUM_ONE(1), x); +} + +bignum_type +s48_bignum_arithmetic_shift(bignum_type arg1, long n) +{ + if (BIGNUM_NEGATIVE_P(arg1) && n < 0) + return + s48_bignum_bitwise_not(bignum_magnitude_ash(s48_bignum_bitwise_not(arg1), + n)); + else + return bignum_magnitude_ash(arg1, n); +} + +/* + * This uses a `long'-returning bignum_length_in_bits() which we don't have. +long +s48_bignum_integer_length(bignum_type arg1) +{ + return((BIGNUM_NEGATIVE_P (arg1)) + ? bignum_length_in_bits (s48_bignum_bitwise_not (arg1)) + : bignum_length_in_bits (arg1)); +} +*/ + +long +s48_bignum_bit_count(bignum_type arg1) +{ + return((BIGNUM_NEGATIVE_P (arg1)) + ? bignum_unsigned_logcount (s48_bignum_bitwise_not (arg1)) + : bignum_unsigned_logcount (arg1)); +} + +#define AND_OP 0 +#define IOR_OP 1 +#define XOR_OP 2 + +bignum_type +s48_bignum_bitwise_and(bignum_type arg1, bignum_type 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) + ); +} + +bignum_type +s48_bignum_bitwise_ior(bignum_type arg1, bignum_type 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) + ); +} + +bignum_type +s48_bignum_bitwise_xor(bignum_type arg1, bignum_type 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) + ); +} + +/* ash for the magnitude */ +/* assume arg1 is a big number, n is a long */ +static bignum_type +bignum_magnitude_ash(bignum_type arg1, long n) +{ + bignum_type result; + bignum_digit_type *scan1; + bignum_digit_type *scanr; + bignum_digit_type *end; + + long digit_offset,bit_offset; + + if (BIGNUM_ZERO_P (arg1)) return (arg1); + + if (n > 0) { + digit_offset = n / BIGNUM_DIGIT_LENGTH; + bit_offset = n % BIGNUM_DIGIT_LENGTH; + + result = bignum_allocate_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 = bignum_allocate_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; + + return (bignum_trim (result)); +} + +static bignum_type +bignum_pospos_bitwise_op(int op, bignum_type arg1, bignum_type arg2) +{ + bignum_type result; + 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); + + result = bignum_allocate(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; + /* + fprintf(stderr, "[pospos op = %d, i = %ld, d1 = %lx, d2 = %lx]\n", + op, endr - scanr, digit1, digit2); + */ + *scanr++ = (op == 0) ? digit1 & digit2 : + (op == 1) ? digit1 | digit2 : + digit1 ^ digit2; + } + return bignum_trim(result); +} + +static bignum_type +bignum_posneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2) +{ + bignum_type result; + 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; + + result = bignum_allocate(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); +} + +static bignum_type +bignum_negneg_bitwise_op(int op, bignum_type arg1, bignum_type arg2) +{ + bignum_type result; + 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; + + result = bignum_allocate(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 == 0) ? digit1 & digit2 : + (op == 1) ? digit1 | digit2 : + digit1 ^ digit2; + } + + if (neg_p) + bignum_negate_magnitude(result); + + return bignum_trim(result); +} + +static void +bignum_negate_magnitude(bignum_type 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; + } +} + +static long +bignum_unsigned_logcount(bignum_type arg) +{ + + bignum_digit_type *scan; + bignum_digit_type *end; + bignum_digit_type digit; + + /* sufficient for any reasonable big number */ + long result; + int i; + + if (BIGNUM_ZERO_P (arg)) return (0L); + + scan = BIGNUM_START_PTR (arg); + end = scan + BIGNUM_LENGTH (arg); + result = 0L; + + while (scan < end) { + digit = *scan++ & BIGNUM_DIGIT_MASK; + for (i = 0; i++ < BIGNUM_DIGIT_LENGTH; digit = digit >> 1L) + result += digit & 1L; + } + + return (result); +} + +int +bignum_logbitp(int shift, bignum_type arg) +{ + return((BIGNUM_NEGATIVE_P (arg)) + ? !bignum_unsigned_logbitp (shift, s48_bignum_bitwise_not (arg)) + : bignum_unsigned_logbitp (shift,arg)); +} + +static int +bignum_unsigned_logbitp(int shift, bignum_type bignum) +{ + bignum_length_type len = (BIGNUM_LENGTH (bignum)); + bignum_digit_type digit; + int index = shift / BIGNUM_DIGIT_LENGTH; + int p; + if (index >= len) + return 0; + digit = (BIGNUM_REF (bignum, index)); + p = shift % BIGNUM_DIGIT_LENGTH; + return digit & (1 << p); +} + diff --git a/native/s48_bignum.h b/native/s48_bignum.h new file mode 100644 index 0000000000..8717f534b0 --- /dev/null +++ b/native/s48_bignum.h @@ -0,0 +1,94 @@ +/* -*-C-*- + +$Id$ + +Copyright (c) 1989-1992 Massachusetts Institute of Technology + +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. */ + +/* External Interface to Bignum Code */ + +/* The `unsigned long' type is used for the conversion procedures + `bignum_to_long' and `long_to_bignum'. Older implementations of C + don't support this type; if you have such an implementation you can + disable these procedures using the following flag (alternatively + you could write alternate versions that don't require this type). */ +/* #define BIGNUM_NO_ULONG */ + +typedef ARRAY * bignum_type; +#define BIGNUM_OUT_OF_BAND ((bignum_type) 0) + +enum bignum_comparison +{ + bignum_comparison_equal = 0, + bignum_comparison_less = -1, + bignum_comparison_greater = 1 +}; + +typedef void * bignum_procedure_context; +extern int s48_bignum_equal_p(bignum_type, bignum_type); +extern enum bignum_comparison s48_bignum_test(bignum_type); +extern enum bignum_comparison s48_bignum_compare(bignum_type, bignum_type); +extern bignum_type s48_bignum_add(bignum_type, bignum_type); +extern bignum_type s48_bignum_subtract(bignum_type, bignum_type); +extern bignum_type s48_bignum_negate(bignum_type); +extern bignum_type s48_bignum_multiply(bignum_type, bignum_type); +extern int s48_bignum_divide(bignum_type numerator, bignum_type denominator, + void * quotient, void * remainder); +extern bignum_type s48_bignum_quotient(bignum_type, bignum_type); +extern bignum_type s48_bignum_remainder(bignum_type, bignum_type); +extern bignum_type s48_long_to_bignum(long); +extern bignum_type s48_ulong_to_bignum(unsigned long); +extern long s48_bignum_to_long(bignum_type); +extern unsigned long s48_bignum_to_ulong(bignum_type); +extern bignum_type s48_double_to_bignum(double); +extern double s48_bignum_to_double(bignum_type); +extern int s48_bignum_fits_in_word_p(bignum_type, long word_length, + int twos_complement_p); +extern bignum_type s48_bignum_length_in_bits(bignum_type); +extern bignum_type s48_bignum_length_upper_limit(void); +extern bignum_type s48_digit_stream_to_bignum + (unsigned int n_digits, + unsigned int (*producer(bignum_procedure_context)), + bignum_procedure_context context, + unsigned int radix, + int negative_p); +extern long s48_bignum_max_digit_stream_radix(void); + +/* Added bitwise operators. */ + +extern bignum_type s48_bignum_bitwise_not(bignum_type), + s48_bignum_arithmetic_shift(bignum_type, long), + s48_bignum_bitwise_and(bignum_type, bignum_type), + s48_bignum_bitwise_ior(bignum_type, bignum_type), + s48_bignum_bitwise_xor(bignum_type, bignum_type); + +extern int s48_bignum_oddp(bignum_type); +extern long s48_bignum_bit_count(bignum_type); diff --git a/native/s48_bignumint.h b/native/s48_bignumint.h new file mode 100644 index 0000000000..a23ebf0ff2 --- /dev/null +++ b/native/s48_bignumint.h @@ -0,0 +1,127 @@ +/* -*-C-*- + +$Id$ + +Copyright (c) 1989-1992 Massachusetts Institute of Technology + +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. */ + +/* Internal Interface to Bignum Code */ +#undef BIGNUM_ZERO_P +#undef BIGNUM_NEGATIVE_P + +/* The memory model is based on the following definitions, and on the + definition of the type `bignum_type'. The only other special + definition is `CHAR_BIT', which is defined in the Ansi C header + file "limits.h". */ + +typedef CELL bignum_digit_type; +typedef CELL bignum_length_type; + +/* BIGNUM_ALLOCATE allocates a (length + 1)-element array of + `bignum_digit_type'; deallocation is the responsibility of the + user (in Factor, the garbage collector handles this). */ +#define BIGNUM_ALLOCATE(length_in_digits) \ + allot_array(BIGNUM_TYPE,length_in_digits + 1) + +/* BIGNUM_TO_POINTER casts a bignum object to a digit array pointer. */ +#define BIGNUM_TO_POINTER(bignum) ((CELL*)AREF(bignum,0)) + +/* BIGNUM_REDUCE_LENGTH allows the memory system to reclaim some + space when a bignum's length is reduced from its original value. */ +#define BIGNUM_REDUCE_LENGTH(target, source, length) \ + target = shrink_array(source, length) +extern ARRAY* shrink_array(ARRAY* array, CELL capacity); + +/* BIGNUM_DEALLOCATE is called when disposing of bignums which are + created as intermediate temporaries; Scheme doesn't need this. */ +#define BIGNUM_DEALLOCATE(bignum) + +/* If BIGNUM_FORCE_NEW_RESULTS is defined, all bignum-valued operations + return freshly-allocated results. This is useful for some kinds of + memory deallocation strategies. */ +/* #define BIGNUM_FORCE_NEW_RESULTS */ + +/* BIGNUM_EXCEPTION is invoked to handle assertion violations. */ +#define BIGNUM_EXCEPTION abort + + +#define BIGNUM_DIGIT_LENGTH (((sizeof (bignum_digit_type)) * CHAR_BIT) - 2) +#define BIGNUM_HALF_DIGIT_LENGTH (BIGNUM_DIGIT_LENGTH / 2) +#define BIGNUM_RADIX (((unsigned long) 1) << BIGNUM_DIGIT_LENGTH) +#define BIGNUM_RADIX_ROOT (((unsigned long) 1) << BIGNUM_HALF_DIGIT_LENGTH) +#define BIGNUM_DIGIT_MASK (BIGNUM_RADIX - 1) +#define BIGNUM_HALF_DIGIT_MASK (BIGNUM_RADIX_ROOT - 1) + +#define BIGNUM_START_PTR(bignum) \ + ((BIGNUM_TO_POINTER (bignum)) + 1) + +#define BIGNUM_LENGTH(bignum) (bignum->capacity - 1) + +#define BIGNUM_NEGATIVE_P(bignum) (AREF(bignum,0) != 0) +#define BIGNUM_SET_NEGATIVE_P(bignum,neg) put(AREF(bignum,0),neg) + +#define BIGNUM_ZERO_P(bignum) \ + ((BIGNUM_LENGTH (bignum)) == 0) + +#define BIGNUM_REF(bignum, index) \ + (* ((BIGNUM_START_PTR (bignum)) + (index))) + +#ifdef BIGNUM_FORCE_NEW_RESULTS +#define BIGNUM_MAYBE_COPY bignum_copy +#else +#define BIGNUM_MAYBE_COPY(bignum) bignum +#endif + +/* These definitions are here to facilitate caching of the constants + 0, 1, and -1. */ +#define BIGNUM_ZERO() s48_bignum_zero +#define BIGNUM_ONE(neg_p) \ + (neg_p ? s48_bignum_pos_one : s48_bignum_neg_one) + +#define HD_LOW(digit) ((digit) & BIGNUM_HALF_DIGIT_MASK) +#define HD_HIGH(digit) ((digit) >> BIGNUM_HALF_DIGIT_LENGTH) +#define HD_CONS(high, low) (((high) << BIGNUM_HALF_DIGIT_LENGTH) | (low)) + +#define BIGNUM_BITS_TO_DIGITS(n) \ + (((n) + (BIGNUM_DIGIT_LENGTH - 1)) / BIGNUM_DIGIT_LENGTH) + +#define BIGNUM_DIGITS_FOR_LONG \ + (BIGNUM_BITS_TO_DIGITS ((sizeof (long)) * CHAR_BIT)) + +#ifndef BIGNUM_DISABLE_ASSERTION_CHECKS + +#define BIGNUM_ASSERT(expression) \ +{ \ + if (! (expression)) \ + BIGNUM_EXCEPTION (); \ +} + +#endif /* not BIGNUM_DISABLE_ASSERTION_CHECKS */ diff --git a/native/types.c b/native/types.c index d14533e944..53acec103e 100644 --- a/native/types.c +++ b/native/types.c @@ -77,6 +77,7 @@ CELL untagged_object_size(CELL pointer) size = CELLS * 2; break; case ARRAY_TYPE: + case BIGNUM_TYPE: size = ASIZE(pointer); break; case VECTOR_TYPE: @@ -88,9 +89,6 @@ CELL untagged_object_size(CELL pointer) case SBUF_TYPE: size = sizeof(SBUF); break; - case BIGNUM_TYPE: - size = sizeof(BIGNUM); - break; case FLOAT_TYPE: size = sizeof(FLOAT); break;