C++ program of (a) power (b) without using builtin function

In a previous post I've criticized a careless, effortless but highly-valued post pretending to contain an algorithm, written in C++, to compute a^b, when indeed it computes the wrong results, since it implements (likely wrongly) a wrong algorithm, hence teaching something wrong.

I've written the algorithm and analysed his code a bit, but I haven't given an actual implementation. Here it is in this post.

First of all, when you provide a functionality, e.g. a piece of code which compute a^b, you put it into a function: it's its natural “home”. Then you can use the function to show it works. This is testing. I'll show you later; now, let's implement our pow function.

Beware! Achtung!

Since C11 was released, I use at least C11 (unless there's legacy code to be maintained).

Let's start…

Let's start with the prototype.

double my_pow(double a, int b);

Exactly. Our power of function can compute a^b when a is a real number (a double in C/C++), but b must be integer. Real exponents are a different story and we leave them alone.

However, even if a were an integer, the result must still be a double, because when b is negative in fact we have a fraction, 1/a^|b|, so the return value of the function must be a double.

Implementation:

#include 

double my_pow(double a, int b)
{
    double result = 1.0;
    bool invert = false;
    if (b < 0) {
        invert = true;
        b = -b;
    }
    for (int n = 0; n < b; ++n) {
        result = result * a;
    }
    if (invert) {
        if (result != 0) {
            result = 1.0/result;
        } else {
            result = std::numeric_limits<double>::infinity();
        }
    }
    return result;
}

You can see that when I must invert result, but this is 0, I put into result a special value, which is the hugest value a double can represent, and can be read a infinity — of course it isn't actually infinity, though. This is provided by the C++ standard. Another option is to throw an exception (then the include stdexcept is needed):

    if (invert) {
        if (result != 0) {
            result = 1.0/result;
        } else {
            throw std::overflow_error("division by zero");
        }
    }

This code throws the standard exception overflow error, with text “division by zero”, because a divide-by-zero exception doesn't exist, and when you do 1/0, you can say it gives infinity, which is a value which cannot fit inside a double, hence the overflow choice.

The exception must then be handled… I prefer the other solution, where result becomes a special value representing “infinity” for a double. In real use cases the choice depends on the application: is it ok the computation goes on with such a “huge” value, or should be thrown an exception because there's a problem?

Testing

We have our function and we want to know if it gives correct results. I hate to check things by hand: let's the computer do it for us!

In order to do so, we must use the standard pow function and compare results. Comparing floating point numbers is tricky. Here I am happy if my poor my_pow() is able to compute a value which is close enough to the value given by pow. How close? We can tune our program changing a constant. In this example the result of my_pow() must be between r - 0.0001 and r + 0.0001, where r is the result given by pow().

The function to check if the values are close enough:

bool close_enough(double a1, double a2, double threshold)
{
    return a1 > (a2 - threshold) && a1 < (a2 + threshold);
}

I want also a function which, given the inputs, computes a^b with my_pow(), then with pow() and at last it outputs the verdict, handling also a little bit of formatting, just to have a “nice” output.

Here it is:

void test_pow(double a, int b)
{
    constexpr double thrs = 0.0001;
    constexpr int w = 10;
    
    double r1 = my_pow(a, b);
    double r2 = pow(a, b);
    std::cout << std::setw(w) << a
              << " ** "
              << std::setw(w) << std::left << b
              << " = "
              << std::right << std::setw(w)
              << std::setprecision(4) << r1
              << " ... "
              << (close_enough(r1, r2, thrs) ? "ok!" : "NO.")
              << std::endl;
}

Now we are ready for the main(). Let's test the input values I've used in my previous post. This is a very boring main, …

int main()
{
    test_pow(0, 0);
    test_pow(1, 10);
    test_pow(2, 10);
    test_pow(2, -1);
    test_pow(2, -2);
    test_pow(2, -3);
    test_pow(-2, 1);
    test_pow(4, 1);
    test_pow(-2, 2);
    test_pow(-2, 3);
    test_pow(-2, -1);
    test_pow(-2, -2);
    test_pow(-2, -3);
    test_pow(10, -1);
    test_pow(11, -1);
    test_pow(12, -1);
    test_pow(10, -2);
    test_pow(12, -2);
    test_pow(10, -3);
    test_pow(100, -1);
    test_pow(100, -2);
    test_pow(100, -3);
    test_pow(100, -4);
    test_pow(-100, -4);
    test_pow(1000, 1);
    return 0;
}

Compiled and run, the output is:

         0 ** 0          =          1 ... ok!
         1 ** 10         =          1 ... ok!
         2 ** 10         =       1024 ... ok!
         2 ** -1         =        0.5 ... ok!
         2 ** -2         =       0.25 ... ok!
         2 ** -3         =      0.125 ... ok!
        -2 ** 1          =         -2 ... ok!
         4 ** 1          =          4 ... ok!
        -2 ** 2          =          4 ... ok!
        -2 ** 3          =         -8 ... ok!
        -2 ** -1         =       -0.5 ... ok!
        -2 ** -2         =       0.25 ... ok!
        -2 ** -3         =     -0.125 ... ok!
        10 ** -1         =        0.1 ... ok!
        11 ** -1         =    0.09091 ... ok!
        12 ** -1         =    0.08333 ... ok!
        10 ** -2         =       0.01 ... ok!
        12 ** -2         =   0.006944 ... ok!
        10 ** -3         =      0.001 ... ok!
       100 ** -1         =       0.01 ... ok!
       100 ** -2         =     0.0001 ... ok!
       100 ** -3         =      1e-06 ... ok!
       100 ** -4         =      1e-08 ... ok!
      -100 ** -4         =      1e-08 ... ok!
      1000 ** 1          =       1000 ... ok!

Let's try with 0 and a negative exponent: we add the line

    test_pow(0, -2);

and the last line of the output now is:

         0 ** -2         =        inf ... NO.

This doesn't mean that pow(0, -2) is something different from inf. The fact is that the special double “huge” value has the property that it can only match exactly. If you look at close_enough(), you see we have excluded the lower and upper bound (< and > instead of <= and >=). But the math of inf is this:

inf + a_finite_value = inf
inf - a_finite_value = inf

Therefore close_enough() return false (inf > inf && inf < inf is clearly false). If you change the code including the extremes of the range, like this:

bool close_enough(double a1, double a2, double threshold)
{
    return a1 >= (a2 - threshold) && a1 <= (a2 + threshold);
}

The the last line of the output will be:

        0 ** -2         =        inf ... ok!

Of course if we did the choice of throwing an exception, the last line wouldn't appear at all. Glad to have chosen the same way pow() does! (Which isn't, indeed, by chance, also considering that it's a plain C function, and C hasn't exceptions, thus…)

The whole code, with a small modification to print the result of pow() in case it doesn't match with my_pow(), can be found in this pastebin and also in this hastebin, whatever you prefer — hastebin is without ads and all the noise pastebin has.