Path: news.mathworks.com!not-for-mail
From: "Tim Davis" <davis@cise.ufl.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: floating point calc differences between Matlab and C?
Date: Tue, 31 Mar 2009 21:45:04 +0000 (UTC)
Organization: University of Florida
Lines: 17
Message-ID: <gqu2t0$r3d$1@fred.mathworks.com>
References: <gqeaq9$56$1@fred.mathworks.com> <gqec1m$c0$1@fred.mathworks.com> <gqepkn$s8r$1@fred.mathworks.com> <gqfcth$o5u$1@fred.mathworks.com> <gqrp7g$nu0$1@fred.mathworks.com> <gqs317$8gs$1@fred.mathworks.com>
Reply-To: "Tim Davis" <davis@cise.ufl.edu>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1238535904 27757 172.30.248.35 (31 Mar 2009 21:45:04 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 31 Mar 2009 21:45:04 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45902
Xref: news.mathworks.com comp.soft-sys.matlab:529181

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message ...

>   In the division of complex numbers, you should not expect the error in the real part of the result to be within a bit or so within that real part.  The same can be true with the imaginary part though not at the same time.  All you can hope for is that each of these errors amounts to only a bit or so within the absolute value of the complex result.  If you look at the expressions for the real and imaginary parts, you can see why:
> 
>  real((a+b*i)/(c+d*i)) = (a*c+b*d)/sqrt(c^2+d^2)
>  imag((a+b*i)/(c+d*i)) = (b*c-a*d)/sqrt(c^2+d^2)
> 
> The quantity a*c+b*d or b*c-a*d can be very small compared with the values in a, b, c, and d if a+b*i and c+d*i are nearly orthogonal or nearly parallel, respectively, in the complex plane and thereby lose relative accuracy if this is taken with respect to the real or imaginary results, respectively.  As I say, the errors should be compared with the result's absolute value instead.


Roger brings up the basic problem.  Those expressions are not the best way to do complex division.  See for example Algorithm 116, ACM Transactions on Mathematical Software, by R. L. Smith (1962).  It tries to avoid underflow and overflow.  You can see my implementation of it here:

http://www.cise.ufl.edu/research/sparse/umfpack/UMFPACK/Source/umfpack_global.c

It's very different than the standard algebraic description.  My guess is that the gcc compiler might have less accurate code for complex division than it could have.

Getting these basic operations (complex multiplication and division) done correctly is a fine art.