diff --git a/bignum.c b/bignum.c index fe22aaeb..942e7200 100644 --- a/bignum.c +++ b/bignum.c @@ -480,7 +480,8 @@ sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) { sexp_bignum_sign(b1) = 1; q = SEXP_ZERO; x = sexp_make_bignum(ctx, sexp_bignum_length(a)); - while (sexp_bignum_compare(a1, b1) > 0) { + while (sexp_bignum_compare_abs(a1, b1) > 0) { + /* guess divisor x */ alen = sexp_bignum_hi(a1); d = sexp_bignum_data(a1)[alen-1] / sexp_bignum_data(b1)[blen-1]; memset(sexp_bignum_data(x), 0, sexp_bignum_length(x)*sizeof(sexp_uint_t)); @@ -491,14 +492,26 @@ sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) { } else { sexp_bignum_data(x)[alen - blen] = d; } + /* update quotient q and remainder a1 estimates */ y = sexp_bignum_mul(ctx, NULL, b1, x); - a1 = sexp_bignum_sub(ctx, NULL, a1, y); - q = (sign < 0) ? sexp_sub(ctx, q, x) : sexp_add(ctx, q, x); - if (sexp_bignum_sign(a1) < 0) { - sexp_bignum_sign(a1) = 1; + if (sign < 0) { + a1 = sexp_bignum_add(ctx, NULL, a1, y); + q = sexp_sub(ctx, q, x); + } else { + a1 = sexp_bignum_sub(ctx, NULL, a1, y); + q = sexp_add(ctx, q, x); + } + /* flip the sign if we overshot in our estimate */ + if (sexp_bignum_sign(a1) != sign) { + sexp_bignum_sign(a1) = -sign; sign *= -1; } } + /* adjust signs */ + if (sign < 0) { + q = sexp_sub(ctx, q, SEXP_ONE); + a1 = sexp_add(ctx, a1, b1); + } if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) { sexp_negate_exact(q); } diff --git a/tests/numeric-tests.scm b/tests/numeric-tests.scm index 80362b9c..a89b0c5f 100644 --- a/tests/numeric-tests.scm +++ b/tests/numeric-tests.scm @@ -176,6 +176,10 @@ (test 1.0 (max 1/2 1.0)) (test 18446744073709551617 (numerator (/ 18446744073709551617 2))) (test "18446744073709551617/2" (number->string (/ 18446744073709551617 2))) + (let ((a 1000000000000000000000000000000000000000) + (b 31622776601683794000)) + (test 31622776601683792639 (quotient a b)) + (test 30922992657207634000 (remainder a b))) (let ((r (/ (expt 2 61) 3))) (test 0 (- r r)) (test 2305843009213693952/3 r))) diff --git a/tests/r7rs-tests.scm b/tests/r7rs-tests.scm index 2fa8311a..2a3bf916 100644 --- a/tests/r7rs-tests.scm +++ b/tests/r7rs-tests.scm @@ -224,6 +224,10 @@ (let*-values (((root rem) (exact-integer-sqrt (expt 2 121)))) (list root rem))) +(test '(31622776601683793319 62545769258890964239) + (let*-values (((root rem) (exact-integer-sqrt (expt 10 39)))) + (list root rem))) + (let*-values (((root rem) (exact-integer-sqrt (expt 2 140)))) (test 0 rem) (test (expt 2 140) (square root)))