Binary Exponentiation by Factoring

This is a specialised trick, worth knowing as an example of how group structure yields algorithms. The task:

compute $a\,x^y \pmod{2^d}$ for an odd $x$, with $d \ge 3$ (typically $d = 32$ or $64$: the natural wrap-around arithmetic of machine integers).

Binary exponentiation needs $O(\log y)$ multiplications. Using the structure of the multiplicative group modulo $2^d$, we can do it with $O(d)$ additions and bit operations and a single multiplication by $y$.

The group modulo $2^d$

Every number $x \equiv 1 \pmod 4$ can be written as a power of a fixed base $b$ with $b \equiv 5 \pmod 8$:

The exponent $L(x)$ is a discrete logarithm. Then $a\,x^y \equiv a\,b^{\,y L(x)}$: we need $L(x)$ (a logarithm), one multiplication by $y$, and then an exponential $b^{(\cdot)}$. The trick is to make both of these cheap.

If $x \equiv 3 \pmod 4$ replace $x \mapsto -x$ and $a \mapsto (-1)^y a$, since $-x \equiv 1 \pmod 4$ and $x^y = (-1)^y(-x)^y$.

We work with $4L(x)$ (modulo $2^d$) rather than $L(x)$; it makes the bit structure of the table cleaner.

Computing $4L(x)$ by peeling off factors $2^n + 1$

Every $x \equiv 1 \pmod 4$ factors as

so $4L(x) = \sum_i 4L(2^{a_i} + 1)$. If we precompute the table $t_n = 4L(2^n + 1)$ for $1 < n < d$, then $4L(x)$ needs at most $d$ table additions.

To find the factors: repeatedly multiply $x$ by $2^n + 1$ for those bits $n$ (increasing) at which $x$ currently has a 1. Since $x \cdot (2^n + 1) = x + (x \ll n)$ and $x$ is odd, this clears bit $n$ and leaves the lower bits untouched. After the loop $x \cdot \prod (2^{n}+1) \equiv 1$, so we found the factorization of $x^{-1}$, and $4L(x) = -\sum t_n$.

The table

Computing $L$ for the numbers $2^n + 1$ can be done by a bit-by-bit search (a specialisation of the Pohlig-Hellman algorithm): the group generated by $b$ has order $2^{d-2}$, so we can determine the bits of $L$ from the lowest one upward by checking whether a suitable power of the remaining quantity is $1$.

D = 32
MASK = (1 << D) - 1

def dlog5(y, d=D):
    """Exponent L with 5^L = y (mod 2^d), for y = 1 (mod 4). Found bit by bit."""
    M = 1 << d
    L, cur = 0, y % M
    for i in range(d - 2):
        if pow(cur, 1 << (d - 3 - i), M) != 1:          # bit i of L is set
            L |= 1 << i
            cur = cur * pow(pow(5, 1 << i, M), -1, M) % M
    return L

assert all(pow(5, dlog5(y), 1 << D) == y for y in (5, 9, 17, 33, 1 + (1 << 20), 0x12345671))

# The published algorithm uses the base b = 0x1998df85 (b = 5 mod 8) instead of 5.
B = 0x1998DF85
L5_of_B = dlog5(B)                                       # b = 5^L5_of_B, and L_b(y) = L_5(y) / L_5(b)
inv = pow(L5_of_B, -1, 1 << (D - 2))
TABLE = [0, 0] + [(4 * dlog5((1 << n) + 1) * inv) & MASK for n in range(2, D)]

published = [
    0x00000000, 0x00000000, 0xd3cfd984, 0x9ee62e18, 0xe83d9070, 0xb59e81e0, 0xa17407c0, 0xce601f80,
    0xf4807f00, 0xe701fe00, 0xbe07fc00, 0xfc1ff800, 0xf87ff000, 0xf1ffe000, 0xe7ffc000, 0xdfff8000,
    0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000, 0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
    0xff000000, 0xfe000000, 0xfc000000, 0xf8000000, 0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
]
assert TABLE == published                     # we regenerated the table from first principles

Our independently computed table equals the constants published in the original article.

Note the pattern in the upper half: for $2n \ge d$ we have $(2^n+1)^2 \equiv 2^{n+1} + 1 \pmod{2^d}$, and the table entries there are simply $-2^n$ (the constants 0xffff0000, 0xfffe0000, ... are $-2^{16}, -2^{17}, \dots$ in two's complement):

assert all(TABLE[n] == (-(1 << n)) & MASK for n in range(16, 32))

The functions

def mbin_log_32(r, x):
    """r + 4*L(x) mod 2^32, for x = 1 (mod 4)."""
    for n in range(2, D):
        if x & (1 << n):
            x = (x + (x << n)) & MASK
            r = (r - TABLE[n]) & MASK
    return r

def mbin_exp_32(r, x):
    """r * b^(x/4) mod 2^32 (x a multiple of 4 given as a 4L value)."""
    for n in range(2, D):
        if x & (1 << n):
            r = (r + (r << n)) & MASK
            x = (x - TABLE[n]) & MASK
    return r

def mbin_power_odd_32(rem, base, exp):
    """rem * base^exp mod 2^32 for odd base."""
    if base & 2:                                         # base = 3 (mod 4): use -base
        base = (-base) & MASK
        if exp & 1:
            rem = (-rem) & MASK
    return mbin_exp_32(rem, (mbin_log_32(0, base) * exp) & MASK)

import random
random.seed(1)
for _ in range(2000):
    a = random.getrandbits(32)
    x = random.getrandbits(32) | 1
    y = random.getrandbits(random.randint(1, 40))
    assert mbin_power_odd_32(a, x, y) == a * pow(x, y, 1 << 32) % (1 << 32)
assert mbin_power_odd_32(1, 3, 5) == 243
assert mbin_power_odd_32(7, 0xFFFFFFFF, 3) == (7 * pow(0xFFFFFFFF, 3, 1 << 32)) & MASK

The loops contain no multiplications: only additions, shifts and mask operations, plus the single product mbin_log_32(0, base) * exp.

Halving the loop

Because the table entries for $2n \ge d$ are just $-2^n$, the loop can stop at $n = d/2$; the remaining high bits of $x$ are handled with a single subtraction (for the logarithm) or a single multiplication (for the exponential) instead of $d/2$ more iterations:

def mbin_log_32_fast(r, x):
    for n in range(2, 16):
        if x & (1 << n):
            x = (x + (x << n)) & MASK
            r = (r - TABLE[n]) & MASK
    return (r - (x & 0xFFFF0000)) & MASK

def mbin_exp_32_fast(r, x):
    for n in range(2, 16):
        if x & (1 << n):
            r = (r + (r << n)) & MASK
            x = (x - TABLE[n]) & MASK
    return (r * (1 - (x & 0xFFFF0000))) & MASK

for _ in range(2000):
    x = random.getrandbits(32) & ~2 | 1                  # x = 1 (mod 4)
    assert mbin_log_32_fast(0, x) == mbin_log_32(0, x)
    e = random.getrandbits(30) << 2
    assert mbin_exp_32_fast(1, e) == mbin_exp_32(1, e)

The exact same code with 64-bit masks and a 64-entry table gives the 64-bit version.

Why it matters

The technique shows that the structure of the multiplicative group (cyclic times $\pm 1$ for powers of two) can turn exponentiation into arithmetic on exponents. The same ideas underlie the discrete logarithm algorithms and primitive roots. Direct uses: fast hash functions with multiplicative mixing modulo $2^{64}$, deterministic pseudo-random jump-ahead for linear congruential generators, and proofs about $x^y \bmod 2^d$.