Re: Machine epsilon, Octave and Matlab

[Thomas Shores]

Craig Stoudt wrote:

> Interesting.  I'm running Octave 2.1.31 under Win2000.
>  Running your script I get:
>
> x2 =  1.1102230246251565404e-16
> y2 = 0
> x4 =  5.5511151231257827021e-17
> y4 = 0
>
> I wonder if this was changed in version 2.1.33
>
> Craig Stoudt
>
> __________________________________________________
> Do You Yahoo!?
> Get email alerts & NEW webcam video instant messaging with Yahoo! Messenger
> http://im.yahoo.com
>
> -------------------------------------------------------------
> 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
> -------------------------------------------------------------

Running Octave 2.1.34 under RedHat Linux 7.0 I get the same results as Craig.
I think the difference might be really a compiler issue.  Pentium architecture
is capable of giving the (nearly) correct result for the following reasons:
for small x,
      1/(1 - x)  =  1 + x + x^2 + x^3 +  ...
      1/(1 + x)  =  1 - x + x^2 -x^3 +  ...
so that
     1/(1 - x) - 1/(1 + x) = 2*(x + x^3 + ... )
Now eps = 1/2^52 in IEEE double precision.  Furthermore, the Pentium stores
numbers in internal registers that give a precision of 64 bits and larger
exponent range, which is more than enough to do the division 1/(1 +/- x) ,
derive the 1 -/+ x  part of the correct answer and store it in an internal
register.  Thus, if these results are stored in internal registers and added,
there is ample range to correctly  subtract and store an approximate answer of
2*x.  Note that when x=eps, 1+x^3 = 1 in machine arithmetic, even with larger
registers for arithmetic, so this higher order part of the answer has no
chance.

So what's the problem?  Well, *if* the compiler is smart enough to take a
compound expression like the one above and store intermediate results in the
Pentium FPU registers, the answer should come out OK.  I'm guessing that for
whatever reason this doesn't happen on the Pentium, but does happen on the SGI
compiler.  Perhaps its a question of optimization on one compiler and not on
the other.

Tom Shores







-------------------------------------------------------------
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
-------------------------------------------------------------