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?

@Hans Passant This is the first time I have heard something like this. What's the issue? Shouldn't the compiler warn me about this kind of stuff? — Byzantian, Mar 11 '13 at 21:00
@HansPassant, never say never. It's OK if the double was assigned *from* an int. In this case though, it wasn't. — Mark Ransom, Mar 11 '13 at 21:03

score 3 · Accepted Answer · 2013-03-12T05:30:41.520

3

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).

edited Mar 12 '13 at 05:30

answered Mar 11 '13 at 21:00

Makes sense so far. But why does my implementation which STILL returns double seem to work? – Byzantian Mar 11 '13 at 21:03
3

@cyberpunk_, your version works because all the doubles you process were converted from integers. The `pow` math library function uses faster math but it requires actual floating point operations that lose a tiny bit of accuracy. – Mark Ransom Mar 11 '13 at 21:05
One way to implement pow(x,y) is exp(y * ln(x)) and both ln() and exp() can cause round off errors. (I am sure the real pow() does not use this method.) – brian beuning Mar 12 '13 at 00:11
Isn't there an underscore missing? – Alexey Frunze Mar 12 '13 at 04:27
@AlexeyFrunze There is, indeed. – Mar 12 '13 at 05:30
How about just adding 0.5 before casting to int? That would actually round the float to the closes int, which is most likely what pow meant by its return value. – notadam Dec 31 '15 at 04:10

score 1 · Answer 2 · answered Mar 11 '13 at 21:01

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.

math.h pow() function not working correctly?

2 Answers2