A python 3 program that print the operation matrix for a CROT operation from the paper at http://arxiv.org/abs/1510.00409 . The printed operation is just a controlled quarter-turn around the Y axis.

walters_bqp_np_CROT.py

# -*- coding: utf-8 -*-
import numpy as np
import math
import cmath

np.set_printoptions(precision=2, suppress=True)
qubit_count = 2
σx = np.mat([[0, 1], [1, 0]])
σy = np.mat([[0, -1j], [1j, 0]])
σz = np.mat([[1, 0], [0, -1]])
π = math.pi

def CROT(c, t, a, θ):
    return R(t, θ/2) * Rc(c, t, a, θ/2, True)

def Rc(c, t, a, θ, on):
    mexp = matrix_lift(lambda λ: cmath.exp((-1 if on else 1) * λ * 1j * θ / 2))
    op = expand(σy, t)*expand(a, c)
    return mexp(op)

def R(t, θ):
    mexp = matrix_lift(lambda λ: cmath.exp(λ * 1j * θ / 2))
    op = expand(σy, t)
    return mexp(op)

def expand(qubit_op_matrix, qubit_index):
    """
    Expands a single-qubit operation to apply to one of the qubits amongst several.
    """
    post = np.identity(1 << qubit_index)
    pre = np.identity(1 << (qubit_count - qubit_index - 1))
    return np.kron(np.kron(pre, qubit_op_matrix), post)

def matrix_lift(f):
    """
    Lifts a function so it modifies matrices instead of numbers, by taking the
    spectral decomposition of the matrix and using the function to transform
    its eigenvalues.
    """

    def matrix_lift_helper(m):
        w, v = np.linalg.eig(m)
        result = np.mat(np.zeros(m.shape, np.complex128))
        for i in range(len(w)):
            eigen_val = w[i]
            eigen_vec = np.mat(v[:, i])
            eigen_mat = np.dot(eigen_vec, eigen_vec.H)
            result += f(eigen_val) * eigen_mat
        return result


    return lambda m: matrix_lift_helper(m)

print(CROT(1, 0, σz, π/2))