6

I need a version of pow for integers. I have 2 problems that I need to solve with pow:

  1. If the result is larger than my integral type I need to clamp to numeric_limits::max()
  2. I need to be able to deal with 41.99999 rounding up to 42, rather than down to 41

Does C++ provide me some kind of inline solution here, or am I stuck writing my own function:

template <typename T>
enable_if_t<is_integral_v<T>, T> mypow(const T base, unsigned int exp) {
    T result = exp == 0U ? base : 1;

    while(exp-- > 1U) {
        if(numeric_limits<T>::max() / result <= base) return numeric_limits<T>::max();
        result *= base;
    }
    return result;
}
Jonathan Mee
  • 37,899
  • 23
  • 129
  • 288
  • Can't you use [`std::pow`](http://en.cppreference.com/w/cpp/numeric/math/pow) and truncate the result if needed? – NathanOliver Jun 15 '17 at 15:10
  • @NathanOliver Yeah... if there was a way to guarantee **1** and **2** wouldn't be violated by that. Is there a way to accomplish that? – Jonathan Mee Jun 15 '17 at 15:19
  • You can use tail recursion and some numeric trickery to improve runtime performance. Further, you can make this method constexpr. https://stackoverflow.com/questions/9348933/using-recursion-to-raise-a-base-to-its-exponent-c – Andrew Jun 15 '17 at 19:25

2 Answers2

3

Does C++ provide me some kind of inline solution here

No, there is no integer pow in the standard library.

or am I stuck writing my own function

Yes, you get to write your own function. Note that your shown multiplication loop may be slower than using std::pow to implement the function, especially since you also have a branch and division in the loop:

template<class I>
I int_pow_no_overflow(I base, I exp)
{
    double max = std::numeric_limits<I>::max();
    double result = std::round(std::pow(base, exp));
    return result >= max
        ? max
        : result;
}

For a more generic approach, you may want to consider underflow as well.

There are as well other, faster algorithms (see for example Exponentiation by squaring) for integer exponentiation than the linear one you showed, but I'm unsure whether it would be worth considering them unless you deal with arbitrary precision arithmetic, or an embedded system with no floating point unit.

eerorika
  • 232,697
  • 12
  • 197
  • 326
  • If you have a better way to avoid the truncation and to clamp using `pow` I'm all ears. What I've got was the best I could come up with, and it is very expensive. – Jonathan Mee Jun 15 '17 at 15:30
  • Why do you need the generic programming? – Malcolm McLean Jun 15 '17 at 15:34
  • 1
    @JonathanMee as a clue, N^8 = N^4 * N^4 so you don't need to calculate N^4 twice. – Klitos Kyriacou Jun 15 '17 at 15:35
  • @JonathanMee how about the one I added to the answer? – eerorika Jun 15 '17 at 15:42
  • @KlitosKyriacou Hmmm... I see what you're saying, I could have a logrithmic reduction in my loop count if I could just figure out how to do this... – Jonathan Mee Jun 15 '17 at 15:45
  • @user2079303 Probably is as good as we can do :/ We could inline too I suppose: `static_cast(min(numeric_limits::max(), round(pow(base, exp))))` – Jonathan Mee Jun 15 '17 at 15:48
  • 4
    Problem with wrapping `std::pow` is that it works with `double`, and (at least on some platforms) a `long long` has some values that cannot be precisely represented by a `double`. – Klitos Kyriacou Jun 15 '17 at 16:03
  • 1
    @KlitosKyriacou good point. If `long long` is needed, I would write a specialization with explicit conversions to `long double` to use that overload of `std::pow`. – eerorika Jun 15 '17 at 16:08
  • @klitos-kyriacou That's a very good point. This could make a very interesting example as a recursive function. – Michaël Roy Jun 15 '17 at 16:13
  • @user2079303 `long double` isn't necessarily able to store `long long`, because on platforms other than x86 it's still mostly the same as `double` – phuclv Jun 16 '17 at 02:51
  • It's actually certain that it can't. not enough mantissa bits available. – Michaël Roy Jun 16 '17 at 02:53
  • @LưuVĩnhPhúc well, you could check numeric_limits whether there are enough significand digits to represent the target values. If you really need to support long longs on all platforms, then perhaps fall-back to purely integer based algorithm. – eerorika Jun 16 '17 at 08:43
  • @MichaëlRoy not at all certain. 80-bit long doubles are very typical, and have just the appropriate number of mantissa bits to represent signed 64 bit integers. – eerorika Jun 16 '17 at 08:44
  • 1
    MS still doesn't support anything longer than the 64 bit double. – Michaël Roy Jun 16 '17 at 10:04
  • @MichaëlRoy or rather: no longer supports :) They used to. – eerorika Jun 16 '17 at 10:28
2

Your code did not compile, you should if possible check your code compiles first, use your compiler, or check it on compiler explorer first.

Also, you forgot to take into account negative values. That's a very important feature of integral powers. Code below is for regular int type. I'll let you explore how you could extend it for other integral types.

#include <type_traits>
#include <iostream>
#include <cmath>
#include <limits>

using namespace std;

template <typename T>
enable_if_t< is_integral<T>::value, T> 
mypow(T base, unsigned int exp) 
{
    T result = T(1);
    bool sign = (base < 0);

    if (sign) base = -base;

    T temp = result;
    while(exp-- != 0) 
    {
        temp *= base;
        if (temp < result)
        {
            return (sign) ? numeric_limits<T>::min() 
                          : numeric_limits<T>::max();
        }
        result = temp;
    }
    return (sign && (exp & 1)) ? -result : result;
}

template <typename T>
enable_if_t< !is_integral<T>::value, int> 
mypow(const T& base, unsigned int exp) 
{
    T result = T(1);
    int i_base = int(floor(base + .5));
    bool sign = (i_base < 0);

    if (sign) i_base = -i_base;

    int temp = result;
    while(exp-- != 0) 
    {
        temp *= i_base;
        if (temp < result)
        {
            return (sign) ? numeric_limits<int>::min() : numeric_limits<int>::max();
        }
        result = temp;
    }
    return (sign && (exp & 1)) ? -result : result;
}

In real life, I would do this note the use of floor, even in the integral case.

  template<typename T>
   enable_if_t< is_integral<T>::value, T> 
   mypow(T x, unsigned int y) { return T(floor(pow(x, y) + .5)); }

   template<typename T>
   enable_if_t< !is_integral<T>::value, int> 
   mypow(T x, unsigned int y) { return int(floor(pow(floor(x + .5), y) + .5)); }
Michaël Roy
  • 6,338
  • 1
  • 15
  • 19
  • I only intended this as an example, but yeah, this is one more reason that I don't want to go with writing my own function :( – Jonathan Mee Jun 15 '17 at 15:54
  • Using floating point pow will generally be faster on CPU that support floating point ops. Some modern CPUs have a floating-point pow instruction. – Michaël Roy Jun 15 '17 at 16:03
  • 1
    @JonathanMee: Maybe you should check this interesting discussion on the subject: https://stackoverflow.com/questions/101439/the-most-efficient-way-to-implement-an-integer-based-power-function-powint-int – Michaël Roy Jun 15 '17 at 16:25
  • If `x` `is_integral`, then why are you doing `floor(x+.5)` rather than `double(x)`? are you worried about 64bit integers not representable as `double`? They screw your code anyway. – Walter Jun 15 '17 at 16:44
  • @Walter sorry, it's a typo in second function, which is for floating points. I'm rounding the value, since that was part of the original question. Thanks for pointing this out.. – Michaël Roy Jun 16 '17 at 02:41