Floating point issue in C when solving a math expression

Viewed 179

I have a problem solving math function in C.

double task1_double() {
    double a = 1000;
    double b = 0.0001;
    double result = (pow((a + b), 2) - (pow(a, 2) + 2 * a * b))/(pow(b, 2));
    return result; }

float task1_float() {
    float a = 1000;
    float b = 0.0001f;
    float result = (powf((a + b), 2) - (powf(a, 2) + 2 * a * b))/(powf(b, 2));
    return result; }

When I use a double data type it returns 1.001172, though with float data type the result is 6250000.000000.

Can somebody explain why this happens? Thanks in advance

2 Answers

Your numerator is exactly

1000000.20000001 - 1000000.2 = 0.00000001

but when you subtract two large numbers which are almost equal, the relative rounding error in the result can blow up. This is what you are seeing. It is a consequence of the fact that float numbers have about 7 decimal digits of precision, whereas double numbers have about 16 decimal digits.

Let's do it step by step:

                           exact            float        double
x = pow((a + b), 2)        1000000.20000001 1000000.25   1000000.200000009965
y = pow(a, 2) + 2 * a * b  1000000.2        1000000.1875 1000000.199999999953
x - y                      0.00000001       0.0625       0.000000010011717677

For such a small difference between the two numbers relative to their magnitude, you would usually get a float result equal to 0.0. But in this case it just happens that 1000000.20000001 and 1000000.2 lie either side of a rounding boundary, causing the former to be rounded up and the latter to be rounded down. So their difference is out by six orders of magnitude.

By outputting intermediate results, it can be found that it is caused by the loss of precision of float


#include <stdio.h>
#include <math.h>

double task1_double() {
    double a = 1000;
    double b = 0.0001;
    double s1, s2, s3;
    s1 = pow((a + b), 2);
    s2 = pow(a, 2) + 2 * a * b;
    s3 = pow(b, 2);
    double result = (pow((a + b), 2) - (pow(a, 2) + 2 * a * b)) / (pow(b, 2));

    printf("_double : %lf - %lf = %.15lf\n", s1, s2, s1 - s2);
    printf("_double : (%lf - %lf) / %.10lf = %lf\n", s1, s2, s3, result);

    return result;
}

float task1_float() {
    float a = 1000;
    float b = 0.0001f;
    float s1, s2, s3;
    s1 = powf((a + b), 2);
    s2 = powf(a, 2) + 2 * a * b;
    s3 = powf(b, 2);
    float result = (powf((a + b), 2) - (powf(a, 2) + 2 * a * b)) / (powf(b, 2));

    printf("_float : %lf - %lf = %.15lf\n", s1, s2, s1 - s2);
    printf("_float : (%lf - %lf) / %.10lf = %lf\n", s1, s2, s3, result);

    return result;
}

int main()
{
    printf("%.10lf\n%.10lf\n", task1_double(), task1_float());
    return 0;
}

Output:

_float : 1000000.250000 - 1000000.187500 = 0.062500000000000
_float : (1000000.250000 - 1000000.187500) / 0.0000000100 = 6250000.500000
_double : 1000000.200000 - 1000000.200000 = 0.000000010011718
_double : (1000000.200000 - 1000000.200000) / 0.0000000100 = 1.001172
1.0011717677
6250000.5000000000

As can be seen: After the float is operated, a slight error occurs, but dividing by the small value of b will cause the error to be magnified many times

Related