avoid potential overflow in complex division

[John W. Eaton]

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:

    (/ <complex>)
    (/ <complex> <complex>)
    (/ <real> <complex>)

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  <complex> / <complex>, case 1
  (/ cn d2)  ;; 19.0+17.0i  <complex> / <complex>, case 2
  (/ 25 d1)  ;; 3.0-4.0i    <real> / <complex>, case 1
  (/ 25 d2)  ;; 4.0-3.0i    <real> / <complex>, 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.

jwe

-- 
www.octave.org        | Unfortunately we were hoplessly optimistic in 1954
www.che.wisc.edu/~jwe | about the problems of debugging FORTRAN programs.
                      |                                       -- J. Backus



2002-03-05  John W. Eaton  <address@hidden>

        * numbers.c (scm_divide): Adapt code from libstdc++/f2c to void
        potential overflow problems.


Index: libguile/numbers.c
===================================================================
RCS file: /cvs/guile/guile-core/libguile/numbers.c,v
retrieving revision 1.157
diff -u -r1.157 numbers.c
--- libguile/numbers.c  1 Mar 2002 00:19:20 -0000       1.157
+++ libguile/numbers.c  6 Mar 2002 03:46:53 -0000
@@ -3691,6 +3691,33 @@
 #undef FUNC_NAME
 
 
+/* The code below for complex division is adapted from the GNU
+   libstdc++, which adapted it from f2c's libF77, and is subject to
+   this copyright:  */
+
+/****************************************************************
+Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
+
+Permission to use, copy, modify, and distribute this software
+and its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the names of AT&T Bell Laboratories or
+Bellcore or any of their entities not be used in advertising or
+publicity pertaining to distribution of the software without
+specific, written prior permission.
+
+AT&T and Bellcore disclaim all warranties with regard to this
+software, including all implied warranties of merchantability
+and fitness.  In no event shall AT&T or Bellcore be liable for
+any special, indirect or consequential damages or any damages
+whatsoever resulting from loss of use, data or profits, whether
+in an action of contract, negligence or other tortious action,
+arising out of or in connection with the use or performance of
+this software.
+****************************************************************/
+
 SCM_GPROC1 (s_divide, "/", scm_tc7_asubr, scm_divide, g_divide);
 /* Divide the first argument by the product of the remaining
    arguments.  If called with one argument @var{z1}, 1/@var{z1} is
@@ -3724,8 +3751,15 @@
     } else if (SCM_COMPLEXP (x)) {
       double r = SCM_COMPLEX_REAL (x);
       double i = SCM_COMPLEX_IMAG (x);
-      double d = r * r + i * i;
-      return scm_make_complex (r / d, -i / d);
+      if (r <= i) {
+       double t = r / i;
+       double d = i * (1.0 + t * t);
+       return scm_make_complex (t / d, -1.0 / d);
+      } else {
+       double t = i / r;
+       double d = r * (1.0 + t * t);
+       return scm_make_complex (1.0 / d, -t / d);
+      }
     } else {
       SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
     }
@@ -3765,8 +3799,15 @@
       {
        double r = SCM_COMPLEX_REAL (y);
        double i = SCM_COMPLEX_IMAG (y);
-       double d = r * r + i * i;
-       return scm_make_complex ((a * r) / d, (-a * i) / d);
+       if (r <= i) {
+         double t = r / i;
+         double d = i * (1.0 + t * t);
+         return scm_make_complex ((a * t) / d,  -a / d);
+       } else {
+         double t = i / r;
+         double d = r * (1.0 + t * t);
+         return scm_make_complex (a / d,  -(a * t) / d);
+       }
       }
     } else {
       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
@@ -3870,9 +3911,15 @@
     } else if (SCM_COMPLEXP (y)) {
       double ry = SCM_COMPLEX_REAL (y);
       double iy = SCM_COMPLEX_IMAG (y);
-      double d = ry * ry + iy * iy;
-      return scm_make_complex ((rx * ry + ix * iy) / d, 
-                              (ix * ry - rx * iy) / d);
+      if (ry <= iy) {
+       double t = ry / iy;
+       double d = iy * (1.0 + t * t);
+       return scm_make_complex ((rx * t + ix) / d, (ix * t - rx) / d);
+      } else {
+       double t = iy / ry;
+       double d = ry * (1.0 + t * t);
+       return scm_make_complex ((rx + ix * t) / d, (ix - rx * t) / d);
+      }
     } else {
       SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }