9

I'm trying to evaluate the machine epsilon of my computer (see below). I wrote this:

eps = 1.0;
p = 0;
While[(1.0 + eps) > 1.0, eps = eps/2.0; p += 1];
eps
p
eps*2

and got:

1.42109*10^-14
46
2.84217*10^-14

which is not what I expected (IEEE754). So I did this

EngineeringForm[1.0 + 2 * eps, 20]
EngineeringForm[1.0 + eps, 20]
EngineeringForm[1.0 + (eps/2.0), 20]
EngineeringForm[1.0 + (eps/4.0), 20]
EngineeringForm[1.0 + (eps/16.0), 20]
EngineeringForm[1.0 + (eps/32.0), 20]

and got that:

1.000000000000028
1.000000000000014
1.000000000000007
1.000000000000004
1.000000000000001
1.

which seems to contradict the termination condition of the 'while' loop. So I would like to understand what's going on under Mathematica's hood. I tried to look for some information on Mathematica's internal float representation without success. Any help to explain these results is welcome. I'm sorry if my question is a bit vague.

Eric

Context:

SystemInformationData[{"Kernel" ->
    {"Version" -> "9.0 for Mac OS X x86 (64-bit) (January 24, 2013)", 
    "ReleaseID" -> "9.0.1.0 (4055646, 4055073)", ...}]

EDIT

It was suggested the following program that works directly for some reason:

eps = 1.0;
n = 0;
While[(1.0 + eps) - 1.0 > 0.0,
  eps = eps/2;
  n++;
];
eps*2
n
Log[2, eps*2]
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Eric
  • 93
  • 3
  • How about $MachineEpsilon? – MikeLimaOscar May 30 '14 at 16:37
  • Hmm, yes. That gives 2.22045*10^-16 which is yet another value. Besides, I was asked to 'calculate' the machine's epsilon not just provide the value. Thx for the feedback. – Eric May 30 '14 at 16:47
  • Try 1. + $MachineEpsilon > 1.... – Michael E2 May 30 '14 at 17:13
  • As per reference doc, this is indeed corresponds to 2^(-n+1) with n=53, which is the mantissa of a double precision IEEE754 float. That's what I expected. I still don't understand why the while loop calculation did not give this result. – Eric May 30 '14 at 17:14
  • @MichaelE2, tried that and got: In[439]:= 1.0 + $MachineEpsilon > 1.0 Out[439]= False. If I'm not wrong, this should have been True, no? – Eric May 30 '14 at 17:16
  • To comment on my last comment above, this test works: (1.0 + $MachineEpsilon) - 1.0 > 0.0 – Eric May 31 '14 at 02:58

1 Answers1

10

There is a tolerance Internal`$EqualTolerance that is applied when testing Real numbers. If the numbers agree up to the last Internal`$EqualTolerance digits, then they are treated as equal. Try this:

eps = 1.0;
p = 0;
Block[{Internal`$EqualTolerance = 0.},
 While[(1.0 + eps) > 1.0, eps = eps/2.0; p += 1];
 ]
eps
p
eps*2
(*
  1.11022*10^-16
  53
  2.22045*10^-16
*)

See How to make the computer consider two numbers equal up to a certain precision and also this answer by Alexey Popkov.

Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747