Issue with pow and round - answers are not equivalent

Question

I'm having an issue creating a function that checks if a root can be simplified. In this example, I'm trying to simplify the cube root of 108, and the first number that this should work for is 27.

In order to do this, I am calling pow() with the number being the index (in this case, 27), and the power being (1/power), which in this instance is 3. I then compare that to the rounded answer of pow(index,(1/power)), which should also be 3.

Included is a picture of my problem, but basically, I am getting two answers that are equivalent to 3, yet my program is not recognizing them as equal. It seems to be working elsewhere in my program, but will not work here. Any suggestions as to why?

int inside = insideVal;
int currentIndex = index;
int coeff = co;
double insideDbl = pow(index, (1/(double)power));
double indexDbl = round(pow(index,(1/(double)power)));
cout<<insideDbl<< " " << indexDbl <<endl;
//double newPow = (1/(double)power);
vector<int> storedInts = storeNum;
if(insideDbl == indexDbl){
    if(inside % currentIndex == 0){
        storedInts.push_back(currentIndex);
        return rootNumerator(inside/currentIndex, currentIndex, coeff, power, storedInts);
    }
    else{
        return rootNumerator(inside, currentIndex + 1, coeff, power, storedInts);
    }
}
else if(currentIndex < inside){
    return rootNumerator(inside, currentIndex + 1, coeff, power, storedInts);
}

I tried to add a picture, but my reputation apparently wasn't high enough. In my console, I am getting "3 3" for the line that reads cout<<insideDbl<< " " << indexDbl <<endl;

EDIT:

Alright, so if the answers aren't exact, why does the same type of code work elsewhere in my program? Taking the 4th Root of 16 (which should equal 2) works using this segment of code:

    else if( pow(initialNumber, (1/initialPower)) == round(pow(initialNumber,(1/initialPower)))){
        int simplifiedNum = pow(initialNumber, (1/initialPower));
        cout<<simplifiedNum;
        Value* simplifiedVal = new RationalNumber(simplifiedNum);
        return simplifiedVal;
    }

despite the fact that the conditions are exactly the same as the ones that I'm having trouble with.

Answer 1

Well you are a victim of finite precision floating point arithmetic.

What happened?

This if(insideDbl == indexDbl), is very dangerous and misleading. It is in fact a question whether (Note: I made up the exact numbers but I can give you precise ones) 3.00000000000001255 is the same as 2.999999999999996234. I put 14 0s and 14 9s. So technically the difference goes beyond 15 most significant places. This is important.

Now if you write insideDbl == indexDbl, the compiler compares the binary representantions of them. Which are clearly different. However, when you simply print them, the default precision is like 5 or 6 significant digits, so they get rounded, and seem to be the same.

How to check it?

Try printing them with:

typedef std::numeric_limits< double > dbl_limits;
cout.precision(dbl::max_digits10);
cout << "Does " << insideDbl  << " == " << indexDbl << "?\n";

This will set the precision, to the number of digits, the are necessary to differentiate two numbers. Please note that this is higher than the guaranteed precision of computation! That is the root of confusion.

I would also encourage reading numeric_limits. Especially about digits10, and max_digits10.

Why sometimes it works?

Because sometimes two algorithms will end up using the same binary representation for the final results, and sometimes they won't.

Also 2 can be a special case, as I believe it can be actually represented exactly in binary form. I think (but won't put my head on it.) all powers of 2 (and their sums) can be, like 0,675 = 0,5+0,125 = 2^-1 + 2^-3. But please don't take it for granted unless someone else confirms it.

What can you do?

Stick to the precise computations. Using integers, or whatever. Or you could assume that everything 3.0 +/- 10^-10 is actually 3.0 (epsilon comparisons), which is very risky, to say the least, when you do care about precise math.

Tl;dr: You can never compare two floats or doubles for equality, even when mathematically you can prove the mentioned equality, because of the finite precision of computations. That is, unless you are actually interested in the same binary representation of the value, as opposed to the value itself. Sometimes this is the case.

Answer 2

I suspect that you'll do better by computing the prime factorisation of insideVal and taking the product of those primes that appear in a multiple of the root.

For example

108 = 2² × 3³

and hence

³√108 = 3 × ³√2²

and

324 = 2² × 3⁴

and hence

³√324 = 3 × ³√(2² × 3)

You can use trial division to construct the factorisation.

Edit A C++ implementation

First we need an integer overload for pow

unsigned long
pow(unsigned long x, unsigned long n)
{
    unsigned long p = 1;

    while(n!=0)
    {
        if(n%2!=0) p *= x;
        n /= 2;
        x *= x;
    }
    return p;
}

Note that this is simply the peasant algorithm applied to powers.
Next we need to compute the prime numbers in sequence

unsigned long
next_prime(const std::vector<unsigned long> &primes)
{
    if(primes.empty()) return 2;

    unsigned long p = primes.back();
    unsigned long i;

    do
    {
        ++p;
        i = 0;
        while(i!=primes.size() && primes[i]*primes[i]<=p && p%primes[i]!=0) ++i;
    }
    while(i!=primes.size() && primes[i]*primes[i]<=p);

    return p;
}

Note that primes is expected to contain all of the prime numbers less than the one we're trying to find and that we can quit checking once we reach a prime greater than the square root of the candidate p since that could not possibly be a factor.
Using these functions, we can calculate the factor that we can take outside the root with

unsigned long
factor(unsigned long x, unsigned long n)
{
    unsigned long f = 1;

    std::vector<unsigned long> primes;
    unsigned long p = next_prime(primes);

    while(pow(p, n)<=x)
    {
        unsigned long i = 0;
        while(x%p==0)
        {
            ++i;
            x /= p;
        }
        f *= pow(p, (i/n));

        primes.push_back(p);
        p = next_prime(primes);
    }
    return f;
}

Applying this to your example

std::cout << factor(108, 3) << std::endl; //output: 3

gives the expected result. For another example, try

std::cout << factor(3333960000UL, 4) << std::endl; //output: 30

which you can confirm is correct by noting that

3333960000 = 30⁴ × 4116

and checking that 4116 doesn't have any factor that is a power of 4.