The Power Algorithm

March 19 2014

In this blog post I would like to show how a very basic idea like raising a number to a certain power could lead us to discover abstractions like Semigroups and Monoids.

There’s a very well known algorithm for calculation powers, that is x to the power of n or simply: x^n. Donald Knuth presents the algorithm in section 4.6.3 Evaluation of Powers of TAOCP.

The naïve way to implement this algorithm would be to multiply x by itself n times, but of course the idea is to provide an algorithm that’s faster than this. The algorithm in question is usually called the binary method, the powering ladder or the repeated-squaring algorithm.

Say we want to calculate 2^23, where x = 2 and n = 23 this algorithm takes into account the binary representation of 23 (10111), scans it and depending if it finds a 0 or a 1 it will either square x or multiply by x.

The problem with this method is that it scans the binary representation of the number from left to right, but for computers it’s usually easier to do it the other way around, so Knuth proposes a alternative algorithm.

A simple implementation of Algorithm A from section 4.6.3 of TAOCP is as follow:

<?php
function power1($x, $n) {
    $y = 1;

    while (true) {
        $t = $n % 2;
        $n = floor($n/2);

        if ($t == 1) {
            $y = $y * $x;
        }

        if ($n == 0) {
            break;
        }

        $x = $x * $x;
    }

    return $y;
}

This function takes two integers, $x and $n and returns $x to the power of $n.

First it sets an auxiliary variable $y to 1 which is the identity of the multiplication.

Then the function scans the binary representation of $n on each iteration of the loop. If it finds a 1 then it multiplies by $x. On every step of the loop it squares $x.

Finding a 1 means that the current value of $n is not divisible by 2, or in other words, $n % 2 == 1.

Also on every iteration of the loop $n gets halved and then we apply floor to the result. When $n equals 0, we terminate the loop and return the value of $y.

The function power can be called like this:

1024 == power1(2, 10);
=> true

At this point I imagine you being like the guy on this gif:

CSB

While the algorithm really looks like a CSB, it is actually quite effective when it comes to computing with very big numbers. For example there are many primality testing algorithms that depend on variations of this algorithm.

Adding some abstraction

So far nothing fancy, but if we notice that raising a number to a power is actually the same as multiplying the number by itself that many times, we could also see that multiplication is the same as adding a number by itself that many times. So for example 2 * 5 could be computed by doing 2 + 2 + 2 + 2 + 2.

Can we convert our algorithm into a more generic form that would work for multiplication and addition? Of course we can. We just need to change a couple of things.

On the current implementation we set $y to the identity element of multiplication, namely 1. If we want to use the algorithm for addition, then we need to set $y to 0. So we could just pass the value of our identity element to the function.

The second step would be to provide a function to our algorithm that would do either multiplication or addition. To do that we will pass a function that will act as a binary operation, i.e.: a function that takes two arguments. This function needs to follow some rules as well, that is, it has to be associative: a⁡• (b • c) = (a • b) •. Also the type of the result value must be the same type as both input parameters.

Luckily for us addition and multiplication are associative operations, so we can just wrap them in a function and pass it to our power algorithm.

Here’s the new implementation of the algorithm:

<?php
function power2($x, $n, $id, $f) {
    $y = $id;

    while (true) {
        $t = $n % 2;
        $n = floor($n/2);

        if ($t == 1) {
            $y = $f($y, $x);
        }

        if ($n == 0) {
            break;
        }

        $x = $f($x, $x);
    }

    return $y;
}

And we can call it like this:

1024 == power2(2, 10, 1, function ($a, $b) { return $a * $b; });
=> true

As you can see there we passed 1 as our identity element, and we provided a function that takes 2 arguments and returns the result of their multiplication. Nothing fancy.

But now say we want to calculate 2 * 10 using our algorithm. We could wrap the addition operator inside a function and then pass 0 as the identity element like this:

20 == power2(2, 10, 0, function ($a, $b) { return $a + $b; });
=> true

Keep in mind that the operation that we pass to our algorithm must be associative, for example, subtraction can’t be used here since for example 10 - (5 - 3) = 8 while (10 - 5) - 3 = 2.

Appending more abstraction

Mathematically speaking this algorithm works with any algebraic structure that has an associative operation (like multiplication and addition in the case of integers), that is, it works on Semigroups. To quote a book on Group Theory:

A semigroup set S equipped with an associative binary operation •; that is, x • (y • z) = (x • y) • z for all x, y, z ∈ S.

Also, that set must have an identity element, which makes it a Monoid:

A monoid is a set M equipped with an associative binary operation •; together with an identity element e ∈ M that satisfies e • x = x • e = x for all x ∈ M.

What other structures that we use every day in programming can work under the previous conditions? If you are a web developer you don’t need to dig so much to find strings. Strings, using string append as the binary operation and the empty string as the identity element would give use similar results as above. If we wanted to repeat a string n times, we could create the following function:

<?php
function repeat($s, $n) {
    return power2($s, $n, "", function ($a, $b) {
               return $a . $b;
           });
}

And to test it:

"aaaaaaaaaa" == repeat("a", 10);
=> true

Now think about arrays (or lists in other languages). We could want to build an array with n copies of the same element. Here the empty array is the identity element, and in the case of PHP array_merge would be our binary operation.

<?php
function repeat_el($el, $n) {
    return power2(array($el), $n, array(), function ($a, $b) {
            return array_merge($a, $b);
           });
}

And here are our results:

$arr = repeat_el("a", 10);
10 == count($arr);
=> true

As you can see here, something that seemed so simple as raising a number to a certain power gave us a neat algorithm that can be used for several things, like repeating strings or elements inside an array.

Did you just say Haskell?

Since I’ve mentioned a Haskell book, here’s an implementation of the power algorithm in Haskell, using a recursive algorithm based on the book Prime Numbers: A Computational Perspective.

power :: (Eq a, Integral b) => (a -> a -> a) -> a -> b -> a
power f a n
  | n == 1 = a
  | even n = square a (n `div` 2)
  | otherwise = f a (square a ((n-1) `div` 2))
  where
    square a' n' = f (power f a' n') (power f a' n')

And here are the results of several invocations of the function:

*Main> :load pow.hs
[1 of 1] Compiling Main             ( pow.hs, interpreted )
Ok, modules loaded: Main.
*Main> power (*) 2 10
1024
*Main> power (+) 2 10
20
*Main> power (++) "a" 10
"aaaaaaaaaa"

As you can see there, the function takes a function (a -> a -> a) which could be either * or + for integers, or ++ for lists, for example.

I hope you found this blog post interesting, and that I have whetted your appetite for learning maths in relation to programming, since I think the more we grasp these ideas, the better the abstractions we would be able to use.

Alvaro VidelaRead my Thoughts. Follow my Leads.

The Power Algorithm

Adding some abstraction

Appending more abstraction

Further Reading

Did you just say Haskell?