Matrix for modular exponentiation in Shor's algorithm
Note that Shor's algorithm doesn't require modular exponentiation. All you need is (controlled) modular multiplication.

There are two steps to getting the matrix for modular multiplication:


Write a function to turn any permutation into a matrix.
Apply it to the permutation $a \rightarrow a \cdot k \pmod{r}$.
(Optional) Instead of making that matrix and multiplying by it, save massive execution time by just applying the permutation directly.


For step 1, just have a function like this:

import numpy as np

def make_permutation_matrix(n, permutation):
    r = np.zeros((n, n))
    for i in range(n):
        r[permutation(i), i] = 1
    return r


Then for step 2, use it:

import math
k = pow(base, exponent, modulus)
n = math.ceil(math.log(modulus, 2))
mul_matrix = make_permutation_matrix(
    2**n,
    permutation=lambda x: x*k % modulus)