import math

def shape(matrix):
    if not matrix or not matrix[0] or any(len(row) != len(matrix[0]) for row in matrix):
        raise ValueError("matrix must be nonempty and rectangular")
    return len(matrix), len(matrix[0])

def matvec(matrix, vector):
    rows, columns = shape(matrix)
    if columns != len(vector):
        raise ValueError("matrix/vector dimension mismatch")
    return tuple(sum(matrix[row][column] * vector[column] for column in range(columns))
                 for row in range(rows))

def matmul(left, right):
    left_rows, shared = shape(left)
    right_rows, columns = shape(right)
    if shared != right_rows:
        raise ValueError("matrix product dimension mismatch")
    return tuple(tuple(sum(left[row][k] * right[k][column] for k in range(shared))
                       for column in range(columns)) for row in range(left_rows))

def max_error(left, right):
    return max(abs(a - b) for row_a, row_b in zip(left, right) for a, b in zip(row_a, row_b))

scale = math.sqrt(0.5)
I = ((1, 0), (0, 1))
X = ((0, 1), (1, 0))
Y = ((0, -1j), (1j, 0))
Z = ((1, 0), (0, -1))
H = ((scale, scale), (scale, -scale))
S = ((1, 0), (0, 1j))
gates = (X, Y, Z, H, S)

identities = (matmul(H, H), matmul(matmul(H, Z), H), matmul(matmul(H, X), H))
targets = (I, X, Z)
identity_error = max(max_error(actual, expected) for actual, expected in zip(identities, targets))
assert identity_error < 1e-12

basis = ((1, 0), (0, 1))
assert matvec(X, basis[0]) == basis[1] and matvec(X, basis[1]) == basis[0]
assert matvec(Z, basis[0]) == basis[0] and matvec(Z, basis[1]) == (0, -1)
assert matvec(Y, basis[0]) == (0j, 1j)

states = ((1, 0), (0, 1), (scale, scale), (0.6, 0.8j))
composition_error = 0.0
for left, right in ((H, Z), (S, H), (X, S), (Y, H)):
    composed = matmul(left, right)
    for state in states:
        sequential = matvec(left, matvec(right, state))
        direct = matvec(composed, state)
        composition_error = max(composition_error, *(abs(a - b) for a, b in zip(sequential, direct)))
assert composition_error < 1e-12

witness = (scale, 1j * scale)
assert any(abs(a - b) > 1e-6 for a, b in zip(matvec(matmul(S, H), witness), matvec(matmul(H, S), witness)))
ragged_rejected = False
try:
    matvec(((1, 0), (0,)), (1, 0))
except ValueError:
    ragged_rejected = True
assert ragged_rejected
print(f"PASS: 11 gate bench verifies {len(gates)} reference matrices and {4 * len(states)} compositions; max identity error={identity_error:.2e}, sequential error={composition_error:.2e}")