Editing Exponentiation by squaring (section)

==Basic method==

===Recursive version===
The method is based on the observation that, for any integer <math>n > 0</math>, one has:
<math display="block"> x^n=
    \begin{cases}
                x \, ( x^{2})^{(n - 1)/2}, & \mbox{if } n \mbox{ is odd} \\
                (x^{2})^{n/2} , & \mbox{if } n \mbox{ is even}
     \end{cases}
</math>

If the exponent {{mvar|n}} is zero then the answer is 1. If the exponent is negative then we can reuse the previous formula by rewriting the value using a positive exponent. That is,
<math display="block">x^n = \left(\frac{1}{x}\right)^{-n}\,.</math>

Together, these may be implemented directly as the following [[recursion (computer science)|recursive algorithm]]:

  In: a real number x; an integer n
  Out: x<sup>n</sup>
  
  exp_by_squaring(x, n)
    if n < 0 then
       return exp_by_squaring(1 / x, -n);
    else if n = 0 then 
       return 1;
    else if n is even then 
       return exp_by_squaring(x * x, n / 2);
    else if n is odd then 
       return x * exp_by_squaring(x * x, (n - 1) / 2);
  end function 

In each recursive call, the least significant digits of the [[binary representation]] of {{mvar|n}} is removed. It follows that the number of recursive calls is <math>\lceil \log_2 n\rceil,</math> the number of [[bit]]s of the binary representation of {{mvar|n}}. So this algorithm computes this number of squares and a lower number of multiplication, which is equal to the number of {{math|1}} in the binary representation of {{mvar|n}}. This logarithmic number of operations is to be compared with the trivial algorithm which requires {{math|''n'' − 1}} multiplications.

This algorithm is not [[Tail call|tail-recursive]]. This implies that it requires an amount of auxiliary memory that is roughly proportional to the number of recursive calls -- or perhaps higher if the amount of data per iteration is increasing.

The algorithms of the next section use a different approach, and the resulting algorithms needs the same number of operations, but use an auxiliary memory that is roughly the same as the memory required to store the result.

===With constant auxiliary memory===

The variants described in this section are based on the formula
:<math> yx^n=
    \begin{cases}
                (yx) \, ( x^{2})^{(n - 1)/2}, & \mbox{if } n \mbox{ is odd} \\
                y\,(x^{2})^{n/2} , & \mbox{if } n \mbox{ is even}.
     \end{cases}
</math>

If one applies recursively this formula, by starting with {{math|1=''y'' = 1}}, one gets eventually an exponent equal to {{math|0}}, and the desired result is then the left factor.

This may be implemented as a tail-recursive function:
<syntaxhighlight lang="pascal">
  Function exp_by_squaring(x, n)
    return exp_by_squaring2(1, x, n)
  Function exp_by_squaring2(y, x, n)
    if n < 0 then return exp_by_squaring2(y, 1 / x, -n);
    else if n = 0 then return y;
    else if n is even then return exp_by_squaring2(y, x * x, n / 2);
    else if n is odd then return exp_by_squaring2(x * y, x * x, (n - 1) / 2).
</syntaxhighlight>

The iterative version of the algorithm also uses a bounded auxiliary space, and is given by

<syntaxhighlight lang="pascal">
  Function exp_by_squaring_iterative(x, n)
    if n < 0 then
      x := 1 / x;
      n := -n;
    if n = 0 then return 1
    y := 1;
    while n > 1 do
      if n is odd then
        y := x * y;
        n := n - 1;
      x := x * x;
      n := n / 2;
    return x * y
</syntaxhighlight>

The correctness of the algorithm results from the fact that <math>yx^n</math> is invariant during the computation; it is <math>1\cdot x^n=x^n</math> at the beginning; and it is <math>yx^1=xy </math> at the end.

These algorithms use exactly the same number of operations as the algorithm of the preceding section, but the multiplications are done in a different order.