summary: complex division sometimes overflows reported: John W. Eaton , 2002-03-05 affects: branch_release-1-6 see-also-dir: 1003.stuff patch: denominator-special-casing.diff test-case: complex-divs.scm John W. Eaton sez: > The textbook formula for complex division > > (rx * ry + ix * iy) + (ix * ry - rx * iy) i > (/ x y) == ------------------------------------------- > ry * ry + iy * iy > > that guile currently uses can cause trouble because overflow can occur > unnecessarily when computing the denominator. A better formula for > this operation doesn't form the sum of the squares of the real and > imaginary parts of the denominator. > > Interestingly, this bug seems to be masked by gcc's optimizations. It > only appeared for me when I compiled without -O2. For example, with > -O2, guile produced > > guile> (/ 1e200+1e200i) > 5.0e-201-5.0e-201i > > which is (mysteriously) correct, but without -O2, it computes > > guile> (/ 1e200+1e200i) > 0.0 > > (the denominator is infinite, 1/inf ==> 0). > > I've appended a patch, relative to the current CVS archive. It > affects three cases in scm_divide: > > (/ ) > (/ ) > (/ ) > > Here are tests for some cases that are easier to think about to make > sure I haven't introduced any new problems. > > (define d1 3+4i) ;; r <= i > (define d2 4+3i) ;; r > i > (define cn 25+125i) > > (/ d1) ;; 0.12-0.16i single arg, case 1 > (/ d2) ;; 0.16-0.12i single arg, case 2 > (/ cn d1) ;; 23.0+11.0i / , case 1 > (/ cn d2) ;; 19.0+17.0i / , case 2 > (/ 25 d1) ;; 3.0-4.0i / , case 1 > (/ 25 d2) ;; 4.0-3.0i / , case 2 > > (There are now separate branches in the code to handle the cases > > real part of the denominator <= imaginary part of the denominator > real part of the denominator > imaginary part of the denominator > > so the tests cover both possibilities.) > > Since this code is ultimately derived from f2c, I've included the > copyright notice, which requires that it appear in the code and > documentation. I'm not sure where it should appear, so I haven't > provided a patch for the docs. ttn sez: other branches need to be checked. we need to rfc from numerical methods gurus. how does this change fit into overall GMP strategy?