/* numbers-c.c * * Copyright 2010-2015 The CHICKEN Team * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are * met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * 2. Redistributions in binary form must reproduce the above * copyright notice, this list of conditions and the following * disclaimer in the documentation and/or other materials provided * with the distribution. * * 3. The name of the author may not be used to endorse or promote * products derived from this software without specific prior * written permission. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. */ #include #include #include /* frexp() */ #define nmax(x, y) ((x) > (y) ? (x) : (y)) /* From runtime.c */ #define nmin(x, y) ((x) < (y) ? (x) : (y)) /* From runtime.c */ #define free_tmp_bignum(b) C_free((void *)(b)) static void *tags; #include "numbers-c.h" /* The bignum digit representation is fullword- little endian, so on * LE machines the halfdigits are numbered in the same order. On BE * machines, we must swap the odd and even positions. */ #ifdef C_BIG_ENDIAN #define C_uhword_ref(x, p) ((C_uhword *)(x))[(p)^1] #else #define C_uhword_ref(x, p) ((C_uhword *)(x))[(p)] #endif #define C_uhword_set(x, p, d) (C_uhword_ref(x,p) = (d)) static C_word init_tags(___scheme_value tagvec); static void bignum_negate_2(C_word c, C_word self, C_word new_big) C_noret; static C_word rat_cmp(C_word x, C_word y); static C_word basic_cmp(C_word x, C_word y, char *loc, int eqp); static void allocate_bignum_2(C_word c, C_word self, C_word bigvec) C_noret; static C_word allocate_tmp_bignum(C_word size, C_word negp, C_word initp); static C_uword bignum_digits_destructive_scale_up_with_carry(C_uword *start, C_uword *end, C_uword factor, C_uword carry); static C_uword bignum_digits_destructive_scale_down(C_uword *start, C_uword *end, C_uword denominator); static C_uword bignum_digits_destructive_shift_right(C_uword *start, C_uword *end, int shift_right, int negp); static C_uword bignum_digits_destructive_shift_left(C_uword *start, C_uword *end, int shift_left); static void bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp) C_noret; static void bignum_plus_unsigned_2(C_word c, C_word self, C_word result) C_noret; static C_word int_flo_cmp(C_word intnum, C_word flonum); static C_word flo_int_cmp(C_word flonum, C_word intnum); static C_word rat_flo_cmp(C_word ratnum, C_word flonum); static C_word flo_rat_cmp(C_word flonum, C_word ratnum); static int bignum_cmp_unsigned(C_word x, C_word y); static void bignum_minus_unsigned(C_word k, C_word x, C_word y) C_noret; static void bignum_minus_unsigned_2(C_word c, C_word self, C_word result) C_noret; static void integer_times_2(C_word c, C_word self, C_word new_big) C_noret; static C_regparm void bignum_digits_multiply(C_word x, C_word y, C_word result); static void bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp) C_noret; static void bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result) C_noret; static void digits_to_integer_2(C_word c, C_word self, C_word result) C_noret; static C_regparm C_word str_to_bignum(C_word bignum, char *str, char *str_end, int radix); static void bignum_to_str_2(C_word c, C_word self, C_word string) C_noret; static void fabs_frexp_to_digits(C_uword exp, double sign, C_uword *start, C_uword *scan); static C_word flo_to_tmp_bignum(C_word x); static void flo_to_int_2(C_word c, C_word self, C_word result) C_noret; static void bignum_actual_shift(C_word c, C_word self, C_word result) C_noret; static void bignum_actual_extraction(C_word c, C_word self, C_word result) C_noret; static void bignum_random_2(C_word c, C_word self, C_word result) C_noret; static void bignum_bitwise_and_2(C_word c, C_word self, C_word result) C_noret; static void bignum_bitwise_ior_2(C_word c, C_word self, C_word result) C_noret; static void bignum_bitwise_xor_2(C_word c, C_word self, C_word result) C_noret; static void bignum_digits_destructive_negate(C_word result); static C_regparm void basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_r, C_word return_q) C_noret; static C_regparm void integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) C_noret; static C_regparm void bignum_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) C_noret; static void divrem_intflo_2(C_word c, C_word self, ...) C_noret; static void bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big) C_noret; static void bignum_negneg_bitwise_op(C_word c, C_word self, C_word result) C_noret; static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result) C_noret; static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result) C_noret; static C_word bignum_remainder_unsigned_halfdigit(C_word num, C_word den); static void bignum_destructive_divide_unsigned_small(C_word c, C_word self, C_word quotient); static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient) C_noret; static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder) C_noret; static void bignum_destructive_divide_normalized(C_word big_u, C_word big_v, C_word big_q); static void barf(int code, char *loc, ...) C_noret; static void try_extended_number(char *ext_proc_name, C_word c, C_word k, ...) C_noret; /* Use high numbers for when CHICKEN 4 grows more error codes! */ #define C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR 99948 #define C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR 99949 /* XXX THIS IS DUPLICATED HERE FROM runtime.c, but should be ripped out */ static void barf(int code, char *loc, ...) { char *msg; int c, i; va_list v; /* Just take a random size that will "always" fit... */ C_word err, ab[C_SIZEOF_STRING(64)], *a = ab; err = C_lookup_symbol(C_intern2(&a, C_text("numbers#@error-hook"))); switch(code) { case C_BAD_ARGUMENT_COUNT_ERROR: msg = C_text("bad argument count"); c = 3; break; case C_BAD_MINIMUM_ARGUMENT_COUNT_ERROR: msg = C_text("too few arguments"); c = 3; break; case C_BAD_ARGUMENT_TYPE_ERROR: msg = C_text("bad argument type"); c = 1; break; case C_DIVISION_BY_ZERO_ERROR: msg = C_text("division by zero"); c = 0; break; case C_OUT_OF_RANGE_ERROR: msg = C_text("out of range"); c = 2; break; case C_CANT_REPRESENT_INEXACT_ERROR: msg = C_text("inexact number cannot be represented as an exact number"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_FIXNUM_ERROR: msg = C_text("bad argument type - not a fixnum"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR: msg = C_text("bad argument type - not a number"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR: msg = C_text("bad argument type - not an integer"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR: msg = C_text("bad argument type - not an unsigned integer"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_FLONUM_ERROR: msg = C_text("bad argument type - not a flonum"); c = 1; break; case C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR: msg = C_text("bad argument type - invalid base"); c = 1; break; case C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR: msg = C_text("bad argument type - not a real"); c = 1; break; case C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR: msg = C_text("bad argument type - complex number has no ordering"); c = 1; break; default: fprintf(stderr, "Unknown error"); abort(); } if(!C_immediatep(err)) { C_save(C_fix(code)); C_save(C_intern2(&a, loc)); va_start(v, loc); i = c; while(i--) C_save(va_arg(v, C_word)); va_end(v); /* No continuation is passed: '##sys#error-hook' may not return: */ C_do_apply(c + 2, C_block_item(err, 0), C_SCHEME_UNDEFINED); } else { fprintf(stderr, "No error hook!"); abort(); } } void C_not_an_integer_error(char *loc, C_word x) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, loc, x); } /* This exists because numbers_div_by_zero_error doesn't exist in CHICKENs < 4.9.0 */ void numbers_div_by_zero_error(char *loc) { barf(C_DIVISION_BY_ZERO_ERROR, loc); } /* Never use extended number hook procedure names longer than this! */ /* Current longest name: numbers#@bignum-2-divrem-burnikel-ziegler */ #define MAX_EXTNUM_HOOK_NAME 64 /* This exists so that we don't have to create any extra closures */ static void try_extended_number(char *ext_proc_name, C_word c, C_word k, ...) { static C_word ab[C_SIZEOF_STRING(MAX_EXTNUM_HOOK_NAME)]; int i; va_list v; C_word ext_proc_sym, ext_proc = C_SCHEME_FALSE, *a = ab; ext_proc_sym = C_lookup_symbol(C_intern2(&a, ext_proc_name)); if(!C_immediatep(ext_proc_sym)) ext_proc = C_block_item(ext_proc_sym, 0); if (!C_immediatep(ext_proc) && C_closurep(ext_proc)) { va_start(v, k); i = c - 1; while(i--) C_save(va_arg(v, C_word)); va_end(v); C_do_apply(c - 1, ext_proc, k); } else { /* TODO: Convert to barf(), add new error code */ fprintf(stderr, "No extended number hook for %s!\n", ext_proc_name); abort(); } } static C_word init_tags(___scheme_value tagvec) { tags = CHICKEN_new_gc_root(); CHICKEN_gc_root_set(tags, tagvec); return C_SCHEME_UNDEFINED; } #define ratnum_type_tag C_block_item(CHICKEN_gc_root_ref(tags), RAT_TAG) #define compnum_type_tag C_block_item(CHICKEN_gc_root_ref(tags), COMP_TAG) C_regparm C_word C_fcall C_i_numbers_numberp(C_word x) { return C_mk_bool((x & C_FIXNUM_BIT) || (!C_immediatep(x) && (C_block_header(x) == C_FLONUM_TAG || C_IS_BIGNUM_TYPE(x) || (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == ratnum_type_tag || C_block_item(x, 0) == compnum_type_tag))))); } /* TODO: Rename to C_i_integerp */ C_regparm C_word C_fcall C_i_numbers_integerp(C_word x) { return C_mk_bool((x & C_FIXNUM_BIT) || (!C_immediatep(x) && ((C_block_header(x) == C_FLONUM_TAG && C_truep(C_u_i_fpintegerp_fixed(x))) || C_IS_BIGNUM_TYPE(x)))); } C_inline C_word basic_eqvp(C_word x, C_word y) { return (x == y || (!C_immediatep(x) && !C_immediatep(y) && C_block_header(x) == C_block_header(y) && ((C_block_header(x) == C_FLONUM_TAG && C_flonum_magnitude(x) == C_flonum_magnitude(y)) || (C_IS_BIGNUM_TYPE(x) && C_u_i_bignum_cmp(x, y) == C_fix(0))))); } C_regparm C_word C_fcall C_i_basic_positivep(C_word x) { if (x & C_FIXNUM_BIT) return C_u_i_fixnum_positivep(x); else if (C_immediatep(x)) barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "positive?", x); else if (C_block_header(x) == C_FLONUM_TAG) return C_mk_bool(C_flonum_magnitude(x) > 0.0); else if (C_IS_BIGNUM_TYPE(x)) return C_mk_nbool(C_bignum_negativep(x)); else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == ratnum_type_tag)) return C_u_i_integer_positivep(C_block_item(x, 1)); else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == compnum_type_tag)) barf(C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR, "positive?", x); else barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "positive?", x); } C_regparm C_word C_fcall C_u_i_integer_positivep(C_word x) { if (x & C_FIXNUM_BIT) return C_u_i_fixnum_positivep(x); else return C_mk_nbool(C_bignum_negativep(x)); } C_regparm C_word C_fcall C_i_basic_negativep(C_word x) { if (x & C_FIXNUM_BIT) return C_u_i_fixnum_negativep(x); else if (C_immediatep(x)) barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "negative?", x); else if (C_block_header(x) == C_FLONUM_TAG) return C_mk_bool(C_flonum_magnitude(x) < 0.0); else if (C_IS_BIGNUM_TYPE(x)) return C_mk_bool(C_bignum_negativep(x)); else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == ratnum_type_tag)) return C_u_i_integer_negativep(C_block_item(x, 1)); else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == compnum_type_tag)) barf(C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR, "negative?", x); else barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "negative?", x); } C_regparm C_word C_fcall C_u_i_integer_negativep(C_word x) { if (x & C_FIXNUM_BIT) return C_u_i_fixnum_negativep(x); else return C_mk_bool(C_bignum_negativep(x)); } /* TODO: Rename to C_i_eqvp */ C_regparm C_word C_fcall C_i_numbers_eqvp(C_word x, C_word y) { return C_mk_bool(basic_eqvp(x, y) || (!C_immediatep(x) && !C_immediatep(y) && (C_block_header(x) == C_block_header(y) && C_header_bits(x) == C_STRUCTURE_TYPE && C_block_item(x, 0) == C_block_item(y, 0) && (C_block_item(x, 0) == ratnum_type_tag || C_block_item(x, 0) == compnum_type_tag) && basic_eqvp(C_block_item(x, 1), C_block_item(y, 1)) && basic_eqvp(C_block_item(x, 2), C_block_item(y, 2))))); } C_regparm C_word C_fcall C_i_nanp(C_word x) { if (x & C_FIXNUM_BIT) { return C_SCHEME_FALSE; } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x); } else if (C_block_header(x) == C_FLONUM_TAG) { return C_u_i_flonum_nanp(x); } else if (C_IS_BIGNUM_TYPE(x)) { return C_SCHEME_FALSE; } else if (C_header_bits(x) == C_STRUCTURE_TYPE) { /* To make this inlineable we don't call try_extended_number */ if (C_block_item(x, 0) == ratnum_type_tag) return C_SCHEME_FALSE; else if (C_block_item(x, 0) == compnum_type_tag) return C_mk_bool(C_truep(C_i_nanp(C_block_item(x, 1))) || C_truep(C_i_nanp(C_block_item(x, 2)))); else barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x); } } /* TODO: Rename to C_i_finitep */ C_regparm C_word C_fcall C_i_numbers_finitep(C_word x) { if (x & C_FIXNUM_BIT) { return C_SCHEME_TRUE; } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x); } else if (C_block_header(x) == C_FLONUM_TAG) { return C_u_i_flonum_finitep(x); } else if (C_IS_BIGNUM_TYPE(x)) { return C_SCHEME_TRUE; } else if (C_header_bits(x) == C_STRUCTURE_TYPE) { /* To make this inlineable we don't use try_extended_number */ if (C_block_item(x, 0) == ratnum_type_tag) return C_SCHEME_TRUE; else if (C_block_item(x, 0) == compnum_type_tag) return C_and(C_i_numbers_finitep(C_block_item(x, 1)), C_i_numbers_finitep(C_block_item(x, 2))); else barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x); } } /* TODO: Rename to C_i_infinitep */ C_regparm C_word C_fcall C_i_numbers_infinitep(C_word x) { if (x & C_FIXNUM_BIT) { return C_SCHEME_FALSE; } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x); } else if (C_block_header(x) == C_FLONUM_TAG) { return C_u_i_flonum_infinitep(x); } else if (C_IS_BIGNUM_TYPE(x)) { return C_SCHEME_FALSE; } else if (C_header_bits(x) == C_STRUCTURE_TYPE) { /* To make this inlineable we don't use try_extended_number */ if (C_block_item(x, 0) == ratnum_type_tag) return C_SCHEME_FALSE; else if (C_block_item(x, 0) == compnum_type_tag) return C_mk_bool(C_truep(C_i_numbers_infinitep(C_block_item(x, 1))) || C_truep(C_i_numbers_infinitep(C_block_item(x, 2)))); else barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x); } } /* TODO: Rename to C_i_zerop */ C_regparm C_word C_fcall C_i_numbers_zerop(C_word x) { if (x & C_FIXNUM_BIT) { return C_mk_bool(x == C_fix(0)); } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "zero?", x); } else if (C_block_header(x) == C_FLONUM_TAG) { return C_mk_bool(C_flonum_magnitude(x) == 0.0); } else if (C_IS_BIGNUM_TYPE(x) || (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == ratnum_type_tag || C_block_item(x, 0) == compnum_type_tag))) { return C_SCHEME_FALSE; } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "zero?", x); } } /* Copy all the digits from source to target, obliterating what was * there. If target is larger than source, the most significant * digits will remain untouched. */ C_inline void bignum_digits_destructive_copy(C_word target, C_word source) { C_memcpy(C_bignum_digits(target), C_bignum_digits(source), C_wordstobytes(C_bignum_size(source))); } void C_ccall C_2_basic_plus(C_word c, C_word self, C_word k, C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, (double)C_unfix(x) + C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_plus(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-plus", 3, k, x, y); } } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); if (y & C_FIXNUM_BIT) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) + (double)C_unfix(y))); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_kontinue(k, C_a_i_flonum_plus(&a, 2, x, y)); } else if (C_IS_BIGNUM_TYPE(y)) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)+C_bignum_to_double(y))); } else { try_extended_number("numbers#@extended-2-plus", 3, k, x, y); } } else if (C_IS_BIGNUM_TYPE(x)) { if (y & C_FIXNUM_BIT) { C_u_2_integer_plus(4, (C_word)NULL, k, x, y); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)+C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_plus(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-plus", 3, k, x, y); } } else { try_extended_number("numbers#@extended-2-plus", 3, k, x, y); } } void C_ccall C_u_2_integer_plus(C_word c, C_word self, C_word k, C_word x, C_word y) { if ((x & y) & C_FIXNUM_BIT) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y)); } else { C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab; if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&a, x); if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&a, y); if (C_bignum_negativep(x)) { if (C_bignum_negativep(y)) bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE); else bignum_minus_unsigned(k, y, x); } else { if (C_bignum_negativep(y)) bignum_minus_unsigned(k, x, y); else bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE); } } } /* Needs C_SIZEOF_FIX_BIGNUM */ C_regparm C_word C_fcall C_a_u_i_2_fixnum_plus(C_word **ptr, C_word n, C_word x, C_word y) { /* Exceptional situation: this will cause a real underflow */ if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) { return C_bignum1(ptr, 1, ((C_uword)-C_MOST_NEGATIVE_FIXNUM) << 1); } else { C_word z = C_unfix(x) + C_unfix(y); if(!C_fitsinfixnump(z)) { /* TODO: function/macro returning either fixnum or bignum from a C int */ /* This should help with the C API/FFI too. */ return C_bignum1(ptr, z < 0, labs(z)); } else { return C_fix(z); } } } void C_ccall C_2_basic_minus(C_word c, C_word self, C_word k, C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, (double)C_unfix(x) - C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_minus(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-minus", 3, k, x, y); } } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); if (y & C_FIXNUM_BIT) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) - (double)C_unfix(y))); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_kontinue(k, C_a_i_flonum_difference(&a, 2, x, y)); /* XXX NAMING! */ } else if (C_IS_BIGNUM_TYPE(y)) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)-C_bignum_to_double(y))); } else { try_extended_number("numbers#@extended-2-minus", 3, k, x, y); } } else if (C_IS_BIGNUM_TYPE(x)) { if (y & C_FIXNUM_BIT) { C_u_2_integer_minus(4, (C_word)NULL, k, x, y); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)-C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_minus(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-minus", 3, k, x, y); } } else { try_extended_number("numbers#@extended-2-minus", 3, k, x, y); } } void C_ccall C_u_2_integer_minus(C_word c, C_word self, C_word k, C_word x, C_word y) { if ((x & y) & C_FIXNUM_BIT) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y)); } else { C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab; if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&a, x); if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&a, y); if (C_bignum_negativep(x)) { if (C_bignum_negativep(y)) bignum_minus_unsigned(k, y, x); else bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE); } else { if (C_bignum_negativep(y)) bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE); else bignum_minus_unsigned(k, x, y); } } } /* Needs C_SIZEOF_FIX_BIGNUM */ C_regparm C_word C_fcall C_a_u_i_2_fixnum_minus(C_word **ptr, C_word n, C_word x, C_word y) { C_word z = C_unfix(x) - C_unfix(y); if(!C_fitsinfixnump(z)) { /* TODO: function/macro returning either fixnum or bignum from a C int */ /* This should help with the C API/FFI too. */ return C_bignum1(ptr, z < 0, labs(z)); } else { return C_fix(z); } } static void bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp) { C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size; if (C_bignum_size(y) > C_bignum_size(x)) { /* Ensure size(y) <= size(x) */ C_word z = x; x = y; y = z; } k2 = C_closure(&ka, 4, (C_word)bignum_plus_unsigned_2, k, x, y); size = C_fix(C_bignum_size(x) + 1); /* One more digit, for possible carry. */ C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } static void bignum_plus_unsigned_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3); C_uword *scan_y = C_bignum_digits(y), *end_y = scan_y + C_bignum_size(y), *scan_r = C_bignum_digits(result), *end_r = scan_r + C_bignum_size(result), sum, digit; int carry = 0; /* Copy x into r so we can operate on two pointers, which is faster * than three, and we can stop earlier after adding y. It's slower * if x and y have equal length. On average it's slightly faster. */ bignum_digits_destructive_copy(result, x); *(end_r-1) = 0; /* Ensure most significant digit is initialised */ /* Move over x and y simultaneously, destructively adding digits w/ carry. */ while (scan_y < end_y) { digit = *scan_r; if (carry) { sum = digit + *scan_y++ + 1; carry = sum <= digit; } else { sum = digit + *scan_y++; carry = sum < digit; } (*scan_r++) = sum; } /* The end of y, the smaller number. Propagate carry into the rest of x. */ while (carry) { sum = (*scan_r) + 1; carry = (sum == 0); (*scan_r++) = sum; } assert(scan_r <= end_r); C_kontinue(k, C_bignum_simplify(result)); } static int bignum_cmp_unsigned(C_word x, C_word y) { C_word xlen = C_bignum_size(x), ylen = C_bignum_size(y); if (xlen < ylen) { return -1; } else if (xlen > ylen) { return 1; } else if (x == y) { return 0; } else { C_uword *startx = C_bignum_digits(x), *scanx = startx + xlen, *scany = C_bignum_digits(y) + ylen; while (startx < scanx) { C_uword xdigit = (*--scanx), ydigit = (*--scany); if (xdigit < ydigit) return -1; if (xdigit > ydigit) return 1; } return 0; } } static void bignum_minus_unsigned(C_word k, C_word x, C_word y) { C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size; switch(bignum_cmp_unsigned(x, y)) { case 0: /* x = y, return 0 */ C_kontinue(k, C_fix(0)); case -1: /* abs(x) < abs(y), return -(abs(y) - abs(x)) */ k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, y, x); size = C_fix(C_bignum_size(y)); /* Maximum size of result is length of y. */ C_allocate_bignum(5, (C_word)NULL, k2, size, C_SCHEME_TRUE, C_SCHEME_FALSE); case 1: /* abs(x) > abs(y), return abs(x) - abs(y) */ default: k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, x, y); size = C_fix(C_bignum_size(x)); /* Maximum size of result is length of x. */ C_allocate_bignum(5, (C_word)NULL, k2, size, C_SCHEME_FALSE, C_SCHEME_FALSE); break; } } static void bignum_minus_unsigned_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3); C_uword *scan_r = C_bignum_digits(result), *end_r = scan_r + C_bignum_size(result), *scan_y = C_bignum_digits(y), *end_y = scan_y + C_bignum_size(y), difference, digit; int borrow = 0; bignum_digits_destructive_copy(result, x); /* See bignum_plus_unsigned_2 */ /* Destructively subtract y's digits w/ borrow from and back into r. */ while (scan_y < end_y) { digit = *scan_r; if (borrow) { difference = digit - *scan_y++ - 1; borrow = difference >= digit; } else { difference = digit - *scan_y++; borrow = difference > digit; } (*scan_r++) = difference; } /* The end of y, the smaller number. Propagate borrow into the rest of x. */ while (borrow) { digit = *scan_r; difference = digit - borrow; borrow = difference >= digit; (*scan_r++) = difference; } assert(scan_r <= end_r); C_kontinue(k, C_bignum_simplify(result)); } C_word C_ccall C_u_i_2_fixnum_gcd(C_word x, C_word y) { x = (x & C_INT_SIGN_BIT) ? -C_unfix(x) : C_unfix(x); y = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y); while(y != 0) { C_word r = x % y; x = y; y = r; } return C_fix(x); } C_word C_ccall C_a_u_i_2_flonum_gcd(C_word **p, C_word n, C_word x, C_word y) { double xub, yub, r; if (!C_truep(C_u_i_fpintegerp_fixed(x))) barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "gcd", x); if (!C_truep(C_u_i_fpintegerp_fixed(y))) barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "gcd", y); xub = C_flonum_magnitude(x); yub = C_flonum_magnitude(y); if (xub < 0.0) xub = -xub; if (yub < 0.0) yub = -yub; while(yub != 0.0) { r = fmod(xub, yub); xub = yub; yub = r; } return C_flonum(p, xub); } void C_ccall C_basic_abs(C_word c, C_word self, C_word k, C_word x) { if (c != 3) { C_bad_argc_2(c, 3, self); } else if (x & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x)); } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "abs", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_a_i_flonum_abs(&a, 1, x)); } else if (C_IS_BIGNUM_TYPE(x)) { C_u_integer_abs(3, (C_word)NULL, k, x); } else { try_extended_number("numbers#@extended-abs", 2, k, x); } } void C_ccall C_u_integer_abs(C_word c, C_word self, C_word k, C_word x) { if (x & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x)); } else if (C_bignum_negativep(x)) { C_u_integer_negate(3, (C_word)NULL, k, x); } else { C_kontinue(k, x); } } void C_ccall C_basic_signum(C_word c, C_word self, C_word k, C_word x) { if (c != 3) { C_bad_argc_2(c, 3, self); } else if (x & C_FIXNUM_BIT) { C_kontinue(k, C_u_i_fixnum_signum(x)); } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "signum", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_a_u_i_flonum_signum(&a, 1, x)); } else if (C_IS_BIGNUM_TYPE(x)) { C_kontinue(k, C_bignum_negativep(x) ? C_fix(-1) : C_fix(1)); } else { try_extended_number("numbers#@extended-signum", 2, k, x); } } C_regparm C_word C_fcall C_u_i_integer_signum(C_word x) { if (x & C_FIXNUM_BIT) return C_u_i_fixnum_signum(x); else return (C_bignum_negativep(x) ? C_fix(-1) : C_fix(1)); } C_regparm C_word C_fcall C_i_basic_evenp(C_word x) { if(x & C_FIXNUM_BIT) { return C_i_fixnumevenp(x); } else if(C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x); } else if (C_block_header(x) == C_FLONUM_TAG) { double val, dummy; val = C_flonum_magnitude(x); if(C_isnan(val) || C_isinf(val) || C_modf(val, &dummy) != 0.0) barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x); else return C_mk_bool(fmod(val, 2.0) == 0.0); } else if (C_IS_BIGNUM_TYPE(x)) { return C_mk_nbool(C_bignum_digits(x)[0] & 1); } else { /* No need to try extended number */ barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x); } } C_regparm C_word C_fcall C_u_i_integer_evenp(C_word x) { if (x & C_FIXNUM_BIT) return C_i_fixnumevenp(x); return C_mk_nbool(C_bignum_digits(x)[0] & 1); } C_regparm C_word C_fcall C_i_basic_oddp(C_word x) { if(x & C_FIXNUM_BIT) { return C_i_fixnumoddp(x); } else if(C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x); } else if(C_block_header(x) == C_FLONUM_TAG) { double val, dummy; val = C_flonum_magnitude(x); if(C_isnan(val) || C_isinf(val) || C_modf(val, &dummy) != 0.0) barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x); else return C_mk_bool(fmod(val, 2.0) != 0.0); } else if (C_IS_BIGNUM_TYPE(x)) { return C_mk_bool(C_bignum_digits(x)[0] & 1); } else { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x); } } C_regparm C_word C_fcall C_u_i_integer_oddp(C_word x) { if (x & C_FIXNUM_BIT) return C_i_fixnumoddp(x); return C_mk_bool(C_bignum_digits(x)[0] & 1); } void C_ccall C_basic_negate(C_word c, C_word self, C_word k, C_word x) { if (x & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x)); } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_a_i_flonum_negate(&a, 1, x)); } else if (C_IS_BIGNUM_TYPE(x)) { C_u_integer_negate(3, (C_word)NULL, k, x); } else { try_extended_number("numbers#@extended-negate", 2, k, x); } } void C_ccall C_u_integer_negate(C_word c, C_word self, C_word k, C_word x) { if (x & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x)); } else { if (C_bignum_negated_fitsinfixnump(x)) { C_kontinue(k, C_fix(C_MOST_NEGATIVE_FIXNUM)); } else { C_word *ka, k2, negp = C_mk_nbool(C_bignum_negativep(x)), size = C_fix(C_bignum_size(x)); ka = C_alloc(C_SIZEOF_CLOSURE(3)); k2 = C_closure(&ka, 3, (C_word)bignum_negate_2, k, x); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } } static void bignum_negate_2(C_word c, C_word self, C_word new_big) { C_word k = C_block_item(self, 1), old_big = C_block_item(self, 2); bignum_digits_destructive_copy(new_big, old_big); C_kontinue(k, C_bignum_simplify(new_big)); } C_regparm C_word C_fcall C_a_u_i_fixnum_negate(C_word **ptr, C_word n, C_word x) { /* Exceptional situation: this will cause an overflow to itself */ if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) /* C_fitsinfixnump(x) */ return C_bignum1(ptr, 0, -C_MOST_NEGATIVE_FIXNUM); else return C_fix(-C_unfix(x)); } /* TODO: Rename to C_nequalp */ void C_ccall C_numbers_nequalp(C_word c, C_word self, C_word k, ...) { C_word x, y, result; va_list v; if (c < 4) C_bad_argc_2(c, 4, self); c -= 2; va_start(v, k); x = va_arg(v, C_word); while(--c) { y = va_arg(v, C_word); result = C_i_2_basic_equalp(x, y); if (result == C_SCHEME_FALSE) break; } va_end(v); C_kontinue(k, result); } /* Compare two numbers as ratnums. Either may be rat-, fix- or bignums */ static C_word rat_cmp(C_word x, C_word y) { C_word ab[C_SIZEOF_FIX_BIGNUM*4], *a = ab, x1, x2, y1, y2, s, t, ssize, tsize, result, negp; C_uword *scan; /* Check for 1 or 0; if x or y is this, the other must be the ratnum */ if (x == C_fix(0)) { /* Only the sign of y1 matters */ return basic_cmp(x, C_block_item(y, 1), "ratcmp", 0); } else if (x == C_fix(1)) { /* x1*y1 <> x2*y2 --> y2 <> y1 | x1/x2 = 1/1 */ return basic_cmp(C_block_item(y, 2), C_block_item(y, 1), "ratcmp", 0); } else if (y == C_fix(0)) { /* Only the sign of x1 matters */ return basic_cmp(C_block_item(x, 1), y, "ratcmp", 0); } else if (y == C_fix(1)) { /* x1*y1 <> x2*y2 --> x1 <> x2 | y1/y2 = 1/1 */ return basic_cmp(C_block_item(x, 1), C_block_item(x, 2), "ratcmp", 0); } /* Extract components x=x1/x2 and y=y1/y2 */ if (x & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(x)) x1 = x, x2 = C_fix(1); else x1 = C_block_item(x, 1), x2 = C_block_item(x, 2); if (y & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(y)) y1 = y, y2 = C_fix(1); else y1 = C_block_item(y, 1), y2 = C_block_item(y, 2); /* We only want to deal with bignums (this is tricky enough) */ if (x1 & C_FIXNUM_BIT) x1 = C_a_u_i_fix_to_big(&a, x1); if (x2 & C_FIXNUM_BIT) x2 = C_a_u_i_fix_to_big(&a, x2); if (y1 & C_FIXNUM_BIT) y1 = C_a_u_i_fix_to_big(&a, y1); if (y2 & C_FIXNUM_BIT) y2 = C_a_u_i_fix_to_big(&a, y2); /* We multiply using schoolbook method, so this will be very slow in * extreme cases. This is a tradeoff we make so that comparisons * are inlineable, which makes a big difference for the common case. */ ssize = C_bignum_size(x1) + C_bignum_size(y2); negp = C_mk_bool(C_bignum_negativep(x1)); s = allocate_tmp_bignum(C_fix(ssize), negp, C_SCHEME_TRUE); bignum_digits_multiply(x1, y2, s); /* Swap args if x1 < y2? */ tsize = C_bignum_size(y1) + C_bignum_size(x2); negp = C_mk_bool(C_bignum_negativep(y1)); t = allocate_tmp_bignum(C_fix(tsize), negp, C_SCHEME_TRUE); bignum_digits_multiply(y1, x2, t); /* Swap args if y1 < x2? */ /* Shorten the numbers if needed */ for (scan = C_bignum_digits(s)+ssize-1; *scan == 0; scan--) ssize--; C_bignum_mutate_size(s, ssize); for (scan = C_bignum_digits(t)+tsize-1; *scan == 0; scan--) tsize--; C_bignum_mutate_size(t, tsize); result = C_u_i_bignum_cmp(s, t); free_tmp_bignum(t); free_tmp_bignum(s); return result; } /* The primitive comparison operator. eqp should be 1 if we're only * interested in equality testing (can speed things up and in case of * compnums, equality checking is the only available operation). This * may return #f, in case there is no answer (for NaNs) or as a quick * and dirty non-zero answer when eqp is true. Ugly but effective :) */ static C_word basic_cmp(C_word x, C_word y, char *loc, int eqp) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { return C_fix((x < y) ? -1 : ((x > y) ? 1 : 0)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } else if (C_block_header(y) == C_FLONUM_TAG) { return int_flo_cmp(x, y); } else if (C_IS_BIGNUM_TYPE(y)) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; return C_u_i_bignum_cmp(C_a_u_i_fix_to_big(&a, x), y); } else if (C_header_bits(y) == C_STRUCTURE_TYPE && C_block_item(y, 0) == ratnum_type_tag) { if (eqp) return C_SCHEME_FALSE; else return rat_cmp(x, y); } else if (C_block_item(y, 0) == compnum_type_tag) { if (eqp) return C_SCHEME_FALSE; else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x); } else if (C_block_header(x) == C_FLONUM_TAG) { if (y & C_FIXNUM_BIT) { return flo_int_cmp(x, y); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } else if (C_block_header(y) == C_FLONUM_TAG) { double a = C_flonum_magnitude(x), b = C_flonum_magnitude(y); if (C_isnan(a) || C_isnan(b)) return C_SCHEME_FALSE; /* "mu" */ else return C_fix((a < b) ? -1 : ((a > b) ? 1 : 0)); } else if (C_IS_BIGNUM_TYPE(y)) { return flo_int_cmp(x, y); } else if (C_header_bits(y) == C_STRUCTURE_TYPE) { if (C_block_item(y, 0) == ratnum_type_tag) { return flo_rat_cmp(x, y); } else if (C_block_item(y, 0) == compnum_type_tag) { if (eqp) return C_SCHEME_FALSE; else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else if (C_IS_BIGNUM_TYPE(x)) { if (y & C_FIXNUM_BIT) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; return C_u_i_bignum_cmp(x, C_a_u_i_fix_to_big(&a, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } else if (C_block_header(y) == C_FLONUM_TAG) { return int_flo_cmp(x, y); } else if (C_IS_BIGNUM_TYPE(y)) { return C_u_i_bignum_cmp(x, y); } else if (C_header_bits(y) == C_STRUCTURE_TYPE && C_block_item(y, 0) == ratnum_type_tag) { if (eqp) return C_SCHEME_FALSE; else return rat_cmp(x, y); } else if (C_block_item(y, 0) == compnum_type_tag) { if (eqp) return C_SCHEME_FALSE; else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == ratnum_type_tag)) { if (y & C_FIXNUM_BIT) { if (eqp) return C_SCHEME_FALSE; else return rat_cmp(x, y); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } else if (C_block_header(y) == C_FLONUM_TAG) { return rat_flo_cmp(x, y); } else if (C_IS_BIGNUM_TYPE(y)) { if (eqp) return C_SCHEME_FALSE; else return rat_cmp(x, y); } else if (C_header_bits(y) == C_STRUCTURE_TYPE && (C_block_item(y, 0) == ratnum_type_tag)) { if (eqp) { return C_and(C_and(C_u_i_2_integer_equalp(C_block_item(x, 1), C_block_item(y, 1)), C_u_i_2_integer_equalp(C_block_item(x, 2), C_block_item(y, 2))), C_fix(0)); } else { return rat_cmp(x, y); } } else if (C_header_bits(y) == C_STRUCTURE_TYPE && (C_block_item(y, 0) == compnum_type_tag)) { if (eqp) return C_SCHEME_FALSE; else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else if (C_header_bits(x) == C_STRUCTURE_TYPE && (C_block_item(x, 0) == compnum_type_tag)) { if (!eqp) { barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, x); } else if (y & C_FIXNUM_BIT) { return C_SCHEME_FALSE; } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } else if (C_block_header(y) == C_FLONUM_TAG || C_IS_BIGNUM_TYPE(y) || (C_header_bits(y) == C_STRUCTURE_TYPE && (C_block_item(y, 0) == ratnum_type_tag))) { return C_SCHEME_FALSE; } else if (C_header_bits(y) == C_STRUCTURE_TYPE && (C_block_item(y, 0) == compnum_type_tag)) { return C_and(C_and(C_i_2_basic_equalp(C_block_item(x, 1), C_block_item(y, 1)), C_i_2_basic_equalp(C_block_item(x, 2), C_block_item(y, 2))), C_fix(0)); } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y); } } else { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x); } } C_regparm C_word C_fcall C_i_2_basic_equalp(C_word x, C_word y) { return C_mk_bool(basic_cmp(x, y, "=", 1) == C_fix(0)); } C_word C_ccall C_u_i_2_integer_equalp(C_word x, C_word y) { if (x & C_FIXNUM_BIT) return C_mk_bool(x == y); else if (y & C_FIXNUM_BIT) return C_SCHEME_FALSE; else return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(0)); } /* TODO: Rename to C_lessp */ void C_ccall C_numbers_lessp(C_word c, C_word self, C_word k, ...) { C_word x, y, result; va_list v; if (c < 4) C_bad_argc_2(c, 4, self); c -= 2; va_start(v, k); x = va_arg(v, C_word); while(--c) { y = va_arg(v, C_word); result = C_i_2_basic_lessp(x, y); if (result == C_SCHEME_FALSE) break; x = y; } va_end(v); C_kontinue(k, result); } /* TODO: Rename to C_less_or_equal_p */ void C_ccall C_numbers_less_or_equal_p(C_word c, C_word self, C_word k, ...) { C_word x, y, result; va_list v; if (c < 4) C_bad_argc_2(c, 4, self); c -= 2; va_start(v, k); x = va_arg(v, C_word); while(--c) { y = va_arg(v, C_word); result = C_i_2_basic_less_or_equalp(x, y); if (result == C_SCHEME_FALSE) break; x = y; } va_end(v); C_kontinue(k, result); } C_regparm C_word C_fcall C_i_2_basic_lessp(C_word x, C_word y) { return C_mk_bool(basic_cmp(x, y, "<", 0) == C_fix(-1)); } C_regparm C_word C_fcall C_i_2_basic_less_or_equalp(C_word x, C_word y) { C_word res = basic_cmp(x, y, "<=", 0); return C_mk_bool(res == C_fix(0) || res == C_fix(-1)); } C_word C_ccall C_u_i_2_integer_lessp(C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { return C_mk_bool(C_unfix(x) < C_unfix(y)); } else { return C_mk_nbool(C_bignum_negativep(y)); } } else if (y & C_FIXNUM_BIT) { return C_mk_bool(C_bignum_negativep(x)); } else { return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(-1)); } } C_word C_ccall C_u_i_2_integer_less_or_equalp(C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { return C_mk_bool(C_unfix(x) <= C_unfix(y)); } else { return C_mk_nbool(C_bignum_negativep(y)); } } else if (y & C_FIXNUM_BIT) { return C_mk_bool(C_bignum_negativep(x)); } else { C_word res = C_u_i_bignum_cmp(x, y); return C_mk_bool(res == C_fix(0) || res == C_fix(-1)); } } /* TODO: Rename to C_greater_p */ void C_ccall C_numbers_greaterp(C_word c, C_word self, C_word k, ...) { C_word x, y, result; va_list v; if (c < 4) C_bad_argc_2(c, 4, self); c -= 2; va_start(v, k); x = va_arg(v, C_word); while(--c) { y = va_arg(v, C_word); result = C_i_2_basic_greaterp(x, y); if (result == C_SCHEME_FALSE) break; x = y; } va_end(v); C_kontinue(k, result); } /* TODO: Rename to C_greater_or_equal_p */ void C_ccall C_numbers_greater_or_equal_p(C_word c, C_word self, C_word k, ...) { C_word x, y, result; va_list v; if (c < 4) C_bad_argc_2(c, 4, self); c -= 2; va_start(v, k); x = va_arg(v, C_word); while(--c) { y = va_arg(v, C_word); result = C_i_2_basic_greater_or_equalp(x, y); if (result == C_SCHEME_FALSE) break; x = y; } va_end(v); C_kontinue(k, result); } C_regparm C_word C_fcall C_i_2_basic_greaterp(C_word x, C_word y) { return C_mk_bool(basic_cmp(x, y, ">", 0) == C_fix(1)); } C_regparm C_word C_fcall C_i_2_basic_greater_or_equalp(C_word x, C_word y) { C_word res = basic_cmp(x, y, ">=", 0); return C_mk_bool(res == C_fix(0) || res == C_fix(1)); } C_word C_ccall C_u_i_2_integer_greaterp(C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { return C_mk_bool(C_unfix(x) > C_unfix(y)); } else { return C_mk_bool(C_bignum_negativep(y)); } } else if (y & C_FIXNUM_BIT) { return C_mk_nbool(C_bignum_negativep(x)); } else { return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(1)); } } C_word C_ccall C_u_i_2_integer_greater_or_equalp(C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { return C_mk_bool(C_unfix(x) >= C_unfix(y)); } else { return C_mk_bool(C_bignum_negativep(y)); } } else if (y & C_FIXNUM_BIT) { return C_mk_nbool(C_bignum_negativep(x)); } else { C_word res = C_u_i_bignum_cmp(x, y); return C_mk_bool(res == C_fix(0) || res == C_fix(1)); } } /* This is a bit weird: We have to compare flonums as bignums due to * precision loss on 64-bit platforms. For simplicity, we convert * fixnums to bignums here. */ static C_word int_flo_cmp(C_word intnum, C_word flonum) { C_word ab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_FLONUM], *a = ab, x, y, res; double i, f; f = C_flonum_magnitude(flonum); if (C_isnan(f)) { return C_SCHEME_FALSE; /* "mu" */ } else if (C_isinf(f)) { return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */ } else { f = modf(f, &i); x = (intnum & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, intnum) : intnum; y = flo_to_tmp_bignum(C_flonum(&a, i)); res = C_u_i_bignum_cmp(x, y); free_tmp_bignum(y); if (res == C_fix(0)) /* Use fraction to break tie. If f > 0, x is smaller */ return C_fix((f > 0.0) ? -1 : ((f < 0.0) ? 1 : 0)); else return res; } } /* For convenience (ie, to reduce the degree of mindfuck) */ static C_word flo_int_cmp(C_word flonum, C_word intnum) { C_word res = int_flo_cmp(intnum, flonum); switch(res) { case C_fix(1): return C_fix(-1); case C_fix(-1): return C_fix(1); default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */ } } /* This code is completely braindead, but at least it allows us to do * inline comparisons! */ static C_word rat_flo_cmp(C_word ratnum, C_word flonum) { C_word ab[C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_FLONUM], *a = ab, num, denom, ibig, res, nscaled, iscaled, negp; C_uword *scan; int shift_amount, ilen, nlen; double i, f; f = C_flonum_magnitude(flonum); if (C_isnan(f)) { return C_SCHEME_FALSE; /* "mu" */ } else if (C_isinf(f)) { return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */ } else { /* Scale up the floating-point number to become a whole integer, * and remember power of two (# of bits) to shift the numerator. */ shift_amount = 0; /* TODO: This doesn't work for denormalized flonums! */ while (modf(f, &i) != 0.0) { f = ldexp(f, 1); shift_amount++; } i = f; /* TODO: split i and f so it'll work for denormalized flonums */ num = C_block_item(ratnum, 1); num = (num & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, num) : num; if (C_bignum_negativep(num) && i >= 0.0) { /* Save time if signs differ */ return C_fix(-1); } else if (!C_bignum_negativep(num) && i <= 0.0) { /* num is never 0 */ return C_fix(1); } else { negp = C_mk_bool(C_bignum_negativep(num)); denom = C_block_item(ratnum, 2); denom = (denom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, denom) : denom; ibig = flo_to_tmp_bignum(C_flonum(&a, i)); nlen = C_bignum_size(num) + C_bignum_size(denom); ilen = C_bignum_size(ibig) + C_bignum_size(denom); /* Now, multiply the scaled flonum by the denominator, so we can * compare it directly to the scaled numerator. Unfortunately, * this won't use Karatsuba multiplication, so for large numbers * it will be slower than it could be if comparisons were done * in CPS context. */ iscaled = allocate_tmp_bignum(C_fix(ilen), negp, C_SCHEME_TRUE); bignum_digits_multiply(denom, ibig, iscaled); /* Swap args if i < d? */ free_tmp_bignum(ibig); nlen += C_BIGNUM_BITS_TO_DIGITS(shift_amount); nscaled = allocate_tmp_bignum(C_fix(nlen), negp, C_SCHEME_TRUE); scan = C_bignum_digits(nscaled) + shift_amount / C_BIGNUM_DIGIT_LENGTH; C_memcpy(scan, C_bignum_digits(num), C_wordstobytes(C_bignum_size(num))); shift_amount = shift_amount % C_BIGNUM_DIGIT_LENGTH; if(shift_amount > 0) { bignum_digits_destructive_shift_left( scan, C_bignum_digits(nscaled) + nlen, shift_amount); } /* Shorten the numbers if needed */ for (scan = C_bignum_digits(iscaled)+ilen-1; *scan == 0; scan--) ilen--; C_bignum_mutate_size(iscaled, ilen); for (scan = C_bignum_digits(nscaled)+nlen-1; *scan == 0; scan--) nlen--; C_bignum_mutate_size(nscaled, nlen); /* Finally, we're ready to compare them! */ res = C_u_i_bignum_cmp(nscaled, iscaled); free_tmp_bignum(nscaled); free_tmp_bignum(iscaled); return res; } } } static C_word flo_rat_cmp(C_word flonum, C_word ratnum) { C_word res = rat_flo_cmp(ratnum, flonum); switch(res) { case C_fix(1): return C_fix(-1); case C_fix(-1): return C_fix(1); default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */ } } C_word C_u_i_bignum_cmp(C_word x, C_word y) { if (C_bignum_negativep(x)) { if (C_bignum_negativep(y)) { /* Largest negative number is smallest */ return C_fix(bignum_cmp_unsigned(y, x)); } else { return C_fix(-1); } } else { if (C_bignum_negativep(y)) { return C_fix(1); } else { return C_fix(bignum_cmp_unsigned(x, y)); } } } void C_ccall C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp) { C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init; k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp); init = C_and(initp, C_make_character('\0')); C_allocate_vector(6, (C_word)NULL, k2, C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */ /* Byte vec, initialization, align at 8 bytes */ C_SCHEME_TRUE, init, C_SCHEME_FALSE); } static void allocate_bignum_2(C_word c, C_word self, C_word bigvec) { C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum, k = C_block_item(self, 1), negp = C_truep(C_block_item(self, 2)) ? 1 : 0, tagvec = CHICKEN_gc_root_ref(tags); C_set_block_item(bigvec, 0, negp); bignum = C_a_i_record2(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec); C_kontinue(k, bignum); } static C_word allocate_tmp_bignum(C_word size, C_word negp, C_word initp) { C_word *mem = malloc(C_wordstobytes(C_SIZEOF_BIGNUM(C_unfix(size)))), bigvec = (C_word)(mem + C_SIZEOF_STRUCTURE(2)), tagvec = CHICKEN_gc_root_ref(tags); if (mem == NULL) abort(); /* TODO: panic */ C_block_header_init(bigvec, (C_STRING_TYPE | C_wordstobytes(C_unfix(size)+1))); C_set_block_item(bigvec, 0, C_truep(negp)); if (C_truep(initp)) { C_memset(((C_uword *)C_data_pointer(bigvec))+1, 0, C_wordstobytes(C_unfix(size))); } return C_a_i_record2(&mem, 2, C_block_item(tagvec, BIG_TAG), bigvec); } /* Simplification: scan trailing zeroes, then return a fixnum if the * value fits, or trim the bignum's length. */ C_word C_ccall C_bignum_simplify(C_word big) { C_uword *start = C_bignum_digits(big), *last_digit = start + C_bignum_size(big) - 1, *scan = last_digit, tmp; int length; while (scan >= start && *scan == 0) scan--; length = scan - start + 1; switch(length) { case 0: return C_fix(0); case 1: tmp = *start; if (C_bignum_negativep(big) ? !(tmp & C_INT_SIGN_BIT) && C_fitsinfixnump(-(C_word)tmp) : C_ufitsinfixnump(tmp)) return C_bignum_negativep(big) ? C_fix(-(C_word)tmp) : C_fix(tmp); /* FALLTHROUGH */ default: if (scan < last_digit) C_bignum_mutate_size(big, length); return big; } } static C_uword bignum_digits_destructive_scale_up_with_carry(C_uword *start, C_uword *end, C_uword factor, C_uword carry) { C_uword digit, p; assert(C_fitsinbignumhalfdigitp(carry)); assert(C_fitsinbignumhalfdigitp(factor)); /* See fixnum_times. Substitute xlo = factor, xhi = 0, y = digit * and simplify the result to reduce variable usage. */ while (start < end) { digit = (*start); p = factor * C_BIGNUM_DIGIT_LO_HALF(digit) + carry; carry = C_BIGNUM_DIGIT_LO_HALF(p); p = factor * C_BIGNUM_DIGIT_HI_HALF(digit) + C_BIGNUM_DIGIT_HI_HALF(p); (*start++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), carry); carry = C_BIGNUM_DIGIT_HI_HALF(p); } return carry; } static C_uword bignum_digits_destructive_scale_down(C_uword *start, C_uword *end, C_uword denominator) { C_uword digit, k = 0; C_uhword q_j_hi, q_j_lo; /* Single digit divisor case from Hacker's Delight, Figure 9-1, * adapted to modify u[] in-place instead of writing to q[]. */ while (start < end) { digit = (*--end); k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_HI_HALF(digit)); /* j */ q_j_hi = k / denominator; k -= q_j_hi * denominator; k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_LO_HALF(digit)); /* j-1 */ q_j_lo = k / denominator; k -= q_j_lo * denominator; *end = C_BIGNUM_DIGIT_COMBINE(q_j_hi, q_j_lo); } return k; } static C_uword bignum_digits_destructive_shift_right(C_uword *start, C_uword *end, int shift_right, int negp) { int shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right; C_uword digit, carry = negp ? ((~(C_uword)0) << shift_left) : 0; assert(shift_right < C_BIGNUM_DIGIT_LENGTH); while (start < end) { digit = *(--end); *end = (digit >> shift_right) | carry; carry = digit << shift_left; } return carry >> shift_left; /* The bits that were shifted out to the right */ } static C_uword bignum_digits_destructive_shift_left(C_uword *start, C_uword *end, int shift_left) { C_uword carry = 0, digit; int shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left; assert(shift_left < C_BIGNUM_DIGIT_LENGTH); while (start < end) { digit = *start; (*start++) = (digit << shift_left) | carry; carry = digit >> shift_right; } return carry; /* This would end up as most significant digit if it fit */ } void C_ccall C_2_basic_times(C_word c, C_word self, C_word k, C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_BIGNUM(2)); C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, (double)C_unfix(x) * C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_times(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-times", 3, k, x, y); } } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x); } else if (C_block_header(x) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); if (y & C_FIXNUM_BIT) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) * (double)C_unfix(y))); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_kontinue(k, C_a_i_flonum_times(&a, 2, x, y)); } else if (C_IS_BIGNUM_TYPE(y)) { C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)*C_bignum_to_double(y))); } else { try_extended_number("numbers#@extended-2-times", 3, k, x, y); } } else if (C_IS_BIGNUM_TYPE(x)) { if (y & C_FIXNUM_BIT) { C_u_2_integer_times(4, (C_word)NULL, k, x, y); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word *a = C_alloc(C_SIZEOF_FLONUM); C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)*C_flonum_magnitude(y))); } else if (C_IS_BIGNUM_TYPE(y)) { C_u_2_integer_times(4, (C_word)NULL, k, x, y); } else { try_extended_number("numbers#@extended-2-times", 3, k, x, y); } } else { try_extended_number("numbers#@extended-2-times", 3, k, x, y); } } /* Needs SIZEOF_BIGNUM(2) */ C_regparm C_word C_fcall C_a_u_i_2_fixnum_times(C_word **ptr, C_word n, C_word x, C_word y) { C_uword negp, xhi, xlo, yhi, ylo, p, rhi, rlo; negp = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT)); x = (x & C_INT_SIGN_BIT) ? -C_unfix(x) : C_unfix(x); y = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y); xhi = C_BIGNUM_DIGIT_HI_HALF(x); xlo = C_BIGNUM_DIGIT_LO_HALF(x); yhi = C_BIGNUM_DIGIT_HI_HALF(y); ylo = C_BIGNUM_DIGIT_LO_HALF(y); /* This is simply bignum_digits_multiply unrolled for 2x2 halfdigits */ p = xlo * ylo; rlo = C_BIGNUM_DIGIT_LO_HALF(p); p = xhi * ylo + C_BIGNUM_DIGIT_HI_HALF(p); rhi = C_BIGNUM_DIGIT_HI_HALF(p); p = xlo * yhi + C_BIGNUM_DIGIT_LO_HALF(p); rlo = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), rlo); rhi = xhi * yhi + C_BIGNUM_DIGIT_HI_HALF(p) + rhi; if (rhi) { return C_bignum2(ptr, negp != 0, rlo, rhi); } else if (negp ? ((rlo & C_INT_SIGN_BIT) || !C_fitsinfixnump(-(C_word)rlo)) : !C_ufitsinfixnump(rlo)) { return C_bignum1(ptr, negp != 0, rlo); } else { return C_fix(negp ? -rlo : rlo); } } void C_ccall C_u_2_integer_times(C_word c, C_word self, C_word k, C_word x, C_word y) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { C_word *a = C_alloc(C_SIZEOF_BIGNUM(2)); C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y)); } else { C_word tmp = x; /* swap to ensure x is a bignum and y a fixnum */ x = y; y = tmp; } } /* Here, we know for sure that X is a bignum */ if (y == C_fix(0)) { C_kontinue(k, C_fix(0)); } else if (y == C_fix(1)) { C_kontinue(k, x); } else if (y == C_fix(-1)) { C_u_integer_negate(3, (C_word)NULL, k, x); } else if (y & C_FIXNUM_BIT) { /* Any other fixnum */ C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y), negp = C_mk_bool((y & C_INT_SIGN_BIT) ? !C_bignum_negativep(x) : C_bignum_negativep(x)); if (C_fitsinbignumhalfdigitp(absy) || (((C_uword)1 << (C_ilen(absy)-1)) == absy && C_fitsinfixnump(absy))) { C_word size, k2, *a = C_alloc(C_SIZEOF_CLOSURE(4)); k2 = C_closure(&a, 4, (C_word)integer_times_2, k, x, C_fix(absy)); size = C_fix(C_bignum_size(x) + 1); /* Needs _at most_ one more digit */ C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } else { C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM); y = C_a_u_i_fix_to_big(&a, y); bignum_times_bignum_unsigned(k, x, y, negp); } } else { C_word negp = C_bignum_negativep(x) ? !C_bignum_negativep(y) : C_bignum_negativep(y); bignum_times_bignum_unsigned(k, x, y, C_mk_bool(negp)); } } static void integer_times_2(C_word c, C_word self, C_word new_big) { C_word k = C_block_item(self, 1), old_bigx = C_block_item(self, 2), absy = C_unfix(C_block_item(self, 3)); C_uword *digits = C_bignum_digits(new_big), *end_digit = digits + C_bignum_size(old_bigx); int shift; bignum_digits_destructive_copy(new_big, old_bigx); /* Scale up, and sanitise the result. */ shift = C_ilen(absy) - 1; if (((C_uword)1 << shift) == absy) { /* Power of two? */ *end_digit = bignum_digits_destructive_shift_left(digits, end_digit, shift); } else { *end_digit = bignum_digits_destructive_scale_up_with_carry(digits, end_digit, absy, 0); } C_kontinue(k, C_bignum_simplify(new_big)); } static void bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp) { if (C_bignum_size(y) < C_bignum_size(x)) { /* Ensure size(x) <= size(y) */ C_word z = x; x = y; y = z; } if (C_bignum_size(x) < C_KARATSUBA_THRESHOLD) { /* Gradebook */ C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size; k2 = C_closure(&ka, 4, (C_word)bignum_times_bignum_unsigned_2, k, x, y); size = C_fix(C_bignum_size(x) + C_bignum_size(y)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE); } else { try_extended_number("numbers#@bignum-2-times-karatsuba", 3, k, x, y); } } static C_regparm void bignum_digits_multiply(C_word x, C_word y, C_word result) { C_uword product, *xd = C_bignum_digits(x), *yd = C_bignum_digits(y), *rd = C_bignum_digits(result); C_uhword carry, yj; /* Lengths in halfwords */ int i, j, length_x = C_bignum_size(x) * 2, length_y = C_bignum_size(y) * 2; /* From Hacker's Delight, Figure 8-1 (top part) */ for (j = 0; j < length_y; ++j) { yj = C_uhword_ref(yd, j); if (yj == 0) continue; carry = 0; for (i = 0; i < length_x; ++i) { product = (C_uword)C_uhword_ref(xd, i) * yj + (C_uword)C_uhword_ref(rd, i + j) + carry; C_uhword_set(rd, i + j, product); carry = C_BIGNUM_DIGIT_HI_HALF(product); } C_uhword_set(rd, j + length_x, carry); } } static void bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3); bignum_digits_multiply(x, y, result); C_kontinue(k, C_bignum_simplify(result)); } void C_ccall C_digits_to_integer(C_word c, C_word self, C_word k, C_word str, C_word start, C_word end, C_word radix, C_word negp) { assert((C_unfix(radix) > 1) && C_fitsinbignumhalfdigitp(C_unfix(radix))); if (start == end) { C_kontinue(k, C_SCHEME_FALSE); } else { C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2, size; size_t nbits; k2 = C_closure(&ka, 6, (C_word)digits_to_integer_2, k, str, start, end, radix); nbits = (C_unfix(end) - C_unfix(start)) * C_ilen(C_unfix(radix)-1); size = C_fix(C_BIGNUM_BITS_TO_DIGITS(nbits)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } C_inline int hex_char_to_digit(int ch) { if (ch == (int)'#') return 0; /* Hash characters in numbers are mapped to 0 */ else if (ch >= (int)'a') return ch - (int)'a' + 10; /* lower hex */ else if (ch >= (int)'A') return ch - (int)'A' + 10; /* upper hex */ else return ch - (int)'0'; /* decimal (OR INVALID; handled elsewhere) */ } static void digits_to_integer_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), str = C_block_item(self, 2), start = C_unfix(C_block_item(self, 3)), end = C_unfix(C_block_item(self, 4)), radix = C_unfix(C_block_item(self, 5)); char *s = C_c_string(str); C_kontinue(k, str_to_bignum(result, s + start, s + end, radix)); } /* Write from digit character stream to bignum. Bignum does not need * to be initialised. Returns the bignum, or a fixnum. Assumes the * string contains only digits that fit within radix (checked by * string->number). */ static C_regparm C_word str_to_bignum(C_word bignum, char *str, char *str_end, int radix) { int radix_shift, str_digit; C_uword *digits = C_bignum_digits(bignum), *end_digits = digits + C_bignum_size(bignum), big_digit = 0; /* Below, we try to save up as much as possible in big_digit, and * only when it exceeds what we would be able to multiply easily, we * scale up the bignum and add what we saved up. */ radix_shift = C_ilen(radix) - 1; if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */ int n = 0; /* Number of bits read so far into current big digit */ /* Read from least to most significant digit to avoid shifting or scaling */ while (str_end > str) { str_digit = hex_char_to_digit((int)*--str_end); big_digit |= (C_uword)str_digit << n; n += radix_shift; if (n >= C_BIGNUM_DIGIT_LENGTH) { n -= C_BIGNUM_DIGIT_LENGTH; *digits++ = big_digit; big_digit = str_digit >> (radix_shift - n); } } assert(n < C_BIGNUM_DIGIT_LENGTH); /* If radix isn't an exact divisor of digit length, write final digit */ if (n > 0) *digits++ = big_digit; assert(digits == end_digits); } else { /* Not a power of two */ C_uword *last_digit = digits, factor; /* bignum starts as zero */ do { factor = radix; while (str < str_end && C_fitsinbignumhalfdigitp(factor)) { str_digit = hex_char_to_digit((int)*str++); factor *= radix; big_digit = radix * big_digit + str_digit; } big_digit = bignum_digits_destructive_scale_up_with_carry( digits, last_digit, factor / radix, big_digit); if (big_digit) { (*last_digit++) = big_digit; /* Move end */ big_digit = 0; } } while (str < str_end); /* Set remaining digits to zero so bignum_simplify can do its work */ assert(last_digit <= end_digits); while (last_digit < end_digits) *last_digit++ = 0; } return C_bignum_simplify(bignum); } /* TODO: Copied from runtime.c */ # define STRING_BUFFER_SIZE 4096 static C_TLS C_char buffer[ STRING_BUFFER_SIZE ]; static char *to_n_nary(C_uword num, C_uword base, int negp, int as_flonum) { static char *digits = "0123456789abcdef"; char *p; C_uword shift = C_ilen(base) - 1; int mask = (1 << shift) - 1; if (as_flonum) { buffer[68] = '\0'; buffer[67] = '0'; buffer[66] = '.'; } else { buffer[66] = '\0'; } p = buffer + 66; if (mask == base - 1) { do { *(--p) = digits [ num & mask ]; num >>= shift; } while (num); } else { do { *(--p) = digits [ num % base ]; num /= base; } while (num); } if (negp) *(--p) = '-'; return p; } void C_ccall C_basic_number_to_string(C_word c, C_word closure, C_word k, C_word num, ...) { C_word radix; if(c == 3) { radix = C_fix(10); } else if(c == 4) { va_list v; va_start(v, num); radix = va_arg(v, C_word); va_end(v); if(!(radix & C_FIXNUM_BIT)) barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix); } else { C_bad_argc(c, 3); } if(num & C_FIXNUM_BIT) { C_u_fixnum_to_string(4, (C_word)NULL, k, num, radix); } else if (C_immediatep(num)) { barf(C_BAD_ARGUMENT_TYPE_ERROR, "number->string", num); } else if(C_block_header(num) == C_FLONUM_TAG) { C_u_flonum_to_string(4, (C_word)NULL, k, num, radix); } else if (C_IS_BIGNUM_TYPE(num)) { C_u_integer_to_string(4, (C_word)NULL, k, num, radix); } else { try_extended_number("numbers#@extended-number->string", 3, k, num, radix); } } void C_ccall C_u_fixnum_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix) { C_char *p; C_word *a, neg = num & C_INT_SIGN_BIT ? 1 : 0; radix = C_unfix(radix); if (radix < 2 || radix > 16) { barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix); } num = neg ? -C_unfix(num) : C_unfix(num); p = to_n_nary(num, radix, neg, 0); num = C_strlen(p); a = C_alloc((C_bytestowords(num) + 1)); C_kontinue(k, C_string(&a, num, p)); } void C_ccall C_u_flonum_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix) { C_word *a; C_char *p; double f; radix = C_unfix(radix); f = C_flonum_magnitude(num); /* XXX TODO: Should inexacts be printable in other bases than 10? * Perhaps output a string starting with #i? * Right now something like (number->string 1e40 16) results in * a string that can't be read back using string->number. */ if((radix < 2) || (radix > 16)){ barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", C_fix(radix)); } if(C_fits_in_unsigned_int_p(num) == C_SCHEME_TRUE) { /* Use fast int code */ if(f < 0) { p = to_n_nary((C_uword)-f, radix, 1, 1); } else { p = to_n_nary((C_uword)f, radix, 0, 1); } } else if(C_isnan(f)) { p = "+nan.0"; } else if(C_isinf(f)) { p = f > 0 ? "+inf.0" : "-inf.0"; } else { /* Doesn't fit an unsigned int and not "special"; use system libc */ C_snprintf(buffer, STRING_BUFFER_SIZE, C_text("%.*g"), /* XXX: flonum_print_precision */ (int)C_unfix(C_get_print_precision()), f); buffer[STRING_BUFFER_SIZE-1] = '\0'; if((p = C_strpbrk(buffer, C_text(".eE"))) == NULL) { /* Already checked for these, so shouldn't happen */ assert(*buffer != 'i'); /* "inf" */ assert(*buffer != 'n'); /* "nan" */ /* Ensure integral flonums w/o expt are always terminated by .0 */ #if defined(HAVE_STRLCAT) || !defined(C_strcat) C_strlcat(buffer, C_text(".0"), sizeof(buffer)); #else C_strcat(buffer, C_text(".0")); #endif } p = buffer; } radix = C_strlen(p); a = C_alloc((C_bytestowords(radix) + 1)); radix = C_string(&a, radix, p); C_kontinue(k, radix); } /* Naming is a little inconsistent, but looks saner. We're not R-O-B-O-T-S! */ void C_ccall C_u_integer_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix) { if (num & C_FIXNUM_BIT) { C_u_fixnum_to_string(4, (C_word)NULL, k, num, radix); } else { int len, radix_shift; size_t nbits; if ((C_unfix(radix) < 2) || (C_unfix(radix) > 16)) { barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix); } /* Approximation of the number of radix digits we'll need. We try * to be as precise as possible to avoid memmove overhead at the end * of the non-powers of two part of the conversion procedure, which * we may need to do because we write strings back-to-front, and * pointers must be aligned (even for byte blocks). */ len = C_bignum_size(num)-1; nbits = (size_t)len * C_BIGNUM_DIGIT_LENGTH; nbits += C_ilen(C_bignum_digits(num)[len]); len = C_ilen(C_unfix(radix))-1; len = (nbits + len - 1) / len; len += C_bignum_negativep(num) ? 1 : 0; /* Add space for negative sign */ radix_shift = C_ilen(C_unfix(radix)) - 1; if (len > C_RECURSIVE_TO_STRING_THRESHOLD && /* The power of two fast path is much faster than recursion */ ((C_uword)1 << radix_shift) != C_unfix(radix)) { try_extended_number("numbers#@integer->string/recursive", 4, k, num, radix, C_fix(len)); } else { C_word k2, *ka; ka = C_alloc(C_SIZEOF_CLOSURE(4)); k2 = C_closure(&ka, 4, (C_word)bignum_to_str_2, k, num, radix); C_allocate_vector(6, (C_word)NULL, k2, C_fix(len), /* Byte vec, no initialization, no align at 8 bytes */ C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE); } } } static void bignum_to_str_2(C_word c, C_word self, C_word string) { static char *characters = "0123456789abcdef"; C_word k = C_block_item(self, 1), bignum = C_block_item(self, 2), radix = C_unfix(C_block_item(self, 3)); char *buf = C_c_string(string), *index = buf + C_header_size(string) - 1; int radix_shift, negp = (C_bignum_negativep(bignum) ? 1 : 0); radix_shift = C_ilen(radix) - 1; if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */ int radix_mask = radix - 1, big_digit_len = 0, radix_digit; C_uword *scan, *end, big_digit = 0; scan = C_bignum_digits(bignum); end = scan + C_bignum_size(bignum); while (scan < end) { /* If radix isn't an exact divisor of digit length, handle overlap */ if (big_digit_len == 0) { big_digit = *scan++; big_digit_len = C_BIGNUM_DIGIT_LENGTH; } else { assert(index >= buf); radix_digit = big_digit; big_digit = *scan++; radix_digit |= (big_digit << big_digit_len) & radix_mask; big_digit >>= (radix_shift - big_digit_len); big_digit_len = C_BIGNUM_DIGIT_LENGTH - big_digit_len; } while(big_digit_len >= radix_shift && index >= buf) { radix_digit = big_digit & radix_mask; *index-- = characters[radix_digit]; big_digit >>= radix_shift; big_digit_len -= radix_shift; } } assert(big_digit < radix); /* Final digit (like overlap at start of while loop) */ if (big_digit) *index-- = characters[big_digit]; if (negp) { /* Loop above might've overwritten sign position with a zero */ if (*(index+1) == '0') *(index+1) = '-'; else *index-- = '-'; } /* Length calculation is always precise for radix powers of two. */ assert(index == buf-1); } else { C_uword base, *start, *scan, big_digit; C_word working_copy; int steps, i; working_copy = allocate_tmp_bignum(C_fix(C_bignum_size(bignum)), C_mk_bool(negp), C_SCHEME_FALSE); bignum_digits_destructive_copy(working_copy, bignum); start = C_bignum_digits(working_copy); scan = start + C_bignum_size(bignum); /* Calculate the largest power of radix that fits a halfdigit: * steps = log10(2^halfdigit_bits), base = 10^steps */ for(steps = 0, base = radix; C_fitsinbignumhalfdigitp(base); base *= radix) steps++; base /= radix; /* Back down: we overshot in the loop */ while (scan > start) { big_digit = bignum_digits_destructive_scale_down(start, scan, base); if (*(scan-1) == 0) scan--; /* Adjust if we exhausted the highest digit */ for(i = 0; i < steps && index >= buf; ++i) { C_word tmp = big_digit / radix; *index-- = characters[big_digit - (tmp*radix)]; /* big_digit % radix */ big_digit = tmp; } } assert(index >= buf-1); free_tmp_bignum(working_copy); /* Move index onto first nonzero digit. We're writing a bignum here: it can't consist of only zeroes. */ while(*++index == '0'); if (negp) *--index = '-'; /* Shorten with distance between start and index. */ if (buf != index) { i = C_header_size(string) - (index - buf); C_memmove(buf, index, i); /* Move start of number to beginning. */ C_block_header(string) = C_STRING_TYPE | i; /* Mutate strlength. */ } } C_kontinue(k, string); } C_regparm double C_bignum_to_double(C_word bignum) { double accumulator = 0; C_uword *start = C_bignum_digits(bignum), *scan = start + C_bignum_size(bignum); while (start < scan) { accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH; accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH; accumulator += (*--scan); } return(C_bignum_negativep(bignum) ? -accumulator : accumulator); } static void fabs_frexp_to_digits(C_uword exp, double sign, C_uword *start, C_uword *scan) { C_uword digit, odd_bits = exp % C_BIGNUM_DIGIT_LENGTH; assert(C_isfinite(sign)); assert(0.5 <= sign && sign < 1); /* Guaranteed by frexp() and fabs() */ assert((scan - start) == C_BIGNUM_BITS_TO_DIGITS(exp)); if (odd_bits > 0) { /* Handle most significant digit first */ sign *= (C_uword)1 << odd_bits; digit = (C_uword)sign; (*--scan) = digit; sign -= (double)digit; } while (start < scan && sign > 0) { sign *= pow(2.0, C_BIGNUM_DIGIT_LENGTH); digit = (C_uword)sign; (*--scan) = digit; sign -= (double)digit; } /* Finish up by clearing any remaining, lower, digits */ while (start < scan) (*--scan) = 0; } static C_word flo_to_tmp_bignum(C_word x) { /* TODO: allocating and initialising the bignum is pointless if we * already know the number of limbs in the comparand. In fact, * bignum_cmp will first check the number of limbs and *then* * compare. Instead, we can check beforehand and check the limbs * directly against the generated limbs, without allocating at all! */ C_word tmp_big, negp = C_mk_bool(C_flonum_magnitude(x) < 0.0); int exponent; double significand = frexp(C_flonum_magnitude(x), &exponent); assert(C_u_i_fpintegerp_fixed(x)); if (exponent <= 0) { tmp_big = allocate_tmp_bignum(C_fix(0), C_SCHEME_FALSE, C_SCHEME_FALSE); } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */ /* Don't use fix_to_big to simplify caller code: it can just free this */ tmp_big = allocate_tmp_bignum(C_fix(1), negp, C_SCHEME_FALSE); C_bignum_digits(tmp_big)[0] = 1; } else { C_uword size, *start, *end; size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent)); tmp_big = allocate_tmp_bignum(size, negp, C_SCHEME_FALSE); start = C_bignum_digits(tmp_big); end = start + C_bignum_size(tmp_big); fabs_frexp_to_digits(exponent, fabs(significand), start, end); } return tmp_big; } void C_ccall C_u_flo_to_int(C_word c, C_word self, C_word k, C_word x) { int exponent; double significand = frexp(C_flonum_magnitude(x), &exponent); assert(C_truep(C_u_i_fpintegerp_fixed(x))); if (exponent <= 0) { C_kontinue(k, C_fix(0)); } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */ C_kontinue(k, significand < 0.0 ? C_fix(-1) : C_fix(1)); } else { C_word kab[C_SIZEOF_CLOSURE(4) + C_SIZEOF_FLONUM], *ka = kab, k2, size, negp = C_mk_bool(C_flonum_magnitude(x) < 0.0), sign = C_flonum(&ka, fabs(significand)); k2 = C_closure(&ka, 4, (C_word)flo_to_int_2, k, C_fix(exponent), sign); size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } static void flo_to_int_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1); C_uword exponent = C_unfix(C_block_item(self, 2)), *start = C_bignum_digits(result), *scan = start + C_bignum_size(result); double significand = C_flonum_magnitude(C_block_item(self, 3)); fabs_frexp_to_digits(exponent, significand, start, scan); C_kontinue(k, C_bignum_simplify(result)); } C_word C_ccall C_u_i_integer_length(C_word x) { if (x & C_FIXNUM_BIT) { return C_u_i_fixnum_length(x); } else { C_uword result = (C_bignum_size(x) - 1) * C_BIGNUM_DIGIT_LENGTH, *last_digit = C_bignum_digits(x) + C_bignum_size(x) - 1, last_digit_length = C_ilen(*last_digit); /* If *only* the highest bit is set, negating will give one less bit */ if (C_bignum_negativep(x) && *last_digit == ((C_uword)1 << (last_digit_length-1))) { C_uword *startx = C_bignum_digits(x); while (startx < last_digit && *startx == 0) ++startx; if (startx == last_digit) result--; } return C_fix(result + last_digit_length); } } void C_ccall /* x is any exact integer but y is _always_ a fixnum */ C_u_integer_shift(C_word c, C_word self, C_word k, C_word x, C_word y) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; y = C_unfix(y); if (y == 0 || x == C_fix(0)) { /* Done (no shift) */ C_kontinue(k, x); } else if (x & C_FIXNUM_BIT) { if (y < 0) { /* Don't shift more than a word's length (that's undefined in C!) */ if (-y < C_WORD_SIZE) { C_kontinue(k, C_fix(C_unfix(x) >> -y)); } else { C_kontinue(k, (x < 0) ? C_fix(-1) : C_fix(0)); } } else if (y > 0 && y < C_WORD_SIZE-2 && /* After shifting, the length still fits a fixnum */ (C_ilen(C_unfix(x)) + y) < C_WORD_SIZE-2) { C_kontinue(k, C_fix(C_unfix(x) << y)); } else { x = C_a_u_i_fix_to_big(&a, x); } } { C_word ab[C_SIZEOF_CLOSURE(6)], *a = ab, k2, size, negp, digit_offset, bit_offset; negp = C_mk_bool(C_bignum_negativep(x)); if (y > 0) { /* y is guaranteed not to be 0 here */ digit_offset = y / C_BIGNUM_DIGIT_LENGTH; bit_offset = y % C_BIGNUM_DIGIT_LENGTH; k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k, x, C_SCHEME_TRUE, C_fix(digit_offset), C_fix(bit_offset)); size = C_fix(C_bignum_size(x) + digit_offset + 1); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } else if (-y >= C_bignum_size(x) * (C_word)C_BIGNUM_DIGIT_LENGTH) { /* All bits are shifted out, just return 0 or -1 */ C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(0)); } else { digit_offset = -y / C_BIGNUM_DIGIT_LENGTH; bit_offset = -y % C_BIGNUM_DIGIT_LENGTH; k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k, x, C_SCHEME_FALSE, C_fix(digit_offset), C_fix(bit_offset)); size = C_fix(C_bignum_size(x) - digit_offset); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } } C_inline C_word maybe_negate_bignum_for_bitwise_op(C_word x, C_word size) { C_word nx = C_SCHEME_FALSE, xsize; if (C_bignum_negativep(x)) { nx = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE); xsize = C_bignum_size(x); /* Copy up until requested size, and init any remaining upper digits */ C_memcpy(C_bignum_digits(nx), C_bignum_digits(x), C_wordstobytes(nmin(size, xsize))); if (size > xsize) C_memset(C_bignum_digits(nx)+xsize, 0, C_wordstobytes(size-xsize)); bignum_digits_destructive_negate(nx); } return nx; } static void bignum_actual_shift(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), shift_left = C_truep(C_block_item(self, 3)), digit_offset = C_unfix(C_block_item(self, 4)), bit_offset = C_unfix(C_block_item(self, 5)); C_uword *startr = C_bignum_digits(result), *startx = C_bignum_digits(x), *endx = startx + C_bignum_size(x), *endr = startr + C_bignum_size(result); if (shift_left) { /* Initialize only the lower digits we're skipping and the MSD */ C_memset(startr, 0, C_wordstobytes(digit_offset)); *(endr-1) = 0; startr += digit_offset; /* Can't use bignum_digits_destructive_copy because it assumes * we want to copy from the start. */ C_memcpy(startr, startx, C_wordstobytes(endx-startx)); if(bit_offset > 0) bignum_digits_destructive_shift_left(startr, endr, bit_offset); } else { C_word nx, size = C_bignum_size(x) + 1; if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) { startx = C_bignum_digits(nx); /* update startx; x and endx are unused */ } startx += digit_offset; /* Can't use bignum_digits_destructive_copy because that assumes * target is at least as big as source. */ C_memcpy(startr, startx, C_wordstobytes(endr-startr)); if(bit_offset > 0) bignum_digits_destructive_shift_right(startr,endr,bit_offset,C_truep(nx)); if (C_truep(nx)) { free_tmp_bignum(nx); bignum_digits_destructive_negate(result); } } C_kontinue(k, C_bignum_simplify(result)); } /* This is currently only used by Karatsuba multiplication and * Burnikel-Ziegler division. It is not intended as a public API! */ void C_ccall C_u_bignum_extract_digits(C_word c, C_word self, C_word k, C_word x, C_word start, C_word end) { if (x & C_FIXNUM_BIT) { /* Needed? */ if (C_unfix(start) == 0 && (end == C_SCHEME_FALSE || C_unfix(end) > 0)) C_kontinue(k, x); else C_kontinue(k, C_fix(0)); } else { C_word negp, size; negp = C_mk_bool(C_bignum_negativep(x)); /* Always false */ start = C_unfix(start); /* We might get passed larger values than actually fits; pad w/ zeroes */ if (end == C_SCHEME_FALSE) end = C_bignum_size(x); else end = nmin(C_unfix(end), C_bignum_size(x)); assert(start >= 0); size = end - start; if (size == 0 || start >= C_bignum_size(x)) { C_kontinue(k, C_fix(0)); } else { C_word k2, kab[C_SIZEOF_CLOSURE(5)], *ka = kab; k2 = C_closure(&ka, 5, (C_word)bignum_actual_extraction, k, x, C_fix(start), C_fix(end)); C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE); } } } static void bignum_actual_extraction(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), start = C_unfix(C_block_item(self, 3)), end = C_unfix(C_block_item(self, 4)); C_uword *result_digits = C_bignum_digits(result), *x_digits = C_bignum_digits(x); /* Can't use bignum_digits_destructive_copy because that assumes * target is at least as big as source. */ C_memcpy(result_digits, x_digits + start, C_wordstobytes(end-start)); C_kontinue(k, C_bignum_simplify(result)); } C_regparm C_word C_ccall C_u_i_integer_randomize(C_word seed) { /* TODO: Rename C_randomize to C_u_i_fixnum_randomize */ if (seed & C_FIXNUM_BIT) { return C_randomize(seed); } else { /* * This random number generator is very simple. Probably too simple... */ C_uword *scan = C_bignum_digits(seed), *end = scan + C_bignum_size(seed), iseed = 0; /* What a cheap way to initialize the random generator. I feel dirty! */ while (scan < end) iseed ^= *scan++; srand(iseed); return C_SCHEME_UNDEFINED; } } void C_ccall C_u_integer_random(C_word c, C_word self, C_word k, C_word max) { /* TODO: for consistency C_random_fixnum should be called C_u_i_fixnum_random */ if (max & C_FIXNUM_BIT) { C_kontinue(k, C_random_fixnum(max)); } else { C_word k2, kab[C_SIZEOF_CLOSURE(4)], *ka = kab, size, max_len, max_bits, max_top_digit, d, negp; max_len = C_bignum_size(max); max_top_digit = d = C_bignum_digits(max)[max_len - 1]; max_bits = (max_len - 1) * C_BIGNUM_DIGIT_LENGTH; while(d) { max_bits++; d >>= 1; } /* Subtract/add one because we don't want zero to be over-represented */ size = ((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1); size = C_fix(C_BIGNUM_BITS_TO_DIGITS(size + 1)); negp = C_mk_bool(C_bignum_negativep(max)); k2 = C_closure(&ka, 4, (C_word)bignum_random_2, k, C_fix(max_top_digit), C_fix(max_len)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } static void bignum_random_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), max_top_digit = C_unfix(C_block_item(self, 2)), max_len = C_unfix(C_block_item(self, 3)); C_uword *scan = C_bignum_digits(result), *end = scan + C_bignum_size(result); /* Go to just before the end. */ while(scan < end) *scan++ = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH); /* * Last word is special when length is max_len: It must be less than * max's most significant digit, instead of 2^{digitlen}. */ if (max_len == C_bignum_size(result)) *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit; else *scan = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH); C_kontinue(k, C_bignum_simplify(result)); } C_word C_ccall C_u_i_integer_bit_setp(C_word n, C_word i) { if (!(i & C_FIXNUM_BIT)) { /* A bit silly, but might be useful */ return C_u_i_integer_negativep(n); } else if (i & C_INT_SIGN_BIT) { barf(C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR, "bit-set?", n, i); } else { i = C_unfix(i); if (n & C_FIXNUM_BIT) { if (i >= C_WORD_SIZE) return C_mk_bool(n & C_INT_SIGN_BIT); else return C_mk_bool((C_unfix(n) & ((C_word)1 << i)) != 0); } else { C_word nn, d; d = i / C_BIGNUM_DIGIT_LENGTH; if (d >= C_bignum_size(n)) return C_mk_bool(C_bignum_negativep(n)); if (C_truep(nn = maybe_negate_bignum_for_bitwise_op(n, d))) n = nn; i %= C_BIGNUM_DIGIT_LENGTH; d = C_mk_bool((C_bignum_digits(n)[d] & (C_uword)1 << i) != 0); if (C_truep(nn)) free_tmp_bignum(nn); return d; } } } void C_ccall C_u_2_integer_bitwise_and(C_word c, C_word self, C_word k, C_word x, C_word y) { if ((x & y) & C_FIXNUM_BIT) { C_kontinue(k, C_u_fixnum_and(x, y)); } else { C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2; if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x); if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y); negp = C_mk_bool(C_bignum_negativep(x) && C_bignum_negativep(y)); /* Allow negative 1-bits to propagate */ if (C_bignum_negativep(x) || C_bignum_negativep(y)) size = nmax(C_bignum_size(x), C_bignum_size(y)) + 1; else size = nmin(C_bignum_size(x), C_bignum_size(y)); ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */ k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_and_2, k, x, y); C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE); } } static void bignum_bitwise_and_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3), size = C_bignum_size(result), nx, ny; C_uword *scanr = C_bignum_digits(result), *endr = scanr + C_bignum_size(result), *scans1, *ends1, *scans2; if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx; if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny; if (C_bignum_size(x) < C_bignum_size(y)) { scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x); scans2 = C_bignum_digits(y); } else { scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y); scans2 = C_bignum_digits(x); } while (scans1 < ends1) *scanr++ = *scans1++ & *scans2++; C_memset(scanr, 0, C_wordstobytes(endr - scanr)); if (C_truep(nx)) free_tmp_bignum(nx); if (C_truep(ny)) free_tmp_bignum(ny); if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result); C_kontinue(k, C_bignum_simplify(result)); } void C_ccall C_u_2_integer_bitwise_ior(C_word c, C_word self, C_word k, C_word x, C_word y) { if ((x & y) & C_FIXNUM_BIT) { C_kontinue(k, C_u_fixnum_or(x, y)); } else { C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2; if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x); if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y); ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */ k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_ior_2, k, x, y); size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1); negp = C_mk_bool(C_bignum_negativep(x) || C_bignum_negativep(y)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } static void bignum_bitwise_ior_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3), size = C_bignum_size(result), nx, ny; C_uword *scanr = C_bignum_digits(result), *endr = scanr + C_bignum_size(result), *scans1, *ends1, *scans2, *ends2; if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx; if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny; if (C_bignum_size(x) < C_bignum_size(y)) { scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x); scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y); } else { scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y); scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x); } while (scans1 < ends1) *scanr++ = *scans1++ | *scans2++; while (scans2 < ends2) *scanr++ = *scans2++; if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */ assert(scanr == endr); if (C_truep(nx)) free_tmp_bignum(nx); if (C_truep(ny)) free_tmp_bignum(ny); if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result); C_kontinue(k, C_bignum_simplify(result)); } void C_ccall C_u_2_integer_bitwise_xor(C_word c, C_word self, C_word k, C_word x, C_word y) { if ((x & y) & C_FIXNUM_BIT) { C_kontinue(k, C_fixnum_xor(x, y)); } else { C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, size, k2, negp; if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x); if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y); ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */ k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_xor_2, k, x, y); size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1); negp = C_mk_bool(C_bignum_negativep(x) != C_bignum_negativep(y)); C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE); } } static void bignum_bitwise_xor_2(C_word c, C_word self, C_word result) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3), size = C_bignum_size(result), nx, ny; C_uword *scanr = C_bignum_digits(result), *endr = scanr + C_bignum_size(result), *scans1, *ends1, *scans2, *ends2; if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx; if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny; if (C_bignum_size(x) < C_bignum_size(y)) { scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x); scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y); } else { scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y); scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x); } while (scans1 < ends1) *scanr++ = *scans1++ ^ *scans2++; while (scans2 < ends2) *scanr++ = *scans2++; if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */ assert(scanr == endr); if (C_truep(nx)) free_tmp_bignum(nx); if (C_truep(ny)) free_tmp_bignum(ny); if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result); C_kontinue(k, C_bignum_simplify(result)); } static void bignum_digits_destructive_negate(C_word result) { C_uword *scan, *end, digit, sum; scan = C_bignum_digits(result); end = scan + C_bignum_size(result); do { digit = ~*scan; sum = digit + 1; *scan++ = sum; } while (sum == 0 && scan < end); for (; scan < end; scan++) { *scan = ~*scan; } } /* This is ugly but really cleans up the code below */ #define RETURN_Q_AND_OR_R(calc_q, calc_r) \ if (C_truep(C_and(return_q, return_r))) { \ C_values(4, C_SCHEME_UNDEFINED, k, calc_q, calc_r); \ } else if (C_truep(return_r)) { \ C_kontinue(k, calc_r); \ } else { \ C_kontinue(k, calc_q); \ } /* Lossy; we could be in "quotient&remainder" or "modulo" */ #define DIVREM_LOC ((C_truep(C_and(return_q, return_r))) ? "/" : \ (C_truep(return_q) ? "quotient" : "remainder")) static C_regparm void basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) { if (x & C_FIXNUM_BIT) { if (y & C_FIXNUM_BIT) { C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC); RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y), C_u_i_fixnum_remainder_checked(x, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word ab[C_SIZEOF_FLONUM*3], *a = ab; if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC); x = C_a_i_fix_to_flo(&a, 1, x); RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y), C_a_i_flonum_remainder_checked(&a, 2, x, y)); } else if (C_IS_BIGNUM_TYPE(y)) { integer_divrem(6, (C_word)NULL, k, x, y, return_q, return_r); } else { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y); } } else if (C_immediatep(x)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, x); } else if (C_block_header(x) == C_FLONUM_TAG) { if (!C_truep(C_u_i_fpintegerp_fixed(x))) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x); } else if (y & C_FIXNUM_BIT) { C_word ab[C_SIZEOF_FLONUM*3], *a = ab; if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC); y = C_a_i_fix_to_flo(&a, 1, y); RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y), C_a_i_flonum_remainder_checked(&a, 2, x, y)); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y); } else if (C_block_header(y) == C_FLONUM_TAG) { C_word ab[C_SIZEOF_FLONUM*2], *a = ab; if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC); RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y), C_a_i_flonum_remainder_checked(&a, 2, x, y)); } else if (C_IS_BIGNUM_TYPE(y)) { C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab; x = flo_to_tmp_bignum(x); k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, x); integer_divrem(6, (C_word)NULL, k2, x, y, return_q, return_r); } else { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y); } } else if (C_IS_BIGNUM_TYPE(x)) { if (y & C_FIXNUM_BIT) { integer_divrem(6, (C_word)NULL, k, x, y, return_q, return_r); } else if (C_immediatep(y)) { barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y); } else if (C_block_header(y) == C_FLONUM_TAG) { if (!C_truep(C_u_i_fpintegerp_fixed(y))) { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y); } else if (C_flonum_magnitude(y) == 0.0) { numbers_div_by_zero_error(DIVREM_LOC); } else { C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab; y = flo_to_tmp_bignum(y); k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, y); integer_divrem(6, (C_word)NULL, k2, x, y, return_q, return_r); } } else if (C_IS_BIGNUM_TYPE(y)) { bignum_divrem(6, (C_word)NULL, k, x, y, return_q, return_r); } else { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y); } } else { barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x); } } static void divrem_intflo_2(C_word c, C_word self, ...) { C_word k = C_block_item(self, 1), x, y; va_list v; free_tmp_bignum(C_block_item(self, 2)); va_start(v, self); if (c == 2) { C_word ab[C_SIZEOF_FLONUM], *a = ab; x = va_arg(v, C_word); va_end(v); if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x); else x = C_a_u_i_big_to_flo(&a, 1, x); C_kontinue(k, x); } else { /* c == 3 */ C_word ab[C_SIZEOF_FLONUM*2], *a = ab; x = va_arg(v, C_word); y = va_arg(v, C_word); va_end(v); if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x); else x = C_a_u_i_big_to_flo(&a, 1, x); if (y & C_FIXNUM_BIT) y = C_a_i_fix_to_flo(&a, 1, y); else y = C_a_u_i_big_to_flo(&a, 1, y); C_values(4, C_SCHEME_UNDEFINED, k, x, y); } } static void bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big) { C_word k = C_block_item(self, 1), return_q = C_block_item(self, 2), return_r = C_block_item(self, 3); RETURN_Q_AND_OR_R(negated_big, C_fix(0)); } static C_regparm void integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) { if (!(y & C_FIXNUM_BIT)) { /* y is bignum. */ if (x & C_FIXNUM_BIT) { /* abs(x) < abs(y), so it will always be [0, x] except for this case: */ if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && C_bignum_negated_fitsinfixnump(y)) { RETURN_Q_AND_OR_R(C_fix(-1), C_fix(0)); } else { RETURN_Q_AND_OR_R(C_fix(0), x); } } else { bignum_divrem(6, (C_word)NULL, k, x, y, return_q, return_r); } } else if (x & C_FIXNUM_BIT) { /* both x and y are fixnum. */ C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC); RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y), C_u_i_fixnum_remainder_checked(x, y)); } else { /* x is bignum, y is fixnum. */ C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y); if (y == C_fix(1)) { RETURN_Q_AND_OR_R(x, C_fix(0)); } else if (y == C_fix(0)) { numbers_div_by_zero_error(DIVREM_LOC); } else if (y == C_fix(-1)) { C_word *ka, k2; ka = C_alloc(C_SIZEOF_CLOSURE(4)); k2 = C_closure(&ka, 4, (C_word)bignum_divrem_fixnum_2, k, return_q, return_r); C_u_integer_negate(3, (C_word)NULL, k2, x); } else if (C_fitsinbignumhalfdigitp(absy) || ((((C_uword)1 << (C_ilen(absy)-1)) == absy) && C_fitsinfixnump(absy))) { if (C_truep(return_q)) { C_word q_negp = C_mk_bool((y & C_INT_SIGN_BIT) ? !(C_bignum_negativep(x)) : C_bignum_negativep(x)), r_negp = C_mk_bool(C_bignum_negativep(x)), *ka, k2, size; ka = C_alloc(C_SIZEOF_CLOSURE(7)); size = C_fix(C_bignum_size(x)); k2 = C_closure(&ka, 7, (C_word)bignum_destructive_divide_unsigned_small, k, x, C_fix(absy), return_q, return_r, r_negp); C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE); } else { C_word rem; C_uword next_power = (C_uword)1 << (C_ilen(absy)-1); if (next_power == absy) { /* Is absy a power of two? */ rem = *(C_bignum_digits(x)) & (next_power - 1); } else { /* Too bad, we have to do some real work */ rem = bignum_remainder_unsigned_halfdigit(x, absy); } C_kontinue(k, C_bignum_negativep(x) ? C_fix(-rem) : C_fix(rem)); } } else { /* Just divide it as two bignums */ C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab; bignum_divrem(6, (C_word)NULL, k, x, C_a_u_i_fix_to_big(&a, y), return_q, return_r); } } } static C_regparm void bignum_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) { C_word q_negp = C_mk_bool(C_bignum_negativep(y) ? !C_bignum_negativep(x) : C_bignum_negativep(x)), r_negp = C_mk_bool(C_bignum_negativep(x)), size; switch(bignum_cmp_unsigned(x, y)) { case 0: RETURN_Q_AND_OR_R(C_truep(q_negp) ? C_fix(-1) : C_fix(1), C_fix(0)); case -1: RETURN_Q_AND_OR_R(C_fix(0), x); case 1: default: size = C_bignum_size(x) - C_bignum_size(y); if (size > C_BURNIKEL_ZIEGLER_DIFF_THRESHOLD) { try_extended_number("numbers#@bignum-2-divrem-burnikel-ziegler", 5, k, x, y, return_q, return_r); } else if (C_truep(return_q)) { C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2; k2 = C_closure(&ka, 6, (C_word)bignum_divide_2_unsigned, k, x, y, return_r, r_negp); size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y)); C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE); } else { /* We can skip bignum_divide_2_unsigned if we need no quotient */ C_word kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2; k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k, x, y, return_q, return_r, C_SCHEME_UNDEFINED); size = C_fix(C_bignum_size(x) + 1); /* May need to be normalized */ C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE); } } } static C_word bignum_remainder_unsigned_halfdigit(C_word num, C_word den) { C_uword *start = C_bignum_digits(num), *scan = start + C_bignum_size(num), rem = 0, two_digits; assert((den > 1) && (C_fitsinbignumhalfdigitp(den))); while (start < scan) { two_digits = (*--scan); rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_HI_HALF(two_digits)) % den; rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_LO_HALF(two_digits)) % den; } return rem; } /* External interface for the above internal divrem functions */ void C_ccall C_basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y) { if (c != 4) C_bad_argc_2(c, 4, self); basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE); } void C_ccall C_u_integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y) { integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE); } void C_ccall C_basic_remainder(C_word c, C_word self, C_word k, C_word x, C_word y) { if (c != 4) C_bad_argc_2(c, 4, self); basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE); } void C_ccall C_u_integer_remainder(C_word c, C_word self, C_word k, C_word x, C_word y) { integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE); } void C_ccall C_basic_quotient(C_word c, C_word self, C_word k, C_word x, C_word y) { if (c != 4) C_bad_argc_2(c, 4, self); basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE); } void C_ccall C_u_integer_quotient(C_word c, C_word self, C_word k, C_word x, C_word y) { integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE); } /* "small" is either a number that fits a halfdigit, or a power of two */ static void bignum_destructive_divide_unsigned_small(C_word c, C_word self, C_word quotient) { C_word k = C_block_item(self, 1), numerator = C_block_item(self, 2), denominator = C_unfix(C_block_item(self, 3)), /* return_quotient = C_block_item(self, 4), */ return_remainder = C_block_item(self, 5), remainder_negp = C_block_item(self, 6); C_uword *start = C_bignum_digits(quotient), *end = start + C_bignum_size(quotient), remainder; int shift_amount; bignum_digits_destructive_copy(quotient, numerator); shift_amount = C_ilen(denominator)-1; if (((C_uword)1 << shift_amount) == denominator) { /* Power of two? Shift! */ remainder = bignum_digits_destructive_shift_right(start,end,shift_amount,0); assert(C_ufitsinfixnump(remainder)); } else { remainder = bignum_digits_destructive_scale_down(start, end, denominator); assert(C_fitsinbignumhalfdigitp(remainder)); } quotient = C_bignum_simplify(quotient); if (C_truep(return_remainder)) { remainder = C_truep(remainder_negp) ? -remainder : remainder; C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder)); } else { C_kontinue(k, quotient); } } /* Full bignum division */ static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient) { C_word k = C_block_item(self, 1), x = C_block_item(self, 2), y = C_block_item(self, 3), size = C_fix(C_bignum_size(x) + 1), return_r = C_block_item(self, 4), r_negp = C_block_item(self, 5), kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2; k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k, x, y, C_SCHEME_TRUE, return_r, quotient); C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE); } /* 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". [Yeah, that's a nice book but that particular section is not helpful at all, which is also pointed out by P. Brinch Hansen's "Multiple-Length Division Revisited: A Tour Of The Minefield". That's a more down-to-earth step-by-step explanation of the algorithm. Add to this the C implementation in Hacker's Delight (section 9-2, p141--142) and you may be able to grok this... ...barely, if you're as math-challenged as I am -- sjamaan] */ static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder) { C_word k = C_block_item(self, 1), numerator = C_block_item(self, 2), denominator = C_block_item(self, 3), return_quotient = C_block_item(self, 4), return_remainder = C_block_item(self, 5), quotient = C_block_item(self, 6), length = C_bignum_size(denominator); C_uword d1 = *(C_bignum_digits(denominator) + length - 1), *startr = C_bignum_digits(remainder), *endr = startr + C_bignum_size(remainder); int shift; shift = C_BIGNUM_DIGIT_LENGTH - C_ilen(d1); /* nlz */ /* We have to work on halfdigits, so we shift out only the necessary * amount in order fill out that halfdigit (base is halved). * This trick is shamelessly stolen from Gauche :) * See below for part 2 of the trick. */ if (shift >= C_BIGNUM_HALF_DIGIT_LENGTH) shift -= C_BIGNUM_HALF_DIGIT_LENGTH; /* Code below won't always set high halfdigit of quotient, so do it here. */ if (quotient != C_SCHEME_UNDEFINED) C_bignum_digits(quotient)[C_bignum_size(quotient)-1] = 0; bignum_digits_destructive_copy(remainder, numerator); *(endr-1) = 0; /* Ensure most significant digit is initialised */ if (shift == 0) { /* Already normalized */ bignum_destructive_divide_normalized(remainder, denominator, quotient); } else { /* Requires normalisation; allocate scratch denominator for this */ C_uword *startnd; C_word ndenom; bignum_digits_destructive_shift_left(startr, endr, shift); ndenom = allocate_tmp_bignum(C_fix(length), C_SCHEME_FALSE, C_SCHEME_FALSE); startnd = C_bignum_digits(ndenom); bignum_digits_destructive_copy(ndenom, denominator); bignum_digits_destructive_shift_left(startnd, startnd+length, shift); bignum_destructive_divide_normalized(remainder, ndenom, quotient); if (C_truep(return_remainder)) /* Otherwise, don't bother shifting back */ bignum_digits_destructive_shift_right(startr, endr, shift, 0); free_tmp_bignum(ndenom); } if (C_truep(return_remainder)) { if (C_truep(return_quotient)) { C_values(4, C_SCHEME_UNDEFINED, k, C_bignum_simplify(quotient), C_bignum_simplify(remainder)); } else { C_kontinue(k, C_bignum_simplify(remainder)); } } else { assert(C_truep(return_quotient)); C_kontinue(k, C_bignum_simplify(quotient)); } } static void bignum_destructive_divide_normalized(C_word big_u, C_word big_v, C_word big_q) { C_uword *v = C_bignum_digits(big_v), *u = C_bignum_digits(big_u), *q = big_q == C_SCHEME_UNDEFINED ? NULL : C_bignum_digits(big_q), p, /* product of estimated quotient & "denominator" */ hat, qhat, rhat, /* estimated quotient and remainder digit */ vn_1, vn_2; /* "cached" values v[n-1], v[n-2] */ C_word t, k; /* Two helpers: temp/final remainder and "borrow" */ /* We use plain ints here, which theoretically may not be enough on * 64-bit for an insanely huge number, but it is a _lot_ faster. */ int n = C_bignum_size(big_v) * 2, /* in halfwords */ m = (C_bignum_size(big_u) * 2) - 2; /* Correct for extra digit */ int i, j; /* loop vars */ /* Part 2 of Gauche's aforementioned trick: */ if (C_uhword_ref(v, n-1) == 0) n--; /* These won't change during the loop, but are used in every step. */ vn_1 = C_uhword_ref(v, n-1); vn_2 = C_uhword_ref(v, n-2); /* See also Hacker's Delight, Figure 9-1. This is almost exactly that. */ for (j = m - n; j >= 0; j--) { hat = C_BIGNUM_DIGIT_COMBINE(C_uhword_ref(u, j+n), C_uhword_ref(u, j+n-1)); if (hat == 0) { if (q != NULL) C_uhword_set(q, j, 0); continue; } qhat = hat / vn_1; rhat = hat % vn_1; /* Two whiles is faster than one big check with an OR. Thanks, Gauche! */ while(qhat >= (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) { qhat--; rhat += vn_1; } while(qhat * vn_2 > C_BIGNUM_DIGIT_COMBINE(rhat, C_uhword_ref(u, j+n-2)) && rhat < (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) { qhat--; rhat += vn_1; } /* Multiply and subtract */ k = 0; for (i = 0; i < n; i++) { p = qhat * C_uhword_ref(v, i); t = C_uhword_ref(u, i+j) - k - C_BIGNUM_DIGIT_LO_HALF(p); C_uhword_set(u, i+j, t); k = C_BIGNUM_DIGIT_HI_HALF(p) - (t >> C_BIGNUM_HALF_DIGIT_LENGTH); } t = C_uhword_ref(u,j+n) - k; C_uhword_set(u, j+n, t); if (t < 0) { /* Subtracted too much? */ qhat--; k = 0; for (i = 0; i < n; i++) { t = (C_uword)C_uhword_ref(u, i+j) + C_uhword_ref(v, i) + k; C_uhword_set(u, i+j, t); k = t >> C_BIGNUM_HALF_DIGIT_LENGTH; } C_uhword_set(u, j+n, (C_uhword_ref(u, j+n) + k)); } if (q != NULL) C_uhword_set(q, j, qhat); } /* end j */ }