19

I am trying to determine the double machine epsilon in Java, using the definition of it being the smallest representable double value x such that 1.0 + x != 1.0, just as in C/C++. According to wikipedia, this machine epsilon is equal to 2^-52 (with 52 being the number of double mantissa bits - 1).

My implementation uses the Math.ulp() function:

double eps = Math.ulp(1.0);
System.out.println("eps = " + eps);
System.out.println("eps == 2^-52? " + (eps == Math.pow(2, -52)));

and the results are what I expected:

eps = 2.220446049250313E-16
eps == 2^-52? true

So far, so good. However, if I check that the given eps is indeed the smallest x such that 1.0 + x != 1.0, there seems to be a smaller one, aka the previous double value according to Math.nextAfter():

double epsPred = Math.nextAfter(eps, Double.NEGATIVE_INFINITY);
System.out.println("epsPred = " + epsPred);
System.out.println("epsPred < eps? " + (epsPred < eps));
System.out.println("1.0 + epsPred == 1.0? " + (1.0 + epsPred == 1.0));

Which yields:

epsPred = 2.2204460492503128E-16
epsPred < eps? true
1.0 + epsPred == 1.0? false

As we see, we have a smaller than machine epsilon such which, added to 1, yields not 1, in contradiction to the definition.

So what is wrong with the commonly accepted value for machine epsilon according to this definition? Or did I miss something? I suspect another esoteric aspect of floating-point maths, but I can't see where I went wrong...

EDIT: Thanks to the commenters, I finally got it. I actually used the wrong definition! eps = Math.ulp(1.0) computes the distance to the smallest representable double > 1.0, but -- and that's the point -- that eps is not the smallest x with 1.0 + x != 1.0, but rather about twice that value: Adding 1.0 + Math.nextAfter(eps/2) is rounded up to 1.0 + eps.

Franz D.
  • 929
  • 7
  • 21

2 Answers2

21

using the definition of it being the smallest representable double value x such that 1.0 + x != 1.0, just as in C/C++

This has never been the definition, not in Java and not in C and not in C++.

The definition is that the machine epsilon is the distance between one and the smallest float/double larger than one.

Your “definition” is wrong by a factor of nearly 2.

Also, the absence of strictfp only allows a larger exponent range and should not have any impact on the empirical measurement of epsilon, since that is computed from 1.0 and its successor, each of which and the difference of which can be represented with the standard exponent range.

Pascal Cuoq
  • 75,447
  • 6
  • 144
  • 260
  • Is `1.0 + x != 1.0 && (1.0 + x) - 1.0 == x` an adequate definition? Or "1.0 minus (the smallest double value greater than 1.0)"? – Random832 Feb 26 '15 at 14:08
  • 2
    @Random832 If you want to “compute” it, use `Math.nextAfter(1.0, Double.POSITIVE_INFINITY) - 1.0`. The first proposal is a property that is true for many numbers including 1.0 (this is true: 1.0 + 1.0 != 1.0 && (1.0 + 1.0) - 1.0 == 1.0) – Pascal Cuoq Feb 26 '15 at 14:15
  • So the proposed definition (1.0 + x != 1.0) is invalid because of rounding? Is that the case? – Brian J Feb 26 '15 at 15:34
  • @BrianJ Yes. `1.0 + x` is generally inexact for values of `x` of the magnitude we are concerned with here. By contrast, though the definition by subtraction leaves this fact implicit, the distance between `0x1.0000000000001` and `0x1.0000000000000` can be represented as a `double` and `0x1.0000000000001 - 0x1.0000000000000`, as a subtraction of close floating-point numbers, is exact and computes exactly this value. – Pascal Cuoq Feb 26 '15 at 15:38
  • @Pascal: I don't think there is "the" definition of machine precision, I just used one of the definitions given in the cited wikipedia article, which according to wikipedia (I'm not a numerics specialist), is the one used for C/C++. However, I don't see the difference between your definition and mine, since yours also boils down to `Math.ulp(1.0)`. Again, I may be wrong -- I just want to know *why*. – Franz D. Feb 26 '15 at 15:52
  • 1
    @FranzD. The wikipedia page claims that the C standard uses that definition. It does not. I **actually cite** the C standard in my blog post. Please **do check the reference for yourself**: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf (the sentence to look for is given in the blog post). The wrong definition gives the wrong result, as anyone can check for themselves (and you have). The correct definition is used in glibc. What else do you need? – Pascal Cuoq Feb 26 '15 at 16:06
  • What else? Maybe the ability to read more attentively :) Got it now, the Java definition is simply not correct. Thanks. – Franz D. Feb 26 '15 at 16:33
