Re: Machine epsilon, Octave and Matlab

[John W. Eaton]

On 30-Aug-2001, Christoph Spiel <address@hidden> wrote:

| It is not the compiler that evaluates the expression, but Octave,
| i.e., the interpreter.

Right, and Octave stores all results of operations (even temporaries)
into octave_value objects, so there is little chance that the width of
floating point registers is a factor here, though I suppose it is
possible if your compiler is really tricky.

The original expression was:

  1/(1-x) - 1/(1+x)

For the case of x == eps/2:

The term 1/(1+x) must evaluate to 1, because eps is the smallest
floating point number such that 1+x is not equal to 1, and eps/2
is smaller than that (i.e., 1+x must be 1, and 1/1 should be too).

The term 1/(1-x) could be diffrerent from 1 if 1-x is different from
1.  On most systems, the smallest floating point number such that 1-x
is not equal to 1 is not the same as eps.  On x86 systems, I believe
it is actually the same as eps/2.  So, for this case, we have 1/y
where y is a floating point number that is just smaller than 1:

  GNU Octave, version 2.0.16.92 (i386-pc-linux-gnu).
  Copyright (C) 1996, 1997, 1998, 1999, 2000 John W. Eaton.
  This is free software with ABSOLUTELY NO WARRANTY.
  For details, type `warranty'.

  octave:1> format hex; x = eps/2; 1, 1-x
  ans = 3ff0000000000000
  ans = 3fefffffffffffff

However, on this platform, dividing 1 by this number produces 1
instead of something just larger than 1.  I'd bet that on the SGI
system, either 1-x produces a different result from the one above,
or, if it is the same, the division produces a number that is not
equal to 1.

Can the original poster please run the following code on their SGI
system?

  format hex
  on = 1
  e2 = eps/2
  yn = 1-e2
  yp = 1+e2
  t1 = 1/yn
  t2 = 1/yp
  df = t1-t2

Here is the result I see on the Intel and AMD systems I have access
to:

  on = 3ff0000000000000
  e2 = 3ca0000000000000
  yn = 3fefffffffffffff
  yp = 3ff0000000000000
  t1 = 3ff0000000000000
  t2 = 3ff0000000000000
  df = 0000000000000000

On the Alpha systems that I have access to (running Debian), I see
the following results:

  on = 3ff0000000000000
  e2 = 3ca0000000000000
  yn = 3fefffffffffffff
  yp = 3ff0000000000000
  t1 = 3ff0000000000001
  t2 = 3ff0000000000000
  df = 3cb0000000000000

I expect the results on the mips-sgi system to be the same as for
the Alpha system.

Thanks,

jwe



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