math.h pow() function not working correctly?

Question

I was working on a personal project in which I needed to determine all prime powers between 0 and 999. Because I am not particularly good at math, I tired the following ham-fisted bruteforce approach.

bool constexpr is_prime(int imp)
{
    return imp == 1 ? false : (imp % imp == 0 && imp % 1 == 0 && [&imp]{ for(int i = 2; i < imp; ++i)  if(imp % i == 0) return false; return true;}());
}

bool is_prime_power(int imp)
{
    for(int i = 1; i < 1000; ++i)
        if (is_prime(i))
            for (int j = 0; j < 100; ++j)
                if (imp == pow(i, j))
                    return true;
    return false;
}

For 0...30 the output should be (according to A000961):

1 2 3 4 5 7 8 9 11 13 16 17 19

However, this is what I get instead:

1 2 3 4 5 7 8 9 11 16 19

Where did 13 and 17 vanish to?

Since I could not find any logical problems with my approach I implemented my own pow() function.

double constexpr _pow(double base, double exp)
{
    return exp == 0 ? 1 : base*pow(base, exp - 1);
}

Now if I call my version of _pow() instead of pow() from math.h the output is displayed as excepted. Is my implementation wrong? If not, pow() from math.h can't be working correctly. Any idea what is causing this?

Answer 1

The problem is that double (and floating-point numbers in general) aren't exact and the math functions in the standard library also use approximate formulae to do the calculations. If you call pow(2, 10), you may get 1023.999375, for example (which is then, depending on the context, may be truncated down to 1023).

So, instead of using the floating-point pow() function from the standard library, you can go with your own, exact implementation for integers:

int constexpr _pow(int base, int exp)
{
    return exp == 0 ? 1 : base * _pow(base, exp - 1);
}

(or change it to whatever integral type you want if you need unsigned or long long).

Answer 2

You never need to verify that a number divides itself and that one divides a number. This is useless: imp % imp == 0 && imp % 1 == 0

Otherwise the problem you have is because pow returns a double or float while you compare it to an integer. The comaparison may fail because of the rounding error. I would suggest you implement your own integer power operation so that you avoid any rounding error. Alternatively convert i to double and use comparison with some small epsylon tolerance.