From: "Michael" <>
Newsgroups: comp.soft-sys.matlab
Subject: Re: floating point calc differences between Matlab and C?
Date: Tue, 31 Mar 2009 18:21:01 +0000 (UTC)
Organization: Circular Logic
Lines: 47
Message-ID: <gqtmud$e35$>
References: <gqeaq9$56$> <gqec1m$c0$> <gqepkn$s8r$> <gqfcth$o5u$> <gqrp7g$nu0$> <gqs317$8gs$>
Reply-To: "Michael" <>
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: 1238523661 14437 (31 Mar 2009 18:21:01 GMT)
NNTP-Posting-Date: Tue, 31 Mar 2009 18:21:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1151890
Xref: comp.soft-sys.matlab:529099

"Roger Stafford" <> wrote in message <gqs317$8gs$>...
> "Michael" <> wrote in message <gqrp7g$nu0$>...
>   Michael, I have time to answer only a few of your questions here.
> Q1:
> > Is this the issue, e.g. with a = -1 and b = 0.01: 
> > sqrt(a^2+b^2) evals to 1.0001, so you've lost some of b's original precision because b^2 is added to a^2 which is   orders of magnitude greater. Then you add -1 so have 0.0001 which greatly reduces the order of mag but you can't regain the original precision of b, so you have a number with a precision less than usual for its order of mag. From then on, your relative precision is much less than what you'd expect for the order of mag of the answer?
>   Yes, that is the gist of it.  The error in x is more or less the same as for other values of a and b, but x itself is very much smaller, so the relative error is larger.
> Q2:
> > I measure like so:
> > 	err = abs( (Matlab_result - C_result) / abs(Matlab_result) )
> > Then I check if err is less than the relative eps:
> > 	err < eps( abs(Matlab_result) ) ? OK : problem
>   No, that doesn't look right.  You should either compare
>  err = abs((Matlab_result-C_result)/abs(Matlab_result))
> with just 'eps', or you should compare
>  err = abs(Matlab_result-C_result)
> with eps(Matlab_result).  This last quantity is simply the value of a single least bit error in 'Matlab_result'.  If you have a sequence of several operations, you could get more than one such bit worth of errors.
> Q3:
> > So my absolute value and divisions algorithms are still giving results different than matlab's by greater than eps.
>   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.
>   I cannot see why there should be trouble with absolute value of complex numbers unless the square root routine for reals is defective.
>   I'm afraid that's all I have time for now.  My alibi is that the April 15 deadline for income tax reports is looming very large.
> Roger Stafford

Thanks again Roger. And please, you don't need any alibi, you've been very helpful with your time. I'll digest what you've written here and best of luck with your tax return. I had a weird fit of motivation this month and have miracurously prepared my return already. Strange days.