#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?
24
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−23 in all C implementations.
6
@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
floatshould be enough (but its not for OP). -
Saving to a
floatvariable and then reading it should force thefloatmath. -
Using
volatileto 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.
2
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.