Re: precision of calculations

[Ted Harding]

On 23-Jul-99 A. Scottedward Hodel wrote:
> double precision arithmetic has about 16 significant digits
> (depending on hardware).  Hence, if (A-B) should be zero, but
> you get
>     norm(A-B)/norm(A) approximately = 1e-16
> you're o.k.  The relative error can creep up as problem dimension
> and problem condition (how sensitive the results are to small
> perturbations in data) increases.

Good reply. In those cases where these tiny discrepancies really matter,
I'd suggest doing the following kind of test at critical points (example
code only -- there'll be better ways of doing it in real life):

  ....

  x = A-B;
  if( abs(x) < 5*eps ), x=0; endif
  ....

In practice, I reckon that this will make the right correction more often
than the wrong one: in real life, (A-B) is much more likely to be
precisely zero, than to be that small but non-zero.

By the way, on my octave eps =  2.2204e-16.

I suggest a factor like 5 (but definitively) because smallish
factors arise very easily indeed; e.g.

  2 - sqrt(2)*sqrt(2) = -4.4409e-16   [ = -2*eps ]

This will heavily complicate what might normally be straightforward code,
and may slow it down a lot; but when it matters, it matters!

Ted.

PS Interestingly,

  2 - (2.^(1/3))*(2.^(1/3))*(2.^(1/3)) = 0

while

  2 - (2.^(1/4))*(2.^(1/4))*(2.^(1/4))*(2.^(1/4)) = -2*eps again!

You never know, do you?





--------------------------------------------------------------------
E-Mail: (Ted Harding) <address@hidden>
Date: 23-Jul-99                                       Time: 16:54:28
------------------------------ XFMail ------------------------------



---------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.  To ensure
that development continues, see www.che.wisc.edu/octave/giftform.html
Instructions for unsubscribing: www.che.wisc.edu/octave/archive.html
---------------------------------------------------------------------