Fcalva/chibi-scheme
1
0
Fork 0
mirror of https://github.com/ashinn/chibi-scheme.git synced 2025-05-19 05:39:18 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

More bignum division fixes.

Browse source
This commit is contained in:
Alex Shinn 2014-01-25 23:52:39 +09:00
parent 48fe3c8014
commit c17a30942f
2 changed files with 39 additions and 11 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

45
bignum.c
View file

@ -468,30 +468,53 @@ sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b) {
}
sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
sexp_uint_t d;
sexp_sint_t alen, blen=sexp_bignum_hi(b), sign=1;
sexp_luint_t d;
sexp_uint_t dlo, dhi;
sexp_sint_t alen, blen=sexp_bignum_hi(b), sign=1, off=0;
sexp_gc_var5(q, x, y, a1, b1);
if (blen == 1 && sexp_bignum_data(b)[0] == 0)
return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
sexp_gc_preserve5(ctx, q, x, y, a1, b1);
a1 = sexp_copy_bignum(ctx, NULL, a, 0);
sexp_bignum_sign(a1) = 1;
/* fast path for single bigit divisor */
if (blen == 1) {
b1 = sexp_make_bignum(ctx, 1);
sexp_bignum_data(b1)[0] = sexp_bignum_fxdiv(ctx, a1, sexp_bignum_data(b)[0], 0);
*rem = sexp_bignum_normalize(b1);
if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
sexp_negate_exact(a1);
}
if (sexp_bignum_sign(a) < 0) {
sexp_negate_exact(*rem);
}
sexp_gc_release5(ctx);
return a1;
}
/* general case */
b1 = sexp_copy_bignum(ctx, NULL, b, 0);
sexp_bignum_sign(b1) = 1;
q = SEXP_ZERO;
x = sexp_make_bignum(ctx, sexp_bignum_length(a));
while (sexp_bignum_compare_abs(a1, b1) > 0) {
while (sexp_bignum_compare_abs(a1, b1) > 0) { /* a1, b1 at least 2 bigits */
/* 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));
sexp_bignum_data(x)[off] = 0;
if (off > 0) sexp_bignum_data(x)[off-1] = 0;
off = alen - blen + 1;
d = (((sexp_luint_t)sexp_bignum_data(a1)[alen-1] << (sizeof(sexp_uint_t)*8))
+ sexp_bignum_data(a1)[alen-2])
/ (((sexp_luint_t)sexp_bignum_data(b1)[blen-1] << (sizeof(sexp_uint_t)*8))
+ sexp_bignum_data(b1)[blen-2]);
if (d == 0) {
sexp_bignum_data(x)[alen - blen - 1] =
(((sexp_luint_t)sexp_bignum_data(a1)[alen-1] << (sizeof(sexp_uint_t)*8))
+ sexp_bignum_data(a1)[alen-2]) / sexp_bignum_data(b1)[blen-1];
} else {
sexp_bignum_data(x)[alen - blen] = d;
d = ((sexp_luint_t)sexp_bignum_data(a1)[alen-1] << (sizeof(sexp_uint_t)*8)) + (sexp_luint_t)sexp_bignum_data(a1)[alen-2];
d /= sexp_bignum_data(b1)[blen-1];
off--;
}
dhi = d >> (sizeof(sexp_uint_t)*8);
dlo = d & (((sexp_luint_t)1<<(sizeof(sexp_uint_t)*8))-1);
sexp_bignum_data(x)[off] = dhi;
if (off > 0) sexp_bignum_data(x)[off-1] = dlo;
/* update quotient q and remainder a1 estimates */
y = sexp_bignum_mul(ctx, NULL, b1, x);
if (sign < 0) {
@ -512,10 +535,10 @@ sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
q = sexp_sub(ctx, q, SEXP_ONE);
a1 = sexp_add(ctx, a1, b1);
}
*rem = a1;
if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
sexp_negate_exact(q);
}
*rem = a1;
if (sexp_bignum_sign(a) < 0) {
sexp_negate_exact(*rem);
}

5
tests/numeric-tests.scm
View file

@ -167,6 +167,8 @@
(test #f (>= -inf.0 +inf.0))
(test #f (> -inf.0 +inf.0))
;;(test (expt 10 154) (sqrt (expt 10 308)))
(cond-expand
(ratios
(test #t (< 1/2 1.0))
@ -180,6 +182,9 @@
(b 31622776601683794000))
(test 31622776601683792639 (quotient a b))
(test 30922992657207634000 (remainder a b)))
(let ((g 18446744073709551616/6148914691236517205))
(test 36893488147419103231/113427455640312821148309287786019553280
(- g (/ 9 g))))
(let ((r (/ (expt 2 61) 3)))
(test 0 (- r r))
(test 2305843009213693952/3 r)))