Greetings, and thanks so much for your report! Something similar going
in now!
Take care,
Daniel Gildea <gil...@cs.rochester.edu> 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 c...@maguirefamily.org
==========================================================================
"The earth is but one country, and mankind its citizens." -- Baha'u'llah
