Float data type uncertainty

Question

I am doing a numerical analysis of a math software I developed. I want to identify what is the uncertainty of my result. Being f() my method and x an input value, I want to identify y of my result as f(x) +/- y. My f() method has multiple operations between float variables. To study the error propagation occurred in f(), I have to apply the Statistical Propagation of Uncertainty formulas and in order to do so I have to know the uncertainty of a float variable.

I do understand the architecture of a float variable as specified in the IEEE 754 standard and the rounding error converting a decimal value to float inherent to the latter.

From what I understood of the literature, the FLT_EPSILON macro in http://www.cplusplus.com/reference/cfloat/ defines my y value but this quick test proves it wrong:

float f1 = 1.234567f;
float f2 = 1.234567f + 1.192092896e-7f;
float f3 = 1.234567f + 1.192092895e-7f;

printf("Inicial:\t%f\n", f1);
printf("Inicial:\t%f\n", f2);
printf("Inicial:\t%f\n\n", f3);

Output:

Inicial:  1.234567
Inicial:  1.234567
Inicial:  1.234567

When the expected output should be:

Inicial:  1.234567
Inicial:  1.234568 <---
Inicial:  1.234567

What is that I am wrong about? Should not the float value of x + FLT_EPSILON and x - FLT_EPSILON be the same?

EDIT: My question is being R the float value of x, what is the y value that x + y || x - y equals the same R float value?

Answer 1

Propagation of uncertainty is from the field of statistics and refers to how uncertainties in inputs affect mathematical functions of them. The analysis of errors that occur in computational arithmetic is numerical analysis.

FLT_EPSILON is not a measure of uncertainty or error in floating-point results. It is the distance between 1 and the next value representable in the float type. Hence, it is the size of steps between representable numbers at the magnitude of 1.

When you convert a decimal numeral to floating-point, the rounding error that results may have a magnitude of up to ½ the step size when the common round-to-nearest mode is used. The reason the bound is ½ the step size is that for any number x (within the finite domain of the floating-point format), there is a representable value within ½ the step size (inclusive). This is because, if there is a representable number more than ½ the step size in one direction, there is a representable number less than ½ the step size in the other direction.

The step size varies with the magnitudes of the numbers. With binary floating-point, it doubles at 2, and again at 4, then 8, and so on. Below 1, it halves, and again at ½, ¼, and so on.

When you perform floating-point arithmetic operations, the rounding that occurs in the computation may compound or cancel previous errors. There is no general formula for the final error.

The two numerals use used in your sample code, 1.192092897e-7f and 1.192092896e-7f, are so close together that they convert to the same float value, 2⁻²³. That is why there is no difference in your f2 and f3.

There is a difference between f1 and f2, but you did not print enough digits to display it.

You ask “Should not the float value of x + FLT_EPSILON and x - FLT_EPSILON be the same?”, but your code does not contain x - FLT_EPSILON.

Re: “My question is being R the float value of x, what is the y value that x + y || x - y equals the same R float value?” This is trivially satisfied by y = 0. Did you mean to ask what is the largest value of y that satisfies the condition? That is a bit complicated.

The step size for a number x is called the ULP of x, which we may consider as a function ULP(x). ULP stands for Unit of Least Precision. It is the place value of the least digit in the floating-point representation of x. It is not a constant; it is a function of x.

For most values representable in a floating-point format, the largest y that satisfies your condition is ½ ULP(x) of the least digit in the floating-point representation of x is even and, if the digit is odd, it is just under ½ ULP(x). This complication arises from the rule that the results of arithmetic are rounded to the nearest representable value and, in case of a tie, the value with the even low digit is chosen. Thus, adding ½ ULP(x) to x will yield a tie that will round to x if the low digit is even, but will not round to x if the low digit is odd.

However, for x that are on the boundary where the ULP changes, the largest y that satisfies your condition is ¼ ULP(x). This is because, just below x (in magnitude), the step size changes, and the next number lower than x is half of x’s step size away instead of the usual full step size. So you can only go halfway toward that value before changing the result of the subtraction, so the most y can be is ¼ ULP(x).

Answer 2

float f1 = 1.234567f;
float f2 = f1 + 1.192092897e-7f;
float f3 = f1 + 1.192092896e-7f;

printf("Inicial:\t%.20f\n", f1);
printf("Inicial:\t%.20f\n", f2);
printf("Inicial:\t%.20f\n\n", f3);

Output:

Inicial:        1.23456704616546630000
Inicial:        1.23456716537475590000
Inicial:        1.23456716537475590000

No, your expectation is wrong
In the first printf call, you're printing the variable f1 with no effect which is just 1.234567f.

Answer 3

Float is a 32 bit IEEE 754 single precision Floating Point Number: 1 bit for the sign, 8 bits for the exponent, and 23* for the value, i.e. float has 7 decimal digits of precision.

Increase the printf number of printed digits to see more but after 7 digits its just noise:

#include <stdio.h>

int main(void) {

 float f1 = 1.234567f;
 float f2 = 1.234567f + 1.192092897e-7f;
 float f3 = 1.234567f + 1.192092896e-7f;

 printf("Inicial:\t%.16f\n", f1);
 printf("Inicial:\t%.16f\n", f2);
 printf("Inicial:\t%.16f\n\n", f3);

 return 0;
}

Output:

Inicial:        1.2345670461654663                                                                                                           
Inicial:        1.2345671653747559                                                                                                           
Inicial:        1.2345671653747559