7

I'm not sure your experimental method / theory is sound. The documentation for the Math class states:

For a given floating-point format, an ulp of a specific real number value is the distance between the two floating-point values bracketing that numerical value

The documentation for the ulp method says:

An ulp of a double value is the positive distance between this floating-point value and the double value next larger in magnitude

So, if you want the smallest eps value such that 1.0 + eps != 1.0, your eps really should generally be less than Math.ulp(1.0), since at least for any value greater than Math.ulp(1.0) / 2, the result will be rounded up.

I think the smallest such value will be given by Math.nextAfter(eps/2, 1.0).

davmac
  • 18,558
  • 1
  • 32
  • 55
  • The mathematical ulp(x) is defined for any real x, so for most values of x, it doesn't matter whether the definition has ≤ or <. a="" agreed="" case="" for="" generally="" http:="" in="" is="" it="" matter="" means="" not="" of="" on="" or="" power="" quite="" there="" thorough="" treatment="" two="" ulp="" what="" when="" whether="" which="" would="" x=""> – Pascal Cuoq Feb 26 '15 at 13:38
  • @davmac: I'm also not sure that my approach is sound, that's why I asked :) But I don't understand: Why isn't `ulp(1.0)`, "the positive distance between 1.0 and the double next larger in magnitude", the smallest value such that `1.0 + eps != 1.0`? There shouldn't be any smaller value with that property, as then `ulp(1.0)` *isn't* the distance to the *next* larger double. But, as @Pascal said, there doesn't seem to be a general definition of `ulp(1.0)`, so maybe my problem lies here. Still I don't get it... – Franz D. Feb 26 '15 at 16:02
  • @FranzD. The value given by the “definition” as “smallest representable double value x such that 1.0 + x != 1.0” is not a power of two, so it is not the result of the ulp function for any definition of the mathematical ulp, nor for the Java method named `ulp`. – Pascal Cuoq Feb 26 '15 at 16:10
  • @Pascal: Huh? Who said ulp(x) returns only powers of two? I only said that the input x = 1.0 is a power of two, so ulp's results may be dubious. – Franz D. Feb 26 '15 at 16:23
  • 1
    @FranzD. I said that ulp returns powers of two. You do not have to believe me. Did you read http://www.ens-lyon.fr/LIP/Pub/Rapports/RR/RR2005/RR2005-09.pdf ? – Pascal Cuoq Feb 26 '15 at 16:30
  • I believe you, but I hadn't read the link... It always amazes me how complex floating-point computation is, if you look closer. – Franz D. Feb 26 '15 at 16:44
  • @FranzD. "I don't understand: Why isn't ulp(1.0), "the positive distance between 1.0 and the double next larger in magnitude", the smallest value such that 1.0 + eps != 1.0?" - because of rounding. (This is hinted at in my answer). `ulp` returns the difference between a value and the next largest (or smallest) value. But, you can take the original value, and add a smaller-than-ulp value, and the result might be _closer_ to the next-largest value than it was to the original value, and so it will be rounded up. – davmac Feb 26 '15 at 17:27
  • Yeah, got my problem now and summarized it in an edit to the question. Thank you all for your patience and insight. – Franz D. Feb 26 '15 at 17:39