improve precision of / on complex numbers with small parts

Closes #2654
This commit is contained in:
Matthew Flatt 2019-05-14 07:56:02 -06:00
parent 1c8d547dbb
commit ddb944d583
2 changed files with 22 additions and 0 deletions

View File

@ -837,6 +837,7 @@
(test 0.2+0.4i / 1.0-2.0i)
(test 0.4-0.2i / 2.0+1.0i)
(test 0.4+0.2i / 2.0-1.0i)
(test 0.0+0.0i / 0.0+0.0i 1+1e-320i)
(test 3 / 1 1/3)
(test -3 / 1 -1/3)
@ -1916,6 +1917,8 @@
(test +nan.0+nan.0i expt -1.5 +inf.0)
(test +nan.0+nan.0i expt -1.5 -inf.0)
(test 1.0-9.9998886718268e-321i expt 1+1.0e-320i -1)
(err/rt-test (expt 0 -1) exn:fail:contract:divide-by-zero?)
(err/rt-test (expt 0 -1.0) exn:fail:contract:divide-by-zero?)
(err/rt-test (expt 0 -inf.0) exn:fail:contract:divide-by-zero?)

View File

@ -258,6 +258,25 @@ Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Objec
r = scheme_bin_div(c, d);
/* If r goes to infinity, try computing a different way to avoid overflow: */
if (SCHEME_FLOATP(r)) {
double v = SCHEME_FLOAT_VAL(r);
if (MZ_IS_POS_INFINITY(v) || MZ_IS_NEG_INFINITY(v)) {
/* From Chez Scheme: a+bi / c+di => (ac+bd)/(cc+dd) + ((bc-ad)/(cc+dd))i
This calculuation does not work as well for complex numbers with
large parts, such as `(/ 1e+300+1e+300i 4e+300+4e+300i)`, but it
works better for small parts, as in `(/ 0.0+0.0i 1+1e-320i)`. */
cm = scheme_bin_plus(scheme_bin_mult(c, c),
scheme_bin_mult(d, d));
return scheme_make_complex(scheme_bin_div(scheme_bin_plus(scheme_bin_mult(a, c),
scheme_bin_mult(b, d)),
cm),
scheme_bin_div(scheme_bin_minus(scheme_bin_mult(b, c),
scheme_bin_mult(a, d)),
cm));
}
}
den = scheme_bin_plus(d, scheme_bin_mult(c, r));
if (swap)