division of complex numbers

[Daniel Gildea]

Here is a patch with a more precision-aware implementation of
division of complex numbers.  I based this code on the implementation
in sbcl.  This patch fixes an error in the maxima test suite:

  Running tests in rtest_elliptic: 
  ***************** rtest_elliptic.mac: Problem 232 (line 1071) 
*****************
  
  Input:
  closeto(jacobi_am(1.5, 1.5 + %i) - (0.9340542168700783
                                            - 0.3723960452146072 %i), 
5.4673E-16)
  
  
  Result:
  1.0990647210786426E-15
  
  This differed from the expected result:
  true
  

diff --git a/gcl/o/num_arith.c b/gcl/o/num_arith.c
index a31743cba..fdc7b4b38 100644
--- a/gcl/o/num_arith.c
+++ b/gcl/o/num_arith.c
@@ -1006,18 +1006,32 @@ number_divide(object x, object y)
 
                x = number_to_complex(x);
                y = number_to_complex(y);
-               z1 = number_times(y->cmp.cmp_real, y->cmp.cmp_real);
-               z2 = number_times(y->cmp.cmp_imag, y->cmp.cmp_imag);
-               z3 = number_plus(z1, z2);
-               /* if (number_zerop(z3 = number_plus(z1, z2))) 
DIVISION_BY_ZERO(sLD,list(2,x,y)); */
-               z1 = number_times(x->cmp.cmp_real, y->cmp.cmp_real);
-               z2 = number_times(x->cmp.cmp_imag, y->cmp.cmp_imag);
-               z1 = number_plus(z1, z2);
-               z = number_times(x->cmp.cmp_imag, y->cmp.cmp_real);
-               z2 = number_times(x->cmp.cmp_real, y->cmp.cmp_imag);
-               z2 = number_minus(z, z2);
-               z1 = number_divide(z1, z3);
-               z2 = number_divide(z2, z3);
+               if ( 1 == number_compare(number_abs(y->cmp.cmp_real),
+                                        number_abs(y->cmp.cmp_imag)) ) {
+                 object r = number_divide(y->cmp.cmp_imag, y->cmp.cmp_real);
+                 object dn = number_plus(y->cmp.cmp_real,
+                                         number_times(r, y->cmp.cmp_imag));
+                 z1 = number_divide(number_plus(x->cmp.cmp_real,
+                                                number_times(x->cmp.cmp_imag,
+                                                             r)),
+                                    dn);
+                 z2 = number_divide(number_minus(x->cmp.cmp_imag,
+                                                 number_times(x->cmp.cmp_real,
+                                                              r)),
+                                    dn);
+               } else {
+                 object r = number_divide(y->cmp.cmp_real, y->cmp.cmp_imag);
+                 object dn = number_plus(y->cmp.cmp_imag,
+                                         number_times(r, y->cmp.cmp_real));
+                 z1 = number_divide(number_plus(number_times(x->cmp.cmp_real,
+                                                             r),
+                                                x->cmp.cmp_imag),
+                                    dn);
+                 z2 = number_divide(number_minus(number_times(x->cmp.cmp_imag,
+                                                              r),
+                                                x->cmp.cmp_real),
+                                    dn);
+               }
                z = make_complex(z1, z2);
                return(z);
        }