How can I solve Float type Machine Epsilon Calculation Error?

Question

#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%.10g\n", 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?

Answer 1

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.

Answer 2

@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("%d\n", 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;
  }
}