We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Randomized testing failure on main
<!-- commit 60c8bba11f7314e376ff2fd7e6f211bad11b702e@master -->
Trying to build `master` at commit 60c8bba11f7314e376ff2fd7e6f211bad11b702e failed.
More info at: https://travis-ci.com/Qiskit/qiskit-terra/jobs/208623536

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `qiskit/quantum_info/operators/symplectic/random.py`
Content:
```
1 # -*- coding: utf-8 -*-
2 
3 # This code is part of Qiskit.
4 #
5 # (C) Copyright IBM 2017, 2020
6 #
7 # This code is licensed under the Apache License, Version 2.0. You may
8 # obtain a copy of this license in the LICENSE.txt file in the root directory
9 # of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
10 #
11 # Any modifications or derivative works of this code must retain this
12 # copyright notice, and modified files need to carry a notice indicating
13 # that they have been altered from the originals.
14 """
15 Random symplectic operator functions
16 """
17 
18 import numpy as np
19 from numpy.random import default_rng
20 
21 from .clifford import Clifford
22 from .stabilizer_table import StabilizerTable
23 from .pauli_table import PauliTable
24 
25 
26 def random_pauli_table(num_qubits, size=1, seed=None):
27     """Return a random PauliTable.
28 
29     Args:
30         num_qubits (int): the number of qubits.
31         size (int): Optional. The number of rows of the table (Default: 1).
32         seed (int or np.random.Generator): Optional. Set a fixed seed or
33                                            generator for RNG.
34 
35     Returns:
36         PauliTable: a random PauliTable.
37     """
38     if seed is None:
39         rng = np.random.default_rng()
40     elif isinstance(seed, np.random.Generator):
41         rng = seed
42     else:
43         rng = default_rng(seed)
44 
45     table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)
46     return PauliTable(table)
47 
48 
49 def random_stabilizer_table(num_qubits, size=1, seed=None):
50     """Return a random StabilizerTable.
51 
52     Args:
53         num_qubits (int): the number of qubits.
54         size (int): Optional. The number of rows of the table (Default: 1).
55         seed (int or np.random.Generator): Optional. Set a fixed seed or
56                                            generator for RNG.
57 
58     Returns:
59         PauliTable: a random StabilizerTable.
60     """
61     if seed is None:
62         rng = np.random.default_rng()
63     elif isinstance(seed, np.random.Generator):
64         rng = seed
65     else:
66         rng = default_rng(seed)
67 
68     table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)
69     phase = rng.integers(2, size=size).astype(np.bool)
70     return StabilizerTable(table, phase)
71 
72 
73 def random_clifford(num_qubits, seed=None):
74     """Return a random Clifford operator.
75 
76     The Clifford is sampled using the method of Reference [1].
77 
78     Args:
79         num_qubits (int): the number of qubits for the Clifford
80         seed (int or np.random.Generator): Optional. Set a fixed seed or
81                                            generator for RNG.
82 
83     Returns:
84         Clifford: a random Clifford operator.
85 
86     Reference:
87         1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the
88            structure of the Clifford group*.
89            `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_
90     """
91     if seed is None:
92         rng = np.random.default_rng()
93     elif isinstance(seed, np.random.Generator):
94         rng = seed
95     else:
96         rng = default_rng(seed)
97 
98     had, perm = _sample_qmallows(num_qubits, rng)
99 
100     gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))
101     gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))
102     delta1 = np.eye(num_qubits, dtype=np.int8)
103     delta2 = delta1.copy()
104 
105     _fill_tril(gamma1, rng, symmetric=True)
106     _fill_tril(gamma2, rng, symmetric=True)
107     _fill_tril(delta1, rng)
108     _fill_tril(delta2, rng)
109 
110     # Compute stabilizer table
111     zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)
112     prod1 = np.matmul(gamma1, delta1) % 2
113     prod2 = np.matmul(gamma2, delta2) % 2
114     inv1 = _inverse_tril(delta1).transpose()
115     inv2 = _inverse_tril(delta2).transpose()
116     table1 = np.block([[delta1, zero], [prod1, inv1]])
117     table2 = np.block([[delta2, zero], [prod2, inv2]])
118 
119     # Apply qubit permutation
120     table = table2[np.concatenate([perm, num_qubits + perm])]
121 
122     # Apply layer of Hadamards
123     inds = had * np.arange(1, num_qubits + 1)
124     inds = inds[inds > 0] - 1
125     lhs_inds = np.concatenate([inds, inds + num_qubits])
126     rhs_inds = np.concatenate([inds + num_qubits, inds])
127     table[lhs_inds, :] = table[rhs_inds, :]
128 
129     # Apply table
130     table = np.mod(np.matmul(table1, table), 2).astype(np.bool)
131 
132     # Generate random phases
133     phase = rng.integers(2, size=2 * num_qubits).astype(np.bool)
134     return Clifford(StabilizerTable(table, phase))
135 
136 
137 def _sample_qmallows(n, rng=None):
138     """Sample from the quantum Mallows distribution"""
139 
140     if rng is None:
141         rng = np.random.default_rng()
142 
143     # Hadmard layer
144     had = np.zeros(n, dtype=np.bool)
145 
146     # Permutation layer
147     perm = np.zeros(n, dtype=int)
148 
149     inds = list(range(n))
150     for i in range(n):
151         m = n - i
152         eps = 4 ** (-m)
153         r = rng.uniform(0, 1)
154         index = -int(np.ceil(np.log2(r + (1 - r) * eps)))
155         had[i] = index < m
156         if index < m:
157             k = index
158         else:
159             k = 2 * m - index - 1
160         perm[i] = inds[k]
161         del inds[k]
162     return had, perm
163 
164 
165 def _fill_tril(mat, rng, symmetric=False):
166     """Add symmetric random ints to off diagonals"""
167     dim = mat.shape[0]
168     # Optimized for low dimensions
169     if dim == 1:
170         return
171 
172     if dim <= 4:
173         mat[1, 0] = rng.integers(2, dtype=np.int8)
174         if symmetric:
175             mat[0, 1] = mat[1, 0]
176         if dim > 2:
177             mat[2, 0] = rng.integers(2, dtype=np.int8)
178             mat[2, 1] = rng.integers(2, dtype=np.int8)
179             if symmetric:
180                 mat[0, 2] = mat[2, 0]
181                 mat[1, 2] = mat[2, 1]
182         if dim > 3:
183             mat[3, 0] = rng.integers(2, dtype=np.int8)
184             mat[3, 1] = rng.integers(2, dtype=np.int8)
185             mat[3, 2] = rng.integers(2, dtype=np.int8)
186             if symmetric:
187                 mat[0, 3] = mat[3, 0]
188                 mat[1, 3] = mat[3, 1]
189                 mat[2, 3] = mat[3, 2]
190         return
191 
192     # Use numpy indices for larger dimensions
193     rows, cols = np.tril_indices(dim, -1)
194     vals = rng.integers(2, size=rows.size, dtype=np.int8)
195     mat[(rows, cols)] = vals
196     if symmetric:
197         mat[(cols, rows)] = vals
198 
199 
200 def _inverse_tril(mat):
201     """Invert a lower-triangular matrix with unit diagonal."""
202     # Optimized inversion function for low dimensions
203     dim = mat.shape[0]
204 
205     if dim <= 2:
206         return mat
207 
208     if dim <= 5:
209         inv = mat.copy()
210         inv[2, 0] = (mat[2, 0] ^ (mat[1, 0] & mat[2, 1]))
211         if dim > 3:
212             inv[3, 1] = (mat[3, 1] ^ (mat[2, 1] & mat[3, 2]))
213             inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])
214         if dim > 4:
215             inv[4, 2] = ((mat[4, 2] ^ (mat[3, 2] & mat[4, 3]))) & 1
216             inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])
217             inv[4, 0] = mat[4, 0] ^ (mat[1, 0] & inv[4, 1]) ^ (
218                 mat[2, 0] & inv[4, 2]) ^ (mat[3, 0] & mat[4, 3])
219         return inv % 2
220 
221     # For higher dimensions we use Numpy's inverse function
222     # however this function tends to fail and result in a non-symplectic
223     # final matrix if n is too large.
224     max_np_inv = 150
225     if dim <= max_np_inv:
226         return np.linalg.inv(mat).astype(np.int8) % 2
227 
228     # For very large matrices  we divide the matrix into 4 blocks of
229     # roughly equal size and use the analytic formula for the inverse
230     # of a block lower-triangular matrix:
231     # inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]
232     # call the inverse function recursively to compute inv(A) and invD
233 
234     dim1 = dim // 2
235     mat_a = _inverse_tril(mat[0:dim1, 0:dim1])
236     mat_d = _inverse_tril(mat[dim1:dim, dim1:dim])
237     mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)
238     inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])
239     return inv % 2
240
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/qiskit/quantum_info/operators/symplectic/random.py b/qiskit/quantum_info/operators/symplectic/random.py
--- a/qiskit/quantum_info/operators/symplectic/random.py
+++ b/qiskit/quantum_info/operators/symplectic/random.py
@@ -107,12 +107,18 @@
     _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).transpose()
-    inv2 = _inverse_tril(delta2).transpose()
+    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]])
 
@@ -197,7 +203,7 @@
         mat[(cols, rows)] = vals
 
 
-def _inverse_tril(mat):
+def _inverse_tril(mat, block_inverse_threshold):
     """Invert a lower-triangular matrix with unit diagonal."""
     # Optimized inversion function for low dimensions
     dim = mat.shape[0]
@@ -221,8 +227,7 @@
     # 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.
-    max_np_inv = 150
-    if dim <= max_np_inv:
+    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
@@ -232,8 +237,8 @@
     # call the inverse function recursively to compute inv(A) and invD
 
     dim1 = dim // 2
-    mat_a = _inverse_tril(mat[0:dim1, 0:dim1])
-    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim])
+    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]])
     return inv % 2

{"golden_diff": "diff --git a/qiskit/quantum_info/operators/symplectic/random.py b/qiskit/quantum_info/operators/symplectic/random.py\n--- a/qiskit/quantum_info/operators/symplectic/random.py\n+++ b/qiskit/quantum_info/operators/symplectic/random.py\n@@ -107,12 +107,18 @@\n     _fill_tril(delta1, rng)\n     _fill_tril(delta2, rng)\n \n+    # For large num_qubits numpy.inv function called below can\n+    # return invalid output leading to a non-symplectic Clifford\n+    # being generated. This can be prevented by manually forcing\n+    # block inversion of the matrix.\n+    block_inverse_threshold = 50\n+\n     # Compute stabilizer table\n     zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)\n     prod1 = np.matmul(gamma1, delta1) % 2\n     prod2 = np.matmul(gamma2, delta2) % 2\n-    inv1 = _inverse_tril(delta1).transpose()\n-    inv2 = _inverse_tril(delta2).transpose()\n+    inv1 = _inverse_tril(delta1, block_inverse_threshold).transpose()\n+    inv2 = _inverse_tril(delta2, block_inverse_threshold).transpose()\n     table1 = np.block([[delta1, zero], [prod1, inv1]])\n     table2 = np.block([[delta2, zero], [prod2, inv2]])\n \n@@ -197,7 +203,7 @@\n         mat[(cols, rows)] = vals\n \n \n-def _inverse_tril(mat):\n+def _inverse_tril(mat, block_inverse_threshold):\n     \"\"\"Invert a lower-triangular matrix with unit diagonal.\"\"\"\n     # Optimized inversion function for low dimensions\n     dim = mat.shape[0]\n@@ -221,8 +227,7 @@\n     # For higher dimensions we use Numpy's inverse function\n     # however this function tends to fail and result in a non-symplectic\n     # final matrix if n is too large.\n-    max_np_inv = 150\n-    if dim <= max_np_inv:\n+    if dim <= block_inverse_threshold:\n         return np.linalg.inv(mat).astype(np.int8) % 2\n \n     # For very large matrices  we divide the matrix into 4 blocks of\n@@ -232,8 +237,8 @@\n     # call the inverse function recursively to compute inv(A) and invD\n \n     dim1 = dim // 2\n-    mat_a = _inverse_tril(mat[0:dim1, 0:dim1])\n-    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim])\n+    mat_a = _inverse_tril(mat[0:dim1, 0:dim1], block_inverse_threshold)\n+    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim], block_inverse_threshold)\n     mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)\n     inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])\n     return inv % 2\n", "issue": "Randomized testing failure on main\n<!-- commit 60c8bba11f7314e376ff2fd7e6f211bad11b702e@master -->\nTrying to build `master` at commit 60c8bba11f7314e376ff2fd7e6f211bad11b702e failed.\nMore info at: https://travis-ci.com/Qiskit/qiskit-terra/jobs/208623536\n", "before_files": [{"content": "# -*- coding: utf-8 -*-\n\n# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2017, 2020\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n\"\"\"\nRandom symplectic operator functions\n\"\"\"\n\nimport numpy as np\nfrom numpy.random import default_rng\n\nfrom .clifford import Clifford\nfrom .stabilizer_table import StabilizerTable\nfrom .pauli_table import PauliTable\n\n\ndef random_pauli_table(num_qubits, size=1, seed=None):\n    \"\"\"Return a random PauliTable.\n\n    Args:\n        num_qubits (int): the number of qubits.\n        size (int): Optional. The number of rows of the table (Default: 1).\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        PauliTable: a random PauliTable.\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)\n    return PauliTable(table)\n\n\ndef random_stabilizer_table(num_qubits, size=1, seed=None):\n    \"\"\"Return a random StabilizerTable.\n\n    Args:\n        num_qubits (int): the number of qubits.\n        size (int): Optional. The number of rows of the table (Default: 1).\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        PauliTable: a random StabilizerTable.\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)\n    phase = rng.integers(2, size=size).astype(np.bool)\n    return StabilizerTable(table, phase)\n\n\ndef random_clifford(num_qubits, seed=None):\n    \"\"\"Return a random Clifford operator.\n\n    The Clifford is sampled using the method of Reference [1].\n\n    Args:\n        num_qubits (int): the number of qubits for the Clifford\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        Clifford: a random Clifford operator.\n\n    Reference:\n        1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the\n           structure of the Clifford group*.\n           `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    had, perm = _sample_qmallows(num_qubits, rng)\n\n    gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))\n    gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))\n    delta1 = np.eye(num_qubits, dtype=np.int8)\n    delta2 = delta1.copy()\n\n    _fill_tril(gamma1, rng, symmetric=True)\n    _fill_tril(gamma2, rng, symmetric=True)\n    _fill_tril(delta1, rng)\n    _fill_tril(delta2, rng)\n\n    # Compute stabilizer table\n    zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)\n    prod1 = np.matmul(gamma1, delta1) % 2\n    prod2 = np.matmul(gamma2, delta2) % 2\n    inv1 = _inverse_tril(delta1).transpose()\n    inv2 = _inverse_tril(delta2).transpose()\n    table1 = np.block([[delta1, zero], [prod1, inv1]])\n    table2 = np.block([[delta2, zero], [prod2, inv2]])\n\n    # Apply qubit permutation\n    table = table2[np.concatenate([perm, num_qubits + perm])]\n\n    # Apply layer of Hadamards\n    inds = had * np.arange(1, num_qubits + 1)\n    inds = inds[inds > 0] - 1\n    lhs_inds = np.concatenate([inds, inds + num_qubits])\n    rhs_inds = np.concatenate([inds + num_qubits, inds])\n    table[lhs_inds, :] = table[rhs_inds, :]\n\n    # Apply table\n    table = np.mod(np.matmul(table1, table), 2).astype(np.bool)\n\n    # Generate random phases\n    phase = rng.integers(2, size=2 * num_qubits).astype(np.bool)\n    return Clifford(StabilizerTable(table, phase))\n\n\ndef _sample_qmallows(n, rng=None):\n    \"\"\"Sample from the quantum Mallows distribution\"\"\"\n\n    if rng is None:\n        rng = np.random.default_rng()\n\n    # Hadmard layer\n    had = np.zeros(n, dtype=np.bool)\n\n    # Permutation layer\n    perm = np.zeros(n, dtype=int)\n\n    inds = list(range(n))\n    for i in range(n):\n        m = n - i\n        eps = 4 ** (-m)\n        r = rng.uniform(0, 1)\n        index = -int(np.ceil(np.log2(r + (1 - r) * eps)))\n        had[i] = index < m\n        if index < m:\n            k = index\n        else:\n            k = 2 * m - index - 1\n        perm[i] = inds[k]\n        del inds[k]\n    return had, perm\n\n\ndef _fill_tril(mat, rng, symmetric=False):\n    \"\"\"Add symmetric random ints to off diagonals\"\"\"\n    dim = mat.shape[0]\n    # Optimized for low dimensions\n    if dim == 1:\n        return\n\n    if dim <= 4:\n        mat[1, 0] = rng.integers(2, dtype=np.int8)\n        if symmetric:\n            mat[0, 1] = mat[1, 0]\n        if dim > 2:\n            mat[2, 0] = rng.integers(2, dtype=np.int8)\n            mat[2, 1] = rng.integers(2, dtype=np.int8)\n            if symmetric:\n                mat[0, 2] = mat[2, 0]\n                mat[1, 2] = mat[2, 1]\n        if dim > 3:\n            mat[3, 0] = rng.integers(2, dtype=np.int8)\n            mat[3, 1] = rng.integers(2, dtype=np.int8)\n            mat[3, 2] = rng.integers(2, dtype=np.int8)\n            if symmetric:\n                mat[0, 3] = mat[3, 0]\n                mat[1, 3] = mat[3, 1]\n                mat[2, 3] = mat[3, 2]\n        return\n\n    # Use numpy indices for larger dimensions\n    rows, cols = np.tril_indices(dim, -1)\n    vals = rng.integers(2, size=rows.size, dtype=np.int8)\n    mat[(rows, cols)] = vals\n    if symmetric:\n        mat[(cols, rows)] = vals\n\n\ndef _inverse_tril(mat):\n    \"\"\"Invert a lower-triangular matrix with unit diagonal.\"\"\"\n    # Optimized inversion function for low dimensions\n    dim = mat.shape[0]\n\n    if dim <= 2:\n        return mat\n\n    if dim <= 5:\n        inv = mat.copy()\n        inv[2, 0] = (mat[2, 0] ^ (mat[1, 0] & mat[2, 1]))\n        if dim > 3:\n            inv[3, 1] = (mat[3, 1] ^ (mat[2, 1] & mat[3, 2]))\n            inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])\n        if dim > 4:\n            inv[4, 2] = ((mat[4, 2] ^ (mat[3, 2] & mat[4, 3]))) & 1\n            inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])\n            inv[4, 0] = mat[4, 0] ^ (mat[1, 0] & inv[4, 1]) ^ (\n                mat[2, 0] & inv[4, 2]) ^ (mat[3, 0] & mat[4, 3])\n        return inv % 2\n\n    # For higher dimensions we use Numpy's inverse function\n    # however this function tends to fail and result in a non-symplectic\n    # final matrix if n is too large.\n    max_np_inv = 150\n    if dim <= max_np_inv:\n        return np.linalg.inv(mat).astype(np.int8) % 2\n\n    # For very large matrices  we divide the matrix into 4 blocks of\n    # roughly equal size and use the analytic formula for the inverse\n    # of a block lower-triangular matrix:\n    # inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]\n    # call the inverse function recursively to compute inv(A) and invD\n\n    dim1 = dim // 2\n    mat_a = _inverse_tril(mat[0:dim1, 0:dim1])\n    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim])\n    mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)\n    inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])\n    return inv % 2\n", "path": "qiskit/quantum_info/operators/symplectic/random.py"}], "after_files": [{"content": "# -*- coding: utf-8 -*-\n\n# This code is part of Qiskit.\n#\n# (C) Copyright IBM 2017, 2020\n#\n# This code is licensed under the Apache License, Version 2.0. You may\n# obtain a copy of this license in the LICENSE.txt file in the root directory\n# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\n#\n# Any modifications or derivative works of this code must retain this\n# copyright notice, and modified files need to carry a notice indicating\n# that they have been altered from the originals.\n\"\"\"\nRandom symplectic operator functions\n\"\"\"\n\nimport numpy as np\nfrom numpy.random import default_rng\n\nfrom .clifford import Clifford\nfrom .stabilizer_table import StabilizerTable\nfrom .pauli_table import PauliTable\n\n\ndef random_pauli_table(num_qubits, size=1, seed=None):\n    \"\"\"Return a random PauliTable.\n\n    Args:\n        num_qubits (int): the number of qubits.\n        size (int): Optional. The number of rows of the table (Default: 1).\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        PauliTable: a random PauliTable.\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)\n    return PauliTable(table)\n\n\ndef random_stabilizer_table(num_qubits, size=1, seed=None):\n    \"\"\"Return a random StabilizerTable.\n\n    Args:\n        num_qubits (int): the number of qubits.\n        size (int): Optional. The number of rows of the table (Default: 1).\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        PauliTable: a random StabilizerTable.\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    table = rng.integers(2, size=(size, 2 * num_qubits)).astype(np.bool)\n    phase = rng.integers(2, size=size).astype(np.bool)\n    return StabilizerTable(table, phase)\n\n\ndef random_clifford(num_qubits, seed=None):\n    \"\"\"Return a random Clifford operator.\n\n    The Clifford is sampled using the method of Reference [1].\n\n    Args:\n        num_qubits (int): the number of qubits for the Clifford\n        seed (int or np.random.Generator): Optional. Set a fixed seed or\n                                           generator for RNG.\n\n    Returns:\n        Clifford: a random Clifford operator.\n\n    Reference:\n        1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the\n           structure of the Clifford group*.\n           `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_\n    \"\"\"\n    if seed is None:\n        rng = np.random.default_rng()\n    elif isinstance(seed, np.random.Generator):\n        rng = seed\n    else:\n        rng = default_rng(seed)\n\n    had, perm = _sample_qmallows(num_qubits, rng)\n\n    gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))\n    gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))\n    delta1 = np.eye(num_qubits, dtype=np.int8)\n    delta2 = delta1.copy()\n\n    _fill_tril(gamma1, rng, symmetric=True)\n    _fill_tril(gamma2, rng, symmetric=True)\n    _fill_tril(delta1, rng)\n    _fill_tril(delta2, rng)\n\n    # For large num_qubits numpy.inv function called below can\n    # return invalid output leading to a non-symplectic Clifford\n    # being generated. This can be prevented by manually forcing\n    # block inversion of the matrix.\n    block_inverse_threshold = 50\n\n    # Compute stabilizer table\n    zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)\n    prod1 = np.matmul(gamma1, delta1) % 2\n    prod2 = np.matmul(gamma2, delta2) % 2\n    inv1 = _inverse_tril(delta1, block_inverse_threshold).transpose()\n    inv2 = _inverse_tril(delta2, block_inverse_threshold).transpose()\n    table1 = np.block([[delta1, zero], [prod1, inv1]])\n    table2 = np.block([[delta2, zero], [prod2, inv2]])\n\n    # Apply qubit permutation\n    table = table2[np.concatenate([perm, num_qubits + perm])]\n\n    # Apply layer of Hadamards\n    inds = had * np.arange(1, num_qubits + 1)\n    inds = inds[inds > 0] - 1\n    lhs_inds = np.concatenate([inds, inds + num_qubits])\n    rhs_inds = np.concatenate([inds + num_qubits, inds])\n    table[lhs_inds, :] = table[rhs_inds, :]\n\n    # Apply table\n    table = np.mod(np.matmul(table1, table), 2).astype(np.bool)\n\n    # Generate random phases\n    phase = rng.integers(2, size=2 * num_qubits).astype(np.bool)\n    return Clifford(StabilizerTable(table, phase))\n\n\ndef _sample_qmallows(n, rng=None):\n    \"\"\"Sample from the quantum Mallows distribution\"\"\"\n\n    if rng is None:\n        rng = np.random.default_rng()\n\n    # Hadmard layer\n    had = np.zeros(n, dtype=np.bool)\n\n    # Permutation layer\n    perm = np.zeros(n, dtype=int)\n\n    inds = list(range(n))\n    for i in range(n):\n        m = n - i\n        eps = 4 ** (-m)\n        r = rng.uniform(0, 1)\n        index = -int(np.ceil(np.log2(r + (1 - r) * eps)))\n        had[i] = index < m\n        if index < m:\n            k = index\n        else:\n            k = 2 * m - index - 1\n        perm[i] = inds[k]\n        del inds[k]\n    return had, perm\n\n\ndef _fill_tril(mat, rng, symmetric=False):\n    \"\"\"Add symmetric random ints to off diagonals\"\"\"\n    dim = mat.shape[0]\n    # Optimized for low dimensions\n    if dim == 1:\n        return\n\n    if dim <= 4:\n        mat[1, 0] = rng.integers(2, dtype=np.int8)\n        if symmetric:\n            mat[0, 1] = mat[1, 0]\n        if dim > 2:\n            mat[2, 0] = rng.integers(2, dtype=np.int8)\n            mat[2, 1] = rng.integers(2, dtype=np.int8)\n            if symmetric:\n                mat[0, 2] = mat[2, 0]\n                mat[1, 2] = mat[2, 1]\n        if dim > 3:\n            mat[3, 0] = rng.integers(2, dtype=np.int8)\n            mat[3, 1] = rng.integers(2, dtype=np.int8)\n            mat[3, 2] = rng.integers(2, dtype=np.int8)\n            if symmetric:\n                mat[0, 3] = mat[3, 0]\n                mat[1, 3] = mat[3, 1]\n                mat[2, 3] = mat[3, 2]\n        return\n\n    # Use numpy indices for larger dimensions\n    rows, cols = np.tril_indices(dim, -1)\n    vals = rng.integers(2, size=rows.size, dtype=np.int8)\n    mat[(rows, cols)] = vals\n    if symmetric:\n        mat[(cols, rows)] = vals\n\n\ndef _inverse_tril(mat, block_inverse_threshold):\n    \"\"\"Invert a lower-triangular matrix with unit diagonal.\"\"\"\n    # Optimized inversion function for low dimensions\n    dim = mat.shape[0]\n\n    if dim <= 2:\n        return mat\n\n    if dim <= 5:\n        inv = mat.copy()\n        inv[2, 0] = (mat[2, 0] ^ (mat[1, 0] & mat[2, 1]))\n        if dim > 3:\n            inv[3, 1] = (mat[3, 1] ^ (mat[2, 1] & mat[3, 2]))\n            inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])\n        if dim > 4:\n            inv[4, 2] = ((mat[4, 2] ^ (mat[3, 2] & mat[4, 3]))) & 1\n            inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])\n            inv[4, 0] = mat[4, 0] ^ (mat[1, 0] & inv[4, 1]) ^ (\n                mat[2, 0] & inv[4, 2]) ^ (mat[3, 0] & mat[4, 3])\n        return inv % 2\n\n    # For higher dimensions we use Numpy's inverse function\n    # however this function tends to fail and result in a non-symplectic\n    # final matrix if n is too large.\n    if dim <= block_inverse_threshold:\n        return np.linalg.inv(mat).astype(np.int8) % 2\n\n    # For very large matrices  we divide the matrix into 4 blocks of\n    # roughly equal size and use the analytic formula for the inverse\n    # of a block lower-triangular matrix:\n    # inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]\n    # call the inverse function recursively to compute inv(A) and invD\n\n    dim1 = dim // 2\n    mat_a = _inverse_tril(mat[0:dim1, 0:dim1], block_inverse_threshold)\n    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim], block_inverse_threshold)\n    mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)\n    inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])\n    return inv % 2\n", "path": "qiskit/quantum_info/operators/symplectic/random.py"}]}