From bug-request at octave dot org Mon Feb 7 10:02:38 2005 Subject: Re: Inconsistency in complex division by zero From: David Bateman To: "John W. Eaton" CC: bug at octave dot org Date: Mon, 07 Feb 2005 17:04:32 +0100 John W. Eaton wrote: > >And I would say both are wrong, because i is pure imaginary, so we >should really be getting just Infi, no? When we have 1/0, we get Inf, >not Inf+NaNi (i.e., there is no imaginary part of 0 that is magically >tacked on to the 1 to screw things up). > >Another example that demonstrates this problem is Inf*i, which you >would expect to produce Infi, not Nan+Infi. > >To solve this correctly, we need a pure imaginary data type or perhaps >special cases in many of the operators in the complex class. > > > As the /0 case only occurs in "div" function in op-{cs,cm}-s.cc and in the "ldiv" function in op-s-{cs,cm}.cc perhaps just a special case there is sufficient. We already have a test there to give the warning message. The Inf*i case only happens in the same places. The CM*M case uses ZGEMM, promoting the real matrix to a complex matrix and so doesn't have this issue... So the fix is pretty easy and isolated. For example DEFBINOP (div, complex_matrix, scalar) { CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_scalar&); double d = v2.double_value (); if (d == 0.0) { gripe_divide_by_zero (); return octave_value (v1.complex_array_value () / Complex (d)); } else return octave_value (v1.complex_array_value () / d); } DEFBINOP (mul, complex_matrix, scalar) { CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_scalar&); return octave_value (v1.complex_array_value () * Complex (d)); } Regards David -- David Bateman David dot Bateman at motorola dot com Motorola Labs - Paris +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax) 91193 Gif-Sur-Yvette FRANCE The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------