Fast Exponentiation, revisited

Quite a while ago, I presented a fast exponentiation algorithm that uses the binary decomposition of the exponent n to perform O(\log_2 n) products to compute x^n.

While discussing this algorithm in class, a student asked a very interesting question: what’s special about base 2? Couldn’t we use another base? Well, yes, yes we can.

Indeed, there’s nothing special about base 2, except maybe that it’s very simple, and it yields a very simple implementation. For example:

unsigned expo_iter(unsigned x, unsigned p)
 {
  unsigned t=x, e=1;

  while (p)
   {
    if (p&1) e*=t;

    t*=t;
    p/=2;
   }

  return e;
 }

This simple procedures raises x to the pth power in O(\log_2 p) steps. And they’re (mostly) inexpensive: a shift right by one bit, a mask, and a multiply (which could be expensive if the objects are complex and we can’t rely on machine-sized multiplications).

But what about other bases? Say, base 7? or 3? or whatever other than 1?

First, we must understand what the binary algorithm does. It exploits the fact that x^a\times{}x^b=x^{a+b} in a very specific way: the exponent will be broken into powers of two, so that x^{25} becomes x^{16}x^{8}x; and indeed the powers to use are given by the binary representation of 25, 11001 (which is 2^{16}+2^{8}+1!). Now, let us see that the above algorithm does. The following table gives you what the algorithm computes at step i:

Step e p t
(init) 25 1 x
1 25 x x^2
2 12 x x^4
3 6 x x^8
4 3 x^8x x^16
5 1 x^{16}x^8x x^32
(done) 0 x^{16}x^8x

So it is clear that the number of iterations is given by the position of the most significant bit in the exponent (assuming it’s integer!), and we do at most two multiply (and one shift) per round.

*
* *

So how do we generalize this algorithm from base 2 to some base b>1? Certainly, the decomposition of the exponent in base b is necessary. Let’s pick b=3 (because it’s not 2 and it’s also not too big) and e=98=10122_3 (also, no special reason other than a good mix of digits and just long enough we can figure something out of them). We will use the digits to multiply each group of powers of 3 so that the total is 98. Indeed, we have 1\cdot{}3^4+0\cdot{}3^3+1\cdot{}3^2+2\cdot{}3^1+2\cdot{}3^0 =81+0+9+6+2=98. Indeed:

0

Step e m p t
(init) 10122_3 1 x
1 1012_3 2 x^2 x^3
2 101_3 2 x^2(x^3)^2=x^2x^6 x^9
3 10_3 1 x^2x^6x^9 x^27
4 1_3 0 x^2x^6x^9 x^{81}
5 0 1 x^2x^6x^9x^{81} x^{243}
(done) 0 x^2x^6x^9x^{81}

Now, let’s turn that into code! Here, Mathematica will help us deal with arbitrary large numbers.

expo[x_, e_, b_: 2] := Module[{tx = x, p = 1, te = e, m},
  While[te != 0,
   m = Mod[te, b];
   If[m != 0, p *= tx^m];

   te = Quotient[te, b];
   If[te != 0, tx = tx^b]; (* if product is expensive *)
   ];
  p
  ]

*
* *

This algorithm trades-off the number of iterations (from \approx\log_2 p to \approx\log_b p to raise to the pth power) for a more complex inner loop: div/mod by something else than 2 could be expensive—even if divmod is likely just one instruction!

2 Responses to Fast Exponentiation, revisited

  1. Kristine says:

    Your “a while ago” link at the start seems to be empty?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: