Exponentiation, or calculating a power of a number, looks like a trivial task. Just multiply a number by itself several times. But, is it the most efficient way? In this article we’ll discuss several approaches to what seems to be a basic mathematical operation. We will only discuss powers that are integers here.

First, let’s take a look at the approaches that can be encountered all over the Internet.

The straightforward approach

First, let’s do it in Python.

def powah_naive(x, k_init):
    k=abs(k_init)
    
    #zero in negative power is undefined
    #https://math.stackexchange.com/a/1009872/62757
    if x == 0 and k_init < 0:
        return float('nan')
    
    #####################
    #the actual algorithm
    #####################
    
    result = 1
    for i in range(k):
        result *= x
        
    ###############################
    #the actual algorithm ends here
    ###############################
            
    #catering for negatives
    if k_init < 0:
        result = 1/result
        
    return result

And a test code. It will be pretty much the same for every approach (except the function name will be different).

from random import randint

for i in range(100000):
    x = randint(-1000,1000)
    k = randint(-1000,1000)
    
    if x == 0 and k<0:
        continue

    #zero in negative power is undefined
    #https://math.stackexchange.com/a/1009872/62757
    if powah_naive(x,k) != x**k:
        print(f"Wrong on x={x} k={k}")
        break

I didn’t want to include negative values of k at first, but since I was playing with those anyway, why not. Also, I learnt that zero in negative power is a complicated matter, so I’ll just leave it be. Python throws ZeroDivisionError: 0.0 cannot be raised to a negative power on 0**-1 anyway.

This approach takes O(n) time and constant space. Not bad, but it can be done faster.

The recursive approach

The fundamental quality of powers is that you can multiply two powers and their exponents sum up: $ x^a * x^b = x^{(a+b)}$. You can go the other way around by dividing exponents evenly until you reach the value of 1. It can be easily done in a recursive way.

def _powah_recursive_algorithm(x, k):
    if k == 0:
        return 1
    
    if k == 1:
        return x
    
    p = _powah_recursive_algorithm(x, int(k/2))
    if k%2 == 1:
        #odd power
        return p * p * x
    else:
        #even power
        return p * p

def powah_recursive(x, k_init):
    k=abs(k_init)
    
    #zero in negative power is undefined
    #https://math.stackexchange.com/a/1009872/62757
    if x == 0 and k_init < 0:
        return float('nan')
    
    #####################
    #the actual algorithm
    #####################
    
    result = _powah_recursive_algorithm(x,k)
        
    ###############################
    #the actual algorithm ends here
    ###############################
            
    #catering for negatives
    if k_init < 0:
        result = 1/result
        
    return result

This approach takes O(log k) time, but it requires more space for recursive calls (probably O(log k) by memory, correct me if I’m wrong).

Iterative approach

As far as I’ve heard, any task that can be performed recursively can be done iteratively as well. Can we go deeper and make the algorithm satisfy the following requirements?

  • O(log k) by time
  • O(1) by memory
  • No modulo or division operations (except the last one for negative powers)

It is possible, with some bit magic.

First of all, let’s look at what an integer actually is from a binary representation perspective.

$ k = (0 \text{ or } 1) \cdot 2^{0} + (0 \text{ or } 1) \cdot 2^{1} + (0 \text{ or } 1) \cdot 2^{2} + (0 \text{ or } 1) \cdot 2^{3} + \cdots $

That’s how an integer is represented in memory, in reverse. The bit that corresponds to $2^{0}$ (or $1$) is the rightmost.

Now, back to our exponent. Let’s create a polynomial (or is it a polynomial? Math-savvy people may correct me).

$ x^{k} = x^{(0 \text{ or } 1) \cdot 2^{0} + (0 \text{ or } 1) \cdot 2^{1} + (0 \text{ or } 1) \cdot 2^{2} + (0 \text{ or } 1) \cdot 2^{3} + \cdots} = x^{(0 \text{ or } 1) \cdot 2^{0}} \cdot x^{(0 \text{ or } 1) \cdot 2^{1}} \cdot x^{(0 \text{ or } 1) \cdot 2^{2}} \cdot x^{(0 \text{ or } 1) \cdot 2^{3}} \cdot \cdots $

The rightmost bit corresponds to $ x^{(0 \text{ or } 1) \cdot 2^{0}} $. Which means, $ x $ can be either $ 1 $ if the bit is 0 or $ x $ if it is 1. This can be done in code very simply:

    #if the power is odd, we start with x^1 in the result, because the rightmost bit is 1.
    result = 1
    if k&1: # another way to determine whether a positive number is odd or even, by the way, along with k%2
        result = x

What about the rest? We can easily iterate ofer the remaining bits using k >>= 1 (or k = k >> 1) and getting the next bit with k&1 on every iteration. The next number in this polynomial will be either $ 1 $ or $ x^{2} $. The one after it is either $ 1 $ or $ x^{2^{2}} = x^{4} $. The following one is either $ 1 $ or $ x^{2^{3}} = x^{8} $. See the pattern? We just have to make a number (which I call cur_pow_x in my program) and intialize it to $ x^{2} $, and multiply it by itself on every iteration. If the current bit is 1, then we just multiply the current result by cur_pow_x.

