division of complex numbers

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);
 	}