suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix complex divide: avoiding overflow/underflow ended up with the sign wrong in one case

Browse Source 
svn: r4606
This commit is contained in:
Matthew Flatt 2006-10-14 22:40:24 +00:00
parent 585d319eac
commit 332195b88c
2 changed files with 29 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/mzscheme/src/complex.c
View File

@ -214,6 +214,7 @@ Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Objec
Scheme_Complex *cn = (Scheme_Complex *)_n;
Scheme_Complex *cd = (Scheme_Complex *)_d;
Scheme_Object *den, *r, *i, *a, *b, *c, *d, *cm, *dm, *aa[1];
int swap;
if ((cn->r == zero) && (cn->i == zero))
return zero;
@ -234,6 +235,21 @@ Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Objec
return scheme_make_complex(r, i);
}
if (!SCHEME_FLOATP(c) && !SCHEME_FLOATP(d)) {
/* The simple way: */
cm = scheme_bin_plus(scheme_bin_mult(c, c),
scheme_bin_mult(d, d));
r = scheme_bin_div(scheme_bin_plus(scheme_bin_mult(c, a),
scheme_bin_mult(d, b)),
cm);
i = scheme_bin_div(scheme_bin_minus(scheme_bin_mult(c, b),
scheme_bin_mult(d, a)),
cm);
return scheme_make_complex(r, i);
}
aa[0] = d;
if (SCHEME_TRUEP(scheme_zero_p(1, aa))) {
/* This is like dividing by a real number, except that
@ -270,12 +286,17 @@ Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Objec
cm = c;
c = d;
d = cm;
}
swap = 1;
} else
swap = 0;
r = scheme_bin_div(c, d);
den = scheme_bin_plus(d, scheme_bin_mult(c, r));
i = scheme_bin_div(scheme_bin_minus(a, scheme_bin_mult(b, r)),
i = scheme_bin_div((swap
? scheme_bin_minus(a, scheme_bin_mult(b, r))
: scheme_bin_minus(scheme_bin_mult(b, r), a)),
den);
r = scheme_bin_div(scheme_bin_plus(b, scheme_bin_mult(a, r)),
den);

5
src/mzscheme/src/number.c
View File

@ -1908,6 +1908,11 @@ scheme_expt(int argc, Scheme_Object *argv[])
if (SCHEME_NUMBERP(e))
return n;
}
if (SCHEME_RATIONALP(e)
&& (((Scheme_Rational *)e)->num == scheme_exact_one)
&& (((Scheme_Rational *)e)->denom == scheme_make_integer(2))) {
return scheme_sqrt(1, argv);
}
if (n == zeroi) {
/* Power of exact zero */