def powah(x, k_init):
    k=abs(k_init)
    
    #zero in negative power is undefined
    #https://math.stackexchange.com/a/1009872/62757
    if x == 0 and k_init < 0:
        return float('nan')
    
    #####################
    #the actual algorithm
    #####################
    
    #if the power is odd, we start with x^1 in the result, because the rightmost bit is 1.
    result = 1
    if k&1: # another way to determine whether a positive number is odd or even, by the way, along with k%2
        result = x
        
    #we have processed the rightmost bit, let's get rid of ot
    k >>= 1
    
    #cur_pow_x will be the power of x corresponding to the current bit
    cur_pow_x = x * x
    
    #iterating over bits until the exponent is zero
    while k:
        if k&1:
        	#if this bit is 1, it influences the result
            result *= cur_pow_x
        
        #moving to the next bit and the current member of the polynomial
        k = k >> 1
        cur_pow_x *= cur_pow_x
    
    ###############################
    #the actual algorithm ends here
    ###############################
            
    #catering for negatives
    if k_init < 0:
        result = 1/result
        
    return result

There you go. Iterative approach.

Benchmarks

I decided to run some benchmarks for Python algorithms on every approach we discussed so far for the $ 9^{9} $ operation (this way I stayed within the boundaries of a 32-bit integer. Python may get slow on higher numbers because it has to increase the number of operations on numbers that are beyond limitations of the CPU architecture). The results were slightly unexpected:

powah_naive: 1.29 µs ± 40.6 ns
powah_recursive: 2 µs ± 25.4 ns
powah_iterative: 1.39 µs ± 31.3 ns

The O(k) algorithm was the fastest. But we have to remember that 9 is likely a pretty small number, and some time is consumed by the overhead of additional operations in the two latter algorithms. So they might be slower on small numbers.

Let’s do the same thing for $ 9^{99} $. The results:

powah_naive: 10.1 µs ± 158 ns
powah_recursive: 3.99 µs ± 160 ns
powah_iterative: 2.68 µs ± 94.8 ns

Now we can finally see the might of O(log k). And the iterative approach is faster. But what if we go for $ 9^{999} $ now?

powah_naive: 177 µs ± 2.95 µs
powah_recursive: 10.3 µs ± 125 ns
powah_iterative: 14.1 µs ± 218 ns

Now the recursive one is faster. And I have tested it with higher values of exponent, the iterative algorithm only gets slower compared to the recursive one. I don’t have an explanation for this at the moment.

In Golang

As a bonus, I’ll include the same code in Golang, since I played around with it as well. Didn’t make any benchmarks though.

package main

import (
	"fmt"
	"math"
	"math/rand"
	"time"
)

func main() {
	s1 := rand.NewSource(time.Now().UnixNano())
	r1 := rand.New(s1)

	for i := 0; i <= 10000000; i++ {
		x := float64(r1.Int63n(100))
		n := r1.Int63n(10)

		//fmt.Println(x, n, powah(x, n), math.Pow(x, float64(n)))

		if powah_naive(x, n) != math.Pow(x, float64(n)) {
			fmt.Println("Nope on naive!", x, n, powah_naive(x, n), math.Pow(x, float64(n)))
			break
		}

		if powah_iterative(x, n) != math.Pow(x, float64(n)) {
			fmt.Println("Nope on iterative!", x, n, powah_iterative(x, n), math.Pow(x, float64(n)))
			break
		}

		if powah_recursive(x, n) != math.Pow(x, float64(n)) {
			fmt.Println("Nope on recursive!", x, n, powah_recursive(x, n), math.Pow(x, float64(n)))
			break
		}
	}
}

func powah_naive(x float64, n int64) float64 {
	var result float64 = 1

	for i := int64(0);i < n; i++{
		result *= x
	}

	return result
}

func powah_iterative(x float64, n int64) float64 {

	var result float64 = 1
	if n&1 == 1 {
		result = x
	}
	n = n >> 1

	var cur_pow_x float64 = x * x

	for n > 0 {
		if n&1 == 1 {
			result *= cur_pow_x
		}
		n = n >> 1

		cur_pow_x = cur_pow_x * cur_pow_x
	}
	return result
}

func powah_recursive(x float64, n int64) float64 {
	//fmt.Println("n", n)

	if n == 0 {
		return 1
	}

	if n == 1 {
		return x
	}

	p := powah_recursive(x, n/2)
	if n%2 == 1 {
		//odd power
		return p * p * x
	}

	//even power
	return p * p
}

Keep in ming that the precision may be lost on big numbers. If you set the upper limit for n to, say, 100, the test code will throw errors, even though the results are seemingly close. Strangely enough, powah_iterative did not have this behaviour, it was the most precise.

Conclusion

There you have it! Several approaches to exponentiation. Sounds like a simple task but I’ve spent hours playing around with it. Of course, the practical approach is by simply using x**k in Python or math.Pow(x,k) in Go, the built-ins are way faster and more available.

Thanks for tuning in!