Learning square-and-multiply algorithm

Question

I've spent some time looking at various algorithms used for square-and-multiply techniques and I've found one that makes more sense to me than others. To put it to use, I am trying to compute the following example:

Compute (31^39 mod 773). 

Let a = 31, b = 39.

So the first step is to convert 39 to binary:

39 = 100111

Now, we look at each bit using the following rules:

If we encounter a 0, we square a. 
If we encounter a 1, we square a, then multiply by a.

Using these rules, the work looks like so:

Bit ------- Work
 1           Starts equation with number 
 0           31
 0           (31^2)
 1           ((31^2)^2 * 31)
 1           (((31^2)^2*31)^2 *31)
 1           (((((31^2)^2*31)^2*31)^2*31)

A side note, I'm not exactly clear on how the first 1 we encounter starts the equation with the number, and why it's blank initially but when we square it (first 0), we get our number.

I'm not sure where to go from here. I understand that what we've just done is break 31^39 down into a much more manageable number, but I don't see how the following helps us compute what we want (31^39 mod 773).

---

Brian and John, that makes complete sense, awesome, thank you! I went ahead and did the computations that looked liked this:

961 mod 773 = 188
188^2 mod 773 = 559
(559*31)^2 mod 773 = 747
(747*31)^2 mod 773 = 316
(316*31) mod 773 = 520

Thus, 31^39 mod 773 = 520.

Answer 1

When you first square 31, you get 961...which you can immediately reduce mod 773 to get 188. You square that...and reduce mod 773 again. And you find that you're never sqauring a number larger than 773, so your work becomes far more manageable.

This uses the idea that $$(ab \bmod n) = ((a \bmod n) (b \bmod n)) \bmod n.$$

Answer 2

It helps because

$$ ab \bmod c = (a \bmod c)(b \bmod c) \bmod c $$

In particular,

$$ (a^m)^n \bmod c = (a^m \bmod c)^n \bmod c $$

Therefore, after each step, you can take the intermediate result modulo $773$, thereby dealing with relatively small numbers (numbers always smaller than $773^2$), instead of constructing a very large number like $31^{39}$, and then taking the result modulo $773$.

Answer 3

For the general modular exponentiation problem $$x=a^b(\text{mod}\ m)$$ this question has, for a computer anyway, very small values for a,b,m.

Indeed it would make some sense to simply evaluate 31*31*31*...*31 (mod 773). but say I had, $$15236232847^{9050829380458}(\text{mod }10748092837)$$

not only would it require 9050829380458 multiplication operations but each multiplication operation would happen on the order of 9050829380458 digits.

instead we rewrite b as the sum $$\sum_{i=0}^{\log_2 b}d_i*2^i$$ where $d_i$ is the $i$-th digit in the binary representation of b (counting from the right), and $a^b$ becomes $$\Pi_{i=0}^{log_2 b}a^{d_i*2^i}$$ if all the $d_i$ are 1, then each term in this product is the square of the previous term, $a^{2^{i+1}}=a^{2*2^i}={a^{2^i}}^2$ and the entire product can be calculated in $\log_2\ {n}$ multiplications. The numbers that we are multiplying are still very large though, but this too is not necessary as $ab(\text{mod }m)=(a(\text{mod }m)*b(\text{mod }m))(\text{mod }m)$ and we can express every term modulo m after we square it. This ensures that that the numbers we are multiplying together are never larger than $m^2$, and the whole process takes $O(\log_2{b}*\log_2 m^2)$ operations (assuming each multiplication is performed in a process similar to the one we just described).

Answer 4

Let's say $d = (b_{k-1}b_{k-2}\ldots b_1 b_0)_2$, where $b_i \in \{0,1\}$ and $2^{k-1} \leq d < 2^k$

$\begin{align} a^d & \equiv a^{\sum\limits_{i=0}^k 2^ib_i} & (mod\; m) \\ & \equiv \prod\limits_{i=0}^k (a^{2^i})^{b_i} & (mod\; m) \end{align}$

This leads to the algorithm and the corresponding python implementation below:

Progressively square the base at each step to get factor starting with least significant bit of exponent

import math
def bin_rev(d):
    b = []
    while d:
        b.append(d % 2)
        d //= 2
    return b
def mult_exp_mod(a, d, m):
    a_mpow = [a % m]
    for i in range(int(math.log2(d))):
        a_mpow.append(a_mpow[-1]*2 % m)
    b, res = bin_rev(d), 1
    for i in range(len(b)):
        if b[i]:
            res = (res  a_mpow[i]) % m
    return res

For example, $47^{69} \equiv 125 \;(mod\; 143)$, which can be computed with the above function:

print(mult_exp_mod(47, 69, 143))
# 125