num_qubits: int, group_phase: bool = False, seed: int | np.random.Generator | None = None

):

"""Return a random Pauli.

Args:

num_qubits (int): the number of qubits.

group_phase (bool): Optional. If True generate random phase.

Otherwise the phase will be set so that the

Pauli coefficient is +1 (default: False).

seed (int or np.random.Generator): Optional. Set a fixed seed or

generator for RNG.

Returns:

Pauli: a random Pauli

"""

if seed is None:

rng = np.random.default_rng()

elif isinstance(seed, np.random.Generator):

rng = seed

else:

rng = default_rng(seed)

z = rng.integers(2, size=num_qubits, dtype=bool)

x = rng.integers(2, size=num_qubits, dtype=bool)

phase = rng.integers(4) if group_phase else 0

pauli = Pauli((z, x, phase))

num_qubits: int,

size: int = 1,

seed: int | np.random.Generator | None = None,

phase: bool = True,

):

"""Return a random PauliList.

Args:

num_qubits (int): the number of qubits.

size (int): Optional. The length of the Pauli list (Default: 1).

seed (int or np.random.Generator): Optional. Set a fixed seed or generator for RNG.

phase (bool): If True the Pauli phases are randomized, otherwise the phases are fixed to 0.

[Default: True]

Returns:

PauliList: a random PauliList.

"""

if seed is None:

rng = np.random.default_rng()

elif isinstance(seed, np.random.Generator):

rng = seed

else:

rng = default_rng(seed)

z = rng.integers(2, size=(size, num_qubits)).astype(bool)

x = rng.integers(2, size=(size, num_qubits)).astype(bool)

if phase:

_phase = rng.integers(4, size=(size))

return PauliList.from_symplectic(z, x, _phase)

"""Return a random Clifford operator.

The Clifford is sampled using the method of Reference [1].

Args:

num_qubits (int): the number of qubits for the Clifford

seed (int or np.random.Generator): Optional. Set a fixed seed or

generator for RNG.

Returns:

Clifford: a random Clifford operator.

Reference:

1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the

structure of the Clifford group*.

`arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_

"""

if seed is None:

rng = np.random.default_rng()

elif isinstance(seed, np.random.Generator):

rng = seed

else:

rng = default_rng(seed)

had, perm = _sample_qmallows(num_qubits, rng)

gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))

gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))

delta1 = np.eye(num_qubits, dtype=np.int8)

delta2 = delta1.copy()

_fill_tril(gamma1, rng, symmetric=True)

_fill_tril(gamma2, rng, symmetric=True)

_fill_tril(delta1, rng)

_fill_tril(delta2, rng)

# For large num_qubits numpy.inv function called below can

# return invalid output leading to a non-symplectic Clifford

# being generated. This can be prevented by manually forcing

# block inversion of the matrix.

block_inverse_threshold = 50

# Compute stabilizer table

zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)

prod1 = np.matmul(gamma1, delta1) % 2

prod2 = np.matmul(gamma2, delta2) % 2

inv1 = _inverse_tril(delta1, block_inverse_threshold).transpose()

inv2 = _inverse_tril(delta2, block_inverse_threshold).transpose()

table1 = np.block([[delta1, zero], [prod1, inv1]])

table2 = np.block([[delta2, zero], [prod2, inv2]])

# Apply qubit permutation

table = table2[np.concatenate([perm, num_qubits + perm])]

# Apply layer of Hadamards

inds = had * np.arange(1, num_qubits + 1)

inds = inds[inds > 0] - 1

lhs_inds = np.concatenate([inds, inds + num_qubits])

rhs_inds = np.concatenate([inds + num_qubits, inds])

table[lhs_inds, :] = table[rhs_inds, :]

# Apply table

tableau = np.zeros((2 * num_qubits, 2 * num_qubits + 1), dtype=bool)

tableau[:, :-1] = np.mod(np.matmul(table1, table), 2)

# Generate random phases

tableau[:, -1] = rng.integers(2, size=2 * num_qubits)

"""Sample from the quantum Mallows distribution"""

if rng is None:

rng = np.random.default_rng()

# Hadmard layer

had = np.zeros(n, dtype=bool)

# Permutation layer

perm = np.zeros(n, dtype=int)

inds = list(range(n))

for i in range(n):

m = n - i

eps = 4 ** (-m)

r = rng.uniform(0, 1)

index = -math.ceil(math.log2(r + (1 - r) * eps))

had[i] = index < m

if index < m:

k = index

else:

k = 2 * m - index - 1

perm[i] = inds[k]

del inds[k]

"""Add symmetric random ints to off diagonals"""

dim = mat.shape[0]

# Optimized for low dimensions

if dim == 1:

return

if dim <= 4:

mat[1, 0] = rng.integers(2, dtype=np.int8)

if symmetric:

mat[0, 1] = mat[1, 0]

if dim > 2:

mat[2, 0] = rng.integers(2, dtype=np.int8)

mat[2, 1] = rng.integers(2, dtype=np.int8)

if symmetric:

mat[0, 2] = mat[2, 0]

mat[1, 2] = mat[2, 1]

if dim > 3:

mat[3, 0] = rng.integers(2, dtype=np.int8)

mat[3, 1] = rng.integers(2, dtype=np.int8)

mat[3, 2] = rng.integers(2, dtype=np.int8)

if symmetric:

mat[0, 3] = mat[3, 0]

mat[1, 3] = mat[3, 1]

mat[2, 3] = mat[3, 2]

return

# Use numpy indices for larger dimensions

rows, cols = np.tril_indices(dim, -1)

vals = rng.integers(2, size=rows.size, dtype=np.int8)

mat[(rows, cols)] = vals

if symmetric:

"""Invert a lower-triangular matrix with unit diagonal."""

# Optimized inversion function for low dimensions

dim = mat.shape[0]

if dim <= 2:

return mat

if dim <= 5:

inv = mat.copy()

inv[2, 0] = mat[2, 0] ^ (mat[1, 0] & mat[2, 1])

if dim > 3:

inv[3, 1] = mat[3, 1] ^ (mat[2, 1] & mat[3, 2])

inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])

if dim > 4:

inv[4, 2] = (mat[4, 2] ^ (mat[3, 2] & mat[4, 3])) & 1

inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])

inv[4, 0] = (

mat[4, 0]

^ (mat[1, 0] & inv[4, 1])

^ (mat[2, 0] & inv[4, 2])

^ (mat[3, 0] & mat[4, 3])

)

return inv % 2

# For higher dimensions we use Numpy's inverse function

# however this function tends to fail and result in a non-symplectic

# final matrix if n is too large.

if dim <= block_inverse_threshold:

return np.linalg.inv(mat).astype(np.int8) % 2

# For very large matrices we divide the matrix into 4 blocks of

# roughly equal size and use the analytic formula for the inverse

# of a block lower-triangular matrix:

# inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]

# call the inverse function recursively to compute inv(A) and invD

dim1 = dim // 2

mat_a = _inverse_tril(mat[0:dim1, 0:dim1], block_inverse_threshold)

mat_d = _inverse_tril(mat[dim1:dim, dim1:dim], block_inverse_threshold)

mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)

inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])