On raising to a power by repeated squaring

Enuntiatio

Propositio IX. Theorema. Let x be any element of a monoid under multiplication (in particular any number), and let n be a non-negative integer. Then the power x^n may be raised, not by n−1 successive multiplications, but by a number of multiplications not exceeding 2·(⌊log₂ n⌋ + 1).

The construction rests upon a single arithmetical identity, which I take as the law beneath the appearance: a power with an even exponent is the square of the power with half the exponent, and a power with an odd exponent is that same square multiplied once by the base. Formally, for n > 0,

x^n = (x^(n/2))² when n is even, and x^n = x · (x^((n−1)/2))² when n is odd,

with the boundary x^0 = 1. Because each appeal to the identity replaces the exponent by its halving, the exponent is exhausted after ⌊log₂ n⌋ + 1 halvings, and the count of multiplications is therefore logarithmic in n, not linear. Q.E.D. follows in the Demonstratio.

Expressio

def power(x, n):
    """Raise x to the non-negative integer power n by repeated squaring.

    Computes x**n using O(log n) multiplications via the identity
        x**n = (x**(n//2))**2            if n is even,
        x**n = x * (x**((n-1)//2))**2    if n is odd,
    with x**0 = 1. Works for any value x that supports multiplication
    (int, float, Fraction, complex, square matrix, ...).

    >>> power(3, 13)
    1594323
    >>> power(2, 10)
    1024
    >>> power(5, 0)
    1
    >>> [power(2, k) for k in range(6)]
    [1, 2, 4, 8, 16, 32]
    """
    if not isinstance(n, int) or n < 0:
        raise ValueError("exponent n must be a non-negative integer")

    result = 1        # invariant accumulator: holds the product of "odd" contributions
    base = x          # holds x**(2**k) for the current bit position k
    e = n
    while e > 0:
        if e & 1:               # current low bit of e is set
            result *= base      # fold in x**(2**k) for this bit
        base *= base            # advance base to x**(2**(k+1)): one squaring
        e >>= 1                 # drop the consumed bit
    return result


# --- concrete demonstrations of correctness (inline assertions) ---
assert power(3, 13) == 1594323          # 3**13, the stated witness
assert power(3, 13) == 3 ** 13
assert power(2, 0) == 1                 # boundary n = 0
assert power(7, 1) == 7                 # boundary n = 1
assert power(2, 10) == 1024
assert power(0, 0) == 1                 # 0**0 taken as the empty product = 1
# Exhaustive cross-check against Python's own exponentiation:
for _x in range(-4, 5):
    for _n in range(0, 30):
        assert power(_x, _n) == _x ** _n, (_x, _n)

if __name__ == "__main__":
    import doctest
    doctest.testmod(verbose=False)
    print("power: all assertions and doctests passed")

Demonstratio

I prove two things: that the procedure returns x^n (correctness), and that it does so within 2·(⌊log₂ n⌋ + 1) multiplications (the bound).

1. The governing identity. For n > 0 write n = 2q + r with r ∈ {0, 1}. By associativity and commutativity of multiplication, x^n = x^(2q + r) = (x^q)² · x^r. When r = 0 (n even, q = n/2) this is (x^(n/2))²; when r = 1 (n odd, q = (n−1)/2) this is x·(x^((n−1)/2))². This is exactly the identity stated in the Enuntiatio. The base case x^0 = 1 is the empty product. This justifies the recursive form; I now show the iterative form computes the same value.

2. Correctness by loop invariant. Let n have binary expansion n = Σ_{k=0}^{m} b_k·2^k with b_k ∈ {0,1} and m = ⌊log₂ n⌋ (for n > 0). The procedure consumes the bits b_0, b_1, … from least significant upward. Index the loop iterations by k = 0, 1, 2, …, and let e_k, base_k, result_k denote the values of e, base, result at the start of iteration k. I claim the invariant:

(I) at the start of iteration k, base_k = x^(2^k), e_k = ⌊n / 2^k⌋, and result_k = x^(Σ_{j<k} b_j·2^j).

Base (k = 0). Before the loop, base_0 = x = x^(2^0), e_0 = n = ⌊n/2^0⌋, and result_0 = 1 = x^(empty sum). The invariant holds.

Step. Suppose (I) holds at the start of iteration k and the loop body executes (so e_k > 0). The low bit of e_k = ⌊n/2^k⌋ is precisely b_k. If b_k = 1 the body multiplies result by base_k = x^(2^k), giving x^(Σ_{j<k} b_j 2^j) · x^(2^k) = x^(Σ_{j≤k} b_j 2^j); if b_k = 0 result is unchanged and equals x^(Σ_{j≤k} b_j 2^j) trivially (the new term contributes 2^k·0 = 0). Either way result_{k+1} = x^(Σ_{j<k+1} b_j 2^j). Next, base *= base gives base_{k+1} = (x^(2^k))² = x^(2^(k+1)). Finally e >>= 1 gives e_{k+1} = ⌊e_k/2⌋ = ⌊⌊n/2^k⌋/2⌋ = ⌊n/2^(k+1)⌋ (the standard identity for iterated integer division). Thus (I) holds at the start of iteration k+1.

Termination and conclusion. e is a non-negative integer strictly decreasing under e ↦ ⌊e/2⌋ while positive, so the loop halts; it halts at the first index K with e_K = ⌊n/2^K⌋ = 0, i.e. K = m+1 = ⌊log₂ n⌋ + 1 for n > 0. By (I) the returned value is result_K = x^(Σ_{j<K} b_j·2^j) = x^(Σ_{k=0}^{m} b_k·2^k) = x^n. For n = 0 the loop body never runs and the procedure returns the initial result = 1 = x^0. Correctness is established in all cases.

3. The multiplication bound. The loop executes exactly K = ⌊log₂ n⌋ + 1 iterations for n > 0 (and none for n = 0). Each iteration performs exactly one squaring (base *= base) and at most one accumulating multiplication (result *= base, only when b_k = 1). Hence the total number of multiplications is at most 2K = 2·(⌊log₂ n⌋ + 1), and at least K (one squaring per iteration). This is Θ(log n), which for all n ≥ 3 is strictly fewer than the n−1 multiplications of the naive successive-product method, and the gap widens without bound. The witness power(3, 13): n = 13 = 1101₂, so K = ⌊log₂ 13⌋ + 1 = 3 + 1 = 4 iterations, three set bits among b_3 b_2 b_1 b_0 = 1101 contributing three accumulations, for 4 + 3 = 7 multiplications against the naive 12 — and the result is 3^13 = 1594323, as asserted in the Expressio.

Q.E.D.