What is the correct way to raise an integer to a positive integer power in C++?

Question

We know that for various reasons, there is no standard integer power function in C++. I'm performing exact arithmetic with rather small integers, what is the correct way to compute powers?

What is wrong with [std::pow](http://en.cppreference.com/w/cpp/numeric/math/pow)? — juanchopanza, Sep 24 '12 at 09:42
How about `pow(x, y) = exp(y * log(x))`? However you get three floating point errors here with weird numbers. Or could you explain what you're up to? — Christian Ivicevic, Sep 24 '12 at 09:42
`std::pow` uses floating point numbers, leading to potential innacuracies. I want the exact value. — static_rtti, Sep 24 '12 at 09:45
What's wrong with repeated multiplication since you're only interested in small integer powers ? — High Performance Mark, Sep 24 '12 at 09:45
Repeated multiplication with a tail-recursive function could be what you want. — Christian Ivicevic, Sep 24 '12 at 09:46
@HighPerformanceMark : if that's indeed the best / fastest solution, nothing. — static_rtti, Sep 24 '12 at 09:46
You might want to look at the individual bits of the exponent... — ltjax, Sep 24 '12 at 09:48
This is a non-dupe in the sense that the question over there was specifically labelled "homework", and it got home-worky answers. This question has a highly-upvoted answer with an iterative solution by-squaring. Whatever we think about the likelihood of a C++ compiler making a tail-call optimization, it's how a *lot* of real programmers would write it in practice. That other question didn't draw any such answer at all. — Steve Jessop, Sep 24 '12 at 11:09
Have to agree with Steve here that this isn't an "exact dupe", so I've re-opened. — Kev, Sep 24 '12 at 23:23

Kerrek SB · Accepted Answer · 2012-09-26T11:34:58.113

45

The standard, fast exponentiation uses repeated squaring:

uint_t power(uint_t base, uint_t exponent)
{
    uint_t result = 1;

    for (uint_t term = base; exponent != 0; term = term * term)
    {
        if (exponent % 2 != 0) { result *= term; }
        exponent /= 2;
    }

    return result;
}

The number of steps is logarithmic in the value of exponent. This algorithm can trivially be extended to modular exponentiation.

Update: Here is a modified version of the algorithm that performs one less multiplication and handles a few trivial cases more efficiently. Moreover, if you know that the exponent is never zero and the base never zero or one, you could even remove the initial checks:

uint_t power_modified(uint_t base, uint_t exponent)
{
    if (exponent == 0) { return 1;    }
    if (base < 2)      { return base; }

    uint_t result = 1;

    for (uint_t term = base; ; term = term * term)
    { 
        if (exponent % 2 != 0) { result *= term; }
        exponent /= 2;
        if (exponent == 0)     { break; }
    }

    return result;
}

edited Sep 26 '12 at 11:34

answered Sep 24 '12 at 09:49

Kerrek SB

464,522
92
875
1,084

2

@PeterWood: It's commonly known, it's general (the same works for multiplication-in-terms-of-addition), it's fast (logarithmic), it's simple, and it exploits some nice machine properties (integer division by two is very fast)... so I'd say it's hard to come up with something better. – Kerrek SB Sep 24 '12 at 09:58
The most voted answer (though, unfortunately, not the accepted one) to http://stackoverflow.com/questions/1505675/power-of-an-integer-in-c proposes a variation on this algorithm, which we can consider more *standard* because it appears in Knuth. But the idea is the same. – Gorpik Sep 24 '12 at 10:27
3

@Gorpik: That answer doesn't look tail-recursive to me. I wouldn't feel comfortable suggesting that for anything but an educational or theoretical work. – Kerrek SB Sep 24 '12 at 10:29
Indeed, it is not tail-recursive, just recursive. And I agree that your solution would work better; in fact, I have upvoted it. – Gorpik Sep 24 '12 at 10:42

score 19 · Answer 2 · edited Sep 25 '12 at 23:42

19

You can use std::pow(double a, double b). There will be no inaccuracies if both a, b and the result fit into a 32-bit integer!

The reason is that the 64-bit double precision fully covers the range of 32-bit integers.

edited Sep 25 '12 at 23:42

eboix

5,113
1
27
38

answered Sep 25 '12 at 22:24

PlinyTheElder

1,454
1
10
15

You still have to round properly here, since the double result might be something like 15.9999999 which could be truncated if carelessly converted back to an int. – NovaDenizen Sep 25 '12 at 23:16
15

@NovaDenizen No, the doubles cover the 32-bit integers exactly. In other words, there is a double with exactly zero fractional part for every integer. And pow(), like other arithmetic operations, must give "the closest possible" double answer, which in this case will equal to an integer. – PlinyTheElder Sep 25 '12 at 23:44
1

Not completely accurate for all 64-bit integers, but interesting answer nevertheless. – static_rtti Sep 26 '12 at 07:32
2

@static_rtti: It's true for all integers that fit entirely into the mantissa, which is 24 bits for a `float`, 53 bits for a `double` and 64 bits for a `long double` in IEEE754. – Kerrek SB Sep 26 '12 at 13:30
@KerrekSB: so which one should be used, in which circumstances ? The two answers are very different methods... – static_rtti Sep 26 '12 at 13:40
@static_rtti: Measure and compare! I think the pure-integer one should be better, but I haven't actually measured it. That said, the pure-integer one can use 64-bit integers all the way... – Kerrek SB Sep 26 '12 at 13:44
2

@user1698678 The C standard requires that all calculations be accurate to within 1 ULP (unit of least precision), which means the true answer can be +/- 1 ULP, so my caution still stands. If C used a 1/2 ULP accuracy standard then you would be correct. – NovaDenizen Sep 26 '12 at 13:49
2

@static_rtti: Quick test, 10M times 7^11: Integer (fast 32 bit): 0.5s. Double: 9.0s. – Kerrek SB Sep 26 '12 at 13:58
@KerrekSB : thanks for testing! Your test was on x86-64, I suppose? – static_rtti Sep 26 '12 at 14:41
If you're worried about rounding, why not overload pow for ints and wrap the return of the double version in an std::lround? http://en.cppreference.com/w/cpp/numeric/math/round – SplinterOfChaos Sep 26 '12 at 17:50

score 4 · Answer 3 · answered Sep 26 '12 at 22:56

4

While Kerrek's answer is correct, there is also a "secret" feature in g++ to do this efficiently. If you look at the SGI power function, it can be easily adapted to do what you want:

http://www.sgi.com/tech/stl/power.html

In g++, this is implemented as __gnu_cxx::power. You probably shouldn't use these things in production code though...

answered Sep 26 '12 at 22:56

Mikola

9,176
2
34
41

Woohoo, awesome, thanks a lot for sharing! I couldn't care less about Windows compatibility for my code :-D – static_rtti Sep 27 '12 at 07:42
1

Good to know - I'll add this to my test and compare. People have been asking for a `power` function in the standard library for some time; I think they have that SGI function in mind. – Kerrek SB Sep 27 '12 at 09:47
1

Testing this on x86_64 with a `uint_fast64_t`, `__gnu_cxx::power` is slightly slower than my algorithm (3.5s vs 2.9s on 1bn computations of 12^15). – Kerrek SB Sep 27 '12 at 11:43
@KerrekSB well duh the sgi version uses templates – pyCthon Sep 28 '12 at 15:50
1

@pyCthon: templates are evaluated at compile-time, they have no impact whatsoever on runtime performance. – static_rtti Oct 09 '12 at 08:50

score 1 · Answer 4 · answered Sep 30 '12 at 03:21

1

In addition to the other answers here, there is also a nice article on Wikipedia explaining various different implementations here -> LINK

answered Sep 30 '12 at 03:21

キキジキ

1,443
1
25
44

What is the correct way to raise an integer to a positive integer power in C++?

4 Answers4

Linked