How can I solve Float type Machine Epsilon Calculation Error?

#include <stdio.h>
void get_eps(float *eps) {
    *eps = 1.0f;

    while ((1.0f + *eps / 2.0f) != 1.0f) {
        *eps /= 2.0f;
    }
}

int main(){
    float eps;

    get_eps(&eps);
    printf("Machine Accuracy (get_eps): t%.10gn", eps);

    return 0;
}

I am trying to calculate the machine epsilon of float(32bit).

Expecting 1.19e-07, which is the known machine epsilon for the float type, I get 1.08e-19, which is the machine epsilon for long double. Please help me with this issue. Is there any points that I am missing?

The C standard permits an implementation to evaluate floating-point expressions with more precision than their nominal type.

However, it requires this excess precision be “discarded” by assignments and casts. So changing the test to (float) (1.0f + *eps / 2.0f) != 1.0f should work.

Supplement

Note that the specific properties of the float type are not specified by the C standard, so the so-called “epsilon” is not required to be 2⁻²³ in all C implementations.

@Eric Postpischil good answer should work.

Since it does not, OP’s non-compliant compiler leaves no way to certainly make it work.

As C allows intermediate results to use wider math.
printf("%dn", FLT_EVAL_METHOD); to see what is being used. Any output other than 0 implies wider math possible.

All approaches comes down to forcing the intermediate result to a float.

Casting to a float should be enough (but its not for OP).
Saving to a float variable and then reading it should force the float math.
Using volatile to defeat certain optimizations may help.

Possible alternative code.

void get_eps(float *eps) {
  *eps = 1.0f;

  while (1) {
    volatile float sum = 1.0f + *eps / 2.0f;
    if (sum == 1.0f) {
      return;  // or break
    }
    *eps /= 2.0f;
  }
}

I added a bit more info to teapot418’s example and got the following:

sizeof(float) == 4
CHAR_BIT == 8
FLT_EPSILON == 1.192093e-07

So confirmed float is 32-bit, but OP’s program is not outputting FLT_EPSILON as the accuracy. As Eric P. alluded to, the compiler is likely converting everything to long double for all computations and using that.

As verification, this code outputs: 1.157921e+77 inf, where it should output inf inf so it’s definitely promoting to a larger type. That it’s ignoring the (float) cast probably qualifies as a significant bug.

    eps = FLT_MAX;
    printf("%e ", (float)(eps * eps));
    eps *= FLT_MAX;
    printf("%en", eps);

Using the alternative code provided by Chux does give the desired answer, even without using volatile, so that’s a viable workaround. For me, I would just #include <float.h> and use FLT_EPSILON as given for the tolerance. That would be the best chance of ensuring portability across compilers and platforms.

I suspect that you are using the MS compiler or another (gcc?) that has its default behaviour to generate legacy x87 code by default and a high optimisation level so that all intermediate results stay on the stack.

It all behaves as expected in debug mode or with optimisation disabled.

However, I cannot quite reproduce your results on MSVC 17.1. Precision is way too good – consistent with double precision 64 bit arithmetic (which has been hardwired into MSVC since v6.0 when they withdrew 80-bit FP support 🙁 ).

The interesting question here is which compiler are you using that has legacy 80bit x87 FP code generation enabled by default. The hardware support on the x87 is for long double 80 bit (64 bit mantissa) by default (some most? compilers set it to 64 bit and 53 bit mantissa – I know MS does).

The sample code when run with default settings for FP code generation (not set), FP:Precise or FP:Strict behave correctly here and in debug mode when temp results are stored and loaded from memory but it misbehaves when the combination of settings is FP legacy x87 code generation and optimised release code generation so that all the intermediate results remain on the x87 stack.

However, under those conditions the MS compiler has default accuracy 64 bit reals with a 53bit mantissa.

This is the generated code that causes this behaviour:

                           get_eps(&eps);
    >00B6186B DD D9                fstp        st(1)  
     00B6186D D9 C9                fxch        st(1)  
     00B6186F D9 C9                fxch        st(1)  
     00B61871 D9 5C 24 30          fstp        dword ptr [esp+30h]  
     00B61875 D9 44 24 30          fld         dword ptr [esp+30h]  
     00B61879 D9 C0                fld         st(0)  
     00B6187B D9 C9                fxch        st(1)  
     00B6187D D8 CB                fmul        st,st(3)  
     00B6187F D9 C0                fld         st(0)  
     00B61881 D8 C3                fadd        st,st(3)  
     00B61883 DD EB                fucomp      st(3)  
     00B61885 DF E0                fnstsw      ax  
     00B61887 F6 C4 44             test        ah,44h  
     00B6188A 7A DF                jp          main+1Bh (0B6186Bh)  
     00B6188C DD D8                fstp        st(0)  
                            printf("Machine Accuracy (get_eps): t%.10gn", eps);

And this is the content of the stack on exit from the loop prior to the fstp.

   ST0 = +1.1102230246251565e-0016   ST1 = +2.2204460492503130e-0016   
   ST2 = +1.0000000000000000e+0000   ST3 = +5.0000000000000000e-0001   
   ST4 = +0.0000000000000000e+0000   ST5 = +0.0000000000000000e+0000   
   ST6 = +0.0000000000000000e+0000   ST7 = +1.0000000000000000e+0000   
   CTRL = 027F STAT = 6020 TAGS = 00FF EIP = 00B61883 EDO = 00000000

Output is:

Machine Accuracy (get_eps):     2.220446049e-16

This is with FP:strict (which is the one that I expected to work OK). I can no longer reproduce this result (although I’m convinced from my notes that was what I observed the first time I tested it).

Things get really interesting when MSVC is allowed either FP:precise or the optimisation freedoms of fastmath when the code becomes a lot tighter and arguably very wrong.

Machine Accuracy (get_eps):     0

The generated code in this case is:

                            get_eps(&eps);
0064184B DD DA                fstp        st(2)  
0064184D D9 C1                fld         st(1)  
0064184F D8 CB                fmul        st,st(3)  
00641851 DD E1                fucom       st(1)  
00641853 DF E0                fnstsw      ax  
00641855 F6 C4 44             test        ah,44h  
00641858 7A F1                jp          main+1Bh (064184Bh)  
0064185A DD D8                fstp        st(0)  
                            printf("Machine Accuracy (get_eps): t%.10gn", eps);

Intel compiler gets the right answer but does it by forcing clunky memory stores of the intermediate results. However, this only seems to be true if the project has default FP (unset).

A quick test showed that on Intel fp:Precise gives the expected answer

Machine Accuracy (get_eps):     1.192092896e-07

Strict gives: Machine Accuracy (get_eps): 2.220446049e-16

Fast gives: Machine Accuracy (get_eps): 1.192092896e-07
(but it seems to have cheated by loading a constant computed at compile time)

The suggested alteration (float) (1.0f + *eps / 2.0f) != 1.0f in the previous answer does not affect the code generation in the failing cases. The temporary results remain on the x87 stack with all the spurious extra precision that goes with it.

Filed under: Kiến thức lập trình - @ 14:05

Thẻ: cfloating-pointfloating-accuracyepsilon

Thiết kế website giá rẻ

Danh mục

How can I solve Float type Machine Epsilon Calculation Error?

Supplement