Re: division of complex numbers

Greetings, and thanks so much for your report!  Something similar going
in now!

Take care,

Daniel Gildea <[email protected]> writes:

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

-- 
Camm Maguire			     		    [email protected]
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah