Exponentiation involving doubles is easily handled using a formula from the logs page.

The case of an integer exponent would seem to be trivial: to calculate aⁿ, just multiply a by itself n, times. We are interested in the number of multiplications required, in this case it is n - 1, so this is a O(n) algorithm. (Notice that, while 2 * 2 is easier for you to do than 43046721 * 81, on a computer all multiplications take the same amount of time.)

There are much better exponentiation algorithms. In fact, these better algorithms are an absolute requirement for the RSA cryptosystem, used for security all over the planet.

Let's look at calculating 3²³ = 94 143 178 827. If we want to jack ourselves up to higher powers of 3, we could try repeated squaring, at least until we get bigger than 3²³:

So where are we? The last usable number is 3¹⁶ = 43 046 721. To get up to 3²³, we need 7 more factors of 3. Well, the table above has 3⁴ = 81, so let's multiply that in:

3¹⁶ * 3⁴ = 3¹⁶⁺⁴ = 3²⁰ = 43 046 721 * 81 = 3 486 784 401

3¹⁶ * 3⁴ * 3² * 3¹ = 3^16+4+3+1 = 3²³ = 43 046 721 * 81 * 9 * 3 = 94 143 178 827.

If you look above, the table shows 3 raised to powers which themselves are successive powers of 2. We were able to get to 3²³ by representing 23 as a sum of powers of two:

23 = 16 + 4 + 2 + 1. But this is just the binary representation of 23 = 1 * 2⁴ + 0 * 2³ + 1 * 2² + 1 * 2¹ + 1 * 2⁰, so that 23 = 10111₂. The 1s in the binary representation of the exponent tell you which table entries to multiply together.

Suppose we have a base a and an exponent b. We want a^b. At each stage we need to raise a to successive powers: 2, 4, 8, 16, 32, 64, .... We will use the parameter a below to hold these powers. At each stage we also need to generate successive bits of the exponent, from the least significant to the most significant. For each 1 bit, we will multiply in the power created at that stage, using, say, a variable d to hold these successive results (d starts out at 1). d is the value returned.

A clever way to produce successive binary bits from right to left is to take b%2 for the low order bit, and then replace b with b / 2 to drop the low order bit; then repeat.

Here's the C code I ended up with. For comparison, I stuck in the slow exponentiation, with a function exps.

// exponen.c: fast exponentiation
#include <stdio.h>
int exponentiation(int, int);

int main() {
   int d;
   int a, b;
   scanf("%i %i", &a, &b);
   d = expon(a, b);
   printf("a: %i, b: %i, d: %i\n",
      a, b, d);
}

// expon: fast exponentiation
int expon(int a, int b) {
   int d = 1;
   while(b > 0) {
      if (b%2 == 1) d = d*a;
      b = b/2;
      if (b > 0) a = a*a;
   }
   return d;
}

// exps: slow exponentiation
int exps(int a, int b) {
   int d = 1;
   while(b > 0) {
      d = d*a;
      b = b--;
   }
   return d;
}
Several Runs:
% cc -o exponen exponen.c
./exponen
3 19
a: 3, b: 19, d: 1162261467
a: 5, b: 13, d: 1220703125
a: 5, b:  2, d: 25
a: 5, b:  4, d: 625
a: 5, b:  5, d: 3125
a: 3, b: 13, d: 1594323

% ruby power.rb 2 127
170141183460469231731687303715884105728
(2^127, 8 loop iters, log2(127) = 6.98868)

% ruby power.rb 3 236
398672347901056050310521584754736038852
144821237913809261660811330548546381908
73316632624557307582586552837323121
(3^236, 9 loop iters, log2(236) = 7.8826)

% ruby power.rb 2 431
55453393882416297191568283682861674068728741
50751633150340959161229242615611251246079948
812208279156194782421922807143657948315648
(2^431, 10 loop iters, log2(431) = 8.7515)

% ruby power.rb 2 814
10924874846830353242652429173917708774841972
03528238392999322879879419861044040838598071
45179079545499532163655162810211465442126301
09167747667322556138831947522306178275761452
58590746363791609811593628243703660868680212
15846014341794428966928384
(2^814, 11 loop iters, log2(814) = 9.668885)

% ruby power.rb 2 2011
23513716639216732656288423960118927032777066
53877647770578212325019169303331574402778438
62797620484079297733558064632376001969963184
75373187603074174293057236813821776959954944
14106196955956118966864424756500198293127484
69547803369364105269330875092932413432951116
59276071860143231654225825006432659524133139
28579256001710194736731848866079926291144222
72878597563815169700909072654819904085527076
62503363911474257186432414060515001147834337
93192991557415606314842532965535407916483361
99847848379659433483019111595470879802058314
58699684763425848480539161635929196298408518
1994635193249304954575153212162048
(2^2011, 12 loop iters, log2(2011) = 10.9737)

