5

I am attempting to calculate the machine epsilon value for doubles and floats in C++ as part of a school assignment. I'm using Cygwin in Windows 7, 64 bit, here is the code:

#include <iostream>

int main() {
    double epsilon = 1;
    while(1 + epsilon > 1)
        epsilon = epsilon / 2;
    epsilon = 2*epsilon;

    std::cout << epsilon << std::endl;

    float epsilon_f = 1;
    while(1 + epsilon_f > 1)
        epsilon_f = epsilon_f / 2;
    epsilon_f = 2*epsilon_f;

    std::cout << epsilon_f << std::endl;

    return 1;
}

When I run the code, I receive 1.0842e-019 for both values. I looked it up and should be getting 2.22e-16 for the double, and 1.19e-07 for the float value. When I run the exact same code on a Macbook, the code returns the correct values. What could be causing the discrepancy on my Windows machine?

JennaLewis
  • 51
  • 2
  • 1
    Try changing -mfpmath to sse or 387, and make sure both are in 64 bit mode with -m64. cygwin may default to something different than the macbook. Also -fno-fast-math – John Meacham Feb 06 '15 at 02:14
  • If you're still around, and if you found my answer helpful, please mark it as accepted - that way, other people won't stumble across this question later and think that it still needs answering. Also, feel free to ask for more information if the answer didn't help. – Aasmund Eldhuset Feb 13 '15 at 23:35

1 Answers1

1

The CPU's floating-point registers typically contain 80 bits, and it looks like the Cygwin compiler chooses to perform the loop calculations entirely in registers (only truncating the result to 32/64 bits when printing the results).

As @Potatoswatter points out, this is entirely legal for the compiler, and your program actually exhibits undefined behavior because it assumes that there is a precision limit. Since there's undefined behavior, the compiler may choose to transform your program into anything it wants (including one that deletes all of your files, but that's fortunately not a common solution...)

P.S. Welcome to StackOverflow, and kudos for a question that (if you read up on the concepts in the answers) will probably make you learn more about processor architecture and compilers than anyone else in your class! :-)

Aasmund Eldhuset
  • 37,289
  • 4
  • 68
  • 81
  • 1
    … and more to the point, the compiler is allowed to do this. C and C++ do not limit the precision of intermediate results at all, even to the precision of some register. So, the loop isn't guaranteed to terminate and the program has undefined behavior. One solution would be to use `volatile epsilon_f`. – Potatoswatter Feb 06 '15 at 02:18
  • … but `volatile` is only a half-soution, which will in practice prevent a static analyzer from finding an endless loop but won't guarantee any single-precision computation. – Potatoswatter Feb 06 '15 at 02:24
  • @Potatoswatter: I don't know that it is even a half-solution. http://coliru.stacked-crooked.com/a/43ccd29b6d973e12 – rici Feb 06 '15 at 02:29
  • @Potatoswatter: I guess `volatile` will force the variable to be reloaded from memory (which actually only contains 32/64 bits) whenever it's referenced in a computation, but the actual computation might still be performed in a register? (So it is now no longer undefined, and guaranteed to terminate, but probably with a smaller-than-expected value?) – Aasmund Eldhuset Feb 06 '15 at 02:29
  • @Potatoswatter: I correct myself. You can use volatile (although I suppose it is not guaranteed), but it's not epsilon which needs to be volatile; it's the trial sum: http://coliru.stacked-crooked.com/a/b0490ddec84525b7 – rici Feb 06 '15 at 02:50
  • How does assuming there is a precision limit make the program behavior *undefined*? While the C++ standard permits `1 + epsilon > 1` to be evaluated with more precision than the nominal type has, the assignment `epsilon = epsilon / 2` is required to discard any excess precision, and therefore the loop must eventually reach the lower (positive) limit of the floating-point type, after which epsilon becomes zero, if round-to-nearest-ties-to-even is in effect, and then `1 + epsilon > 1` is false even if “infinite precision” is used, and the loop terminates. – Eric Postpischil Dec 24 '20 at 12:08
  • @Potatoswatter: If the arithmetic is performed with round-to-nearest-ties-to-even or round-toward-zero or round-downward, the loop must terminate, and there is no undefined behavior. The assignment `epsilon = epsilon / 2` is required to discard excess precision, so `epsilon` becomes smaller in each iteration until it reaches the lower bound of positive representable numbers, after which it becomes zero, and then `1 + epsilon > 1` is false even if “infinite precision” is used, and the loop terminates. – Eric Postpischil Dec 24 '20 at 19:06