Beginner C Conundrum: why doesn't 25 == 25?

Question

I have a simple task for my first lab in my intermediate C programming class. I'm taking an array of 8 doubles from the user, and then 1 additional double. I then check the square of one double from the array plus the square of other doubles in the array to see if they are equivalent to the square of the last input given to the program (the additional double).

My problem is that for some reason, when my two inputs squared equal the additional input squared, my compiler doesn't think so.

Please let me know what I'm doing wrong here; I'm using Codelite with the gnu gdb debugger and gcc compiler.

Sample input: 4 3 3 3 3 3 3 3 5

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

int main(int argc, char **argv)
{
    int input[8];
    double inputSquared[8];
int pairs[16];
int i, j, k;  //i and j are for loop workers.
double x;
int numPairs = 0;
printf("Welcome to Lab 1.\nEnter 8 integers.  After each entry, press the enter button.\n");
printf("---------------------------------------\n");
for(i=0; i<8; i++){
    printf("Enter integer %d:", i+1);
    scanf("%d", &input[i]);
    printf("\n");
}

//printf("Now enter one more integer.\n  The sum of the squares of the following o this integer squared.\n");
printf("Enter an integer: ");
scanf("%lf", &x);

for(k = 0; k<8; k++){
    inputSquared[k] = pow((double)input[k], 2);
}

for(i = 0; i<8; i++){
    for(j = i + 1; j<8-1; j++){  //does not check for pairs reflexively. If 1 is in the array, it does not check 1^2 + 1^2.
        printf("%lf, %lf; %lf; %lf, %d \n", inputSquared[i], inputSquared[j], pow(x, 2.0), inputSquared[i] + inputSquared[j], ((inputSquared[i] + inputSquared[j]) == ((pow(x, 2.0)))));
        if(inputSquared[i] + inputSquared[j] == pow(x, 2.0)){
            pairs[2 * numPairs] = input[i];
            pairs[2 * numPairs + 1] = input[j];
            numPairs++;
        }
    }
}
if(numPairs == 1)
    printf("\nYou have %d pair:", numPairs);  // grammar condition for having 1 pair
else
    printf("\nYou have %d pairs:\n", numPairs);

for(i = 0; i < numPairs; i++)
    printf("(%d,%d)", pairs[2 * i], pairs[2 * i + 1]);

scanf("%lf", &x);

return 0;
}

Answer 1

When you're using double you are invoking what is called "Floating Point Arithmetics". In practice, this is an imprecise way to represent numbers, simply because not every real number can be precisely represented in a finite number of bits.

As a result, some numbers, for instance, can never be represented precisely, and instead of "25" your computer (not the compiler), has a result of, say, 25.00000000000071 - close enough for all practical purposes, but clearly not equal to 25.

There are few other quirks related to IEEE-754 floating point - adding 0.1 ten times to 0.0 doesn't necessarily reach 1.0, etc.

Which is why, in practice, you should never compare for equality, but you should compare whether the absolute value of the difference between two numbers is smaller than an (small, but large enough) epsilon.

const double espilon = 0.00000001
if (fabs(x - y) < epsilon) {
    // Assume approximately equal
}

A good (MATLAB based, but the concepts apply throughout everything that uses IEEE-754, nowadays, that means about everywhere) introduction to what is going on is here:

https://www2.bc.edu/~quillen/sp11/mt414/pres/floatarith.pdf

Also, pow() is not the most reasonable way to do squares - pow(x, m) internally computes exp(m log(x)) in order to handle arbitrary real-number bases and exponents. This particular conversion is exactly why it looses precision.

This also means pow() is ridiculously slow when compared with multiplication-based approach.

Answer 2

If you computed the square of x as:

x * x

or even

(double)x * (double)x

then you would get an exact square. In other words,

4 * 4 + 3 * 3 == 5 * 5             => 1 (true)
4.0 * 4.0 + 3.0 * 3.0 == 5.0 * 5.0 => 1 (true)

In fact,

5 * 5 == 5.0 * 5.0                 => 1 (true)

But,

5 * 5 == pow(5.0, 2.0)             => 0 (false)

because the math library (not the compiler) does not bother checking to see if the second argument of pow happens to be a small integer. It just goes ahead and works out the value by computing pow(x, p) = ⅇp*ln(x). Unfortunately, that value is not, in general, representable as a double -- in fact, it is usually not even a rational number, so it doesn't really have a finite representation at all -- and the math library settles for a reasonably close approximation computed by a Taylor series (or something similar).

So not only is pow(5.0, 2.0) slightly inaccurate, it is also quite complicated to compute. Whereas 5 * 5 is just a single machine instruction. You can probably draw your own conclusion.

---

Many people will tell you that you should "never" compare floating point values for equality. That advice tends to lead to (or come from) a mental model where a floating point number is a kind of fuzzy, out-of-focus thing which can fluctuate unpredictably, almost as though it were subject to Heisenberg uncertainty. This is not actually a very good way of thinking about floating point numbers. Every floating point number is a precise representation of some rational value; in fact, it is a precise representation of a rational value whose denominator is a power of 2 and whose numerator is 0 or an integer between 1 and 2^k-1 for some smallish k (52 for doubles on most modern CPUs).

Consequently, it is quite predictable under which circumstances floating point arithmetic will be precise. For example, if x is known to be an integer whose absolute value is less than the number of seconds in a year, then (double)x * x is precisely correct. On the other hand, fractions whose denominators (when reduced to lowest terms) are odd (or have an odd factor) can never be represented precisely as floating point numbers.

Answer 3

You should never compare double values for exactness.

Use an epsilon (here an epsilon value of 0.00000001 is used):

if (abs(x - y) < 0.00000001)
{
}

so,

if ( abs(inputSquared[i] + inputSquared[j] - pow(x, 2.0)) < 0.0000001)
{ 

[Bbtw: there are a quite a few similar questions on SO]

Answer 4

The core problem here is that you have a low-quality math library in which pow fails to return good results. You can work around this problem in this particular case by using x*x in place of pow(x, 2).

Although floating-point cannot represent all numbers exactly (and neither can integer arithmetic, of course) and some mathematical routines are difficult to implement accuracy, a good math library strives to return accurate results.

For the pow function, there are a number of cases in which accuracy should be especially sought, including those for which the result is exactly representable. This includes cases such as pow(x, 2), where x*x is exactly representable in the floating-point format.

The workaround above will fail for values of x that have enough significant bits that x*x cannot be represented exactly in the floating-point format (or when you attempt to read values from input that are not exactly representable). In these cases, other answers have advised you to compare floating-point numbers using a tolerance. However, this introduces additional bugs. It reduces false negatives (rejecting numbers as unequal even though they would be equal if computed exactly) at the expense of increasing false positives (accepting numbers as equal even though they are truly unequal).