This algorithm has a loop iteration for each bit in the binary form of the exponent. If n is the exponent, then log₂ n gives the number of bits in its binary representation (see logarithms). We do one or two multiplications during each loop iteration except the last, which has zero or one. So the number of multiplications is between log₂ n and 2 * log₂ n. Thus this is an O(log n) algorithm.

As n gets large, the difference between it and log₂ n becomes phenomenal. For example, if the exponent is 1 000 000, the simple algorithm would require that many mulitiplications, while the better algorithm here would need only log₂ 10⁶ = 6 * log₂ 10 = 6 * 3.321928095 = 19.93156857 loop interations (that is, 20), and up to twice that many multiplications. Here the comparison is between forty and a million multiplications. But who would use such large exponents? The RSA cryptosystem routinely uses exponents that are 200 decimal digits (about 666 bits) long. So exponentiation needs only 666 loop iterations, rather than 10²⁰⁰ iterations, a number that will forever be impossible for one of our computers to carry out.

From the way this program has been developed, it should be plausible that it is correct. However, a formal correctness proof is important, when we can give it. Such a proof for a simple program sometimes has a loop invariant that should remain true at the beginning and end of each loop. For this proof, it also helps to introduce separate variables for the formal parameters. The execution is the same as before.

// expon: fast exponentiation
int expon(int A, int B) {
   int a = A, b = B, d = 1;
   while(b > 0) {
      if (b%2 == 1) d = d*a;
      b = b/2;
      if (b > 0) a = a*a;
   }
   return d;
}

The loop invariant is:

d * a^b = A^B

Case 1: At the beginning of the loop, the rightmost bit of b is a 1.

Then d' = d * a, a' = a * a, and b' = (b-1)/2 (since b is odd). Substituting, d' * a'^b' = d * a * (a * a)^(b-1)/2 = d * a * (a²)^(b-1)/2 = d * a * (a)^2*((b-1)/2) = d * a * (a)^(b-1) = d * a^b = A^B.

Case 2: At the beginning of the loop, rightmost bit of b is 0, but b > 0.

Then d' = d, a' = a * a, and b' = b/2, since b is even. Substituting, d' * a'^b' = d * (a * a)^b/2 = d * (a²)^b/2 = d * (a)^2*(b/2) = d * a^b = A^B.

Coming into the last loop iteration, b will have value 1. At the beginning of this last iteration loop, d * a^b = d * a¹ = d * a = A^B. At the end of the loop, we have d' = d * a, a' = a * a, and b' = 0, so that d' * a'^b' = d * a * (a * a)⁰ = d * a * (a²)⁰ = d * a * 1 = d * a = A^B. So at the end of this final loop, we have d' = A^B, and that's what is returned. (In the final loop iteration, the computation a = a * a is not needed, so I have cut it out.)

Finally, we have to argue that the loop terminates. At each iteration, the loop drops one bit of the variable b, so assuming b > 0 initially, it will terminate. Notice that if b is zero initially, the program will terminate with d = 1 and return 1 as it ought to.

Your text has a similar-looking algorithm in their presentation of the RSA cryptosystem toward the end of the book, but it has significant differences from the one here. [See 3rd Ed., page 957; 2nd Ed., page 879.] With RSA, one is constantly taking the remainder on dividion by a large number, so the size of the numbers never gets out of hand.

Ruby supports arbitrarily large integers, as do many of the scripting languages.

// Exp: Java version of exponentiation
// Use long type to get 64-bit integers
public class Exp {

   public static long exp2(long X, long Y) {
      long x = X, y = Y, z = 1;
      while (y > 0) {
         if (y%2 == 1) z = z*x;
         y = y/2;
         if (y > 0) x = x*x;
      }
      return z;
   }

   public static void main(String[] args) {
      long x = Long.parseLong(args[0]);
      long y = Long.parseLong(args[1]);
      long z = Exp.exp2(x, y);
      System.out.println("exp(" + x + ", " +
         y + ") = " + z);
   }
}
% java Exp  3  19
exp(3, 19) = 1162261467

% java Exp  3  23
exp(3, 23) = 94143178827

% java Exp  3  38
exp(3, 38) = 1350851717672992089

% ruby exp.rb 3 38
1350851717672992089

def expn(a, b)
   d = 1
   iter = 1
   while b > 0
      iter = iter + 1
      if b%2 == 1 # right 1 bit
         d = d*a
      end
      b = b/2 # drop right bit
      if b > 0
         a = a*a
      end
   end
   print "Iterations: ", iter, "\n"
   d
end
   
x = ARGV[0].to_i
y = ARGV[1].to_i
print "exp(", x, ",", y, "):"
print expn(x, y), "\n"
% ruby exp.rb 3 45
Iterations: 7
exp(3,45):2954312706550833698643

% ruby exp.rb 17 17
Iterations: 6
exp(17,17):827240261886336764177

% ruby exp.rb 47 1
Iterations: 2
exp(47,1):47

% ruby exp.rb 47 0
Iterations: 1
exp(47,0):1

%ruby exp.rb 29 29
Iterations: 6
exp(29,29):256768615316121113456
1828214731016126483469