Repository: Qiskit/qiskit-terra (python)

Issue:
random_circuit fails if conditional=True for large circuits
<!-- ⚠️ If you do not respect this template, your issue will be closed -->
<!-- ⚠️ Make sure to browse the opened and closed issues -->

### Information

- **Qiskit Terra version**: 0.19.0.dev0+07299e1
- **Python version**:
- **Operating system**:

### What is the current behavior?
When attempting to generate a random circuit using `qiskit.circuit.random.utils.random_circuit` with ~60 qubits with mid-circuit conditionals an error is generated;
```
_bounded_integers.pyx in numpy.random._bounded_integers._rand_int64()

ValueError: low >= high
```


### Steps to reproduce the problem
```python
from qiskit.circuit.random.utils import random_circuit
random_circuit(64, 64, conditional=True)
```
### What is the expected behavior?



### Suggested solutions




Pull request description:
Speed up `random_circuit`
### Summary

There is no need to use the slower `append` (which adds a lot of checking overhead for large circuits) since we fully control the construction of the circuit.

This also overhauls all the randomisation components to produce all the necessary randomness for a given layer in a single go.  This massively reduces the Python-space and Numpy-dispatch overhead of the nested loop.

The changes to the randomisation in this new form mean that on average, more gates will be generated per circuit, because the output will choose 1q gates more frequently than it chooses 2q, which in turn will be more frequent than 3q, but each layer of the "depth" still has to involve every qubit.  The change is because the likelihood of choosing a gate with a given number of qubits is now proportional to how many different gates of that number of qubits there are in the options.

For sample timings, the call `random_circuit(433, 1000, seed=1)` went from 15.2s on my machine (macOS Python 3.10, Intel i7 @ 2.3 GHz) down to 2.3s, so a speed-up of about 6-7x.

<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add a release note file using reno if the change needs to be
  documented in the release notes.
⚠️ If your pull request fixes an open issue, please link to the issue.

- [ ] I have added the tests to cover my changes.
- [ ] I have updated the documentation accordingly.
- [ ] I have read the CONTRIBUTING document.
-->


### Details and comments

If we want to make the circuit generate as many 3q gates as it did before, it's pretty trivial to adjust the weights in the `rng.choice` call to restore the same sort of behaviour as before.

cc: @kdk, @ismaelfaro.  Ismael: I'm tagging you just because it might be of interest, since I know you were using this function and finding it slow before.

Fix #6994.
Close #8425. (obsoleting a stale PR)

Patch:
diff --git a/qiskit/circuit/random/utils.py b/qiskit/circuit/random/utils.py
index ca144cb05..b63f04c74 100644
--- a/qiskit/circuit/random/utils.py
+++ b/qiskit/circuit/random/utils.py
@@ -14,41 +14,14 @@
 
 import numpy as np
 
-from qiskit.circuit import QuantumRegister, ClassicalRegister, QuantumCircuit
+from qiskit.circuit import ClassicalRegister, QuantumCircuit, CircuitInstruction
 from qiskit.circuit import Reset
-from qiskit.circuit.library.standard_gates import (
-    IGate,
-    U1Gate,
-    U2Gate,
-    U3Gate,
-    XGate,
-    YGate,
-    ZGate,
-    HGate,
-    SGate,
-    SdgGate,
-    TGate,
-    TdgGate,
-    RXGate,
-    RYGate,
-    RZGate,
-    CXGate,
-    CYGate,
-    CZGate,
-    CHGate,
-    CRZGate,
-    CU1Gate,
-    CU3Gate,
-    SwapGate,
-    RZZGate,
-    CCXGate,
-    CSwapGate,
-)
+from qiskit.circuit.library import standard_gates
 from qiskit.circuit.exceptions import CircuitError
 
 
 def random_circuit(
-    num_qubits, depth, max_operands=3, measure=False, conditional=False, reset=False, seed=None
+    num_qubits, depth, max_operands=4, measure=False, conditional=False, reset=False, seed=None
 ):
     """Generate random circuit of arbitrary size and form.
 
@@ -65,7 +38,7 @@ def random_circuit(
     Args:
         num_qubits (int): number of quantum wires
         depth (int): layers of operations (i.e. critical path length)
-        max_operands (int): maximum operands of each gate (between 1 and 3)
+        max_operands (int): maximum qubit operands of each gate (between 1 and 4)
         measure (bool): if True, measure all qubits at the end
         conditional (bool): if True, insert middle measurements and conditionals
         reset (bool): if True, insert middle resets
@@ -77,81 +50,157 @@ def random_circuit(
     Raises:
         CircuitError: when invalid options given
     """
-    if max_operands < 1 or max_operands > 3:
-        raise CircuitError("max_operands must be between 1 and 3")
-
-    one_q_ops = [
-        IGate,
-        U1Gate,
-        U2Gate,
-        U3Gate,
-        XGate,
-        YGate,
-        ZGate,
-        HGate,
-        SGate,
-        SdgGate,
-        TGate,
-        TdgGate,
-        RXGate,
-        RYGate,
-        RZGate,
+    if num_qubits == 0:
+        return QuantumCircuit()
+    if max_operands < 1 or max_operands > 4:
+        raise CircuitError("max_operands must be between 1 and 4")
+    max_operands = max_operands if num_qubits > max_operands else num_qubits
+
+    gates_1q = [
+        # (Gate class, number of qubits, number of parameters)
+        (standard_gates.IGate, 1, 0),
+        (standard_gates.SXGate, 1, 0),
+        (standard_gates.XGate, 1, 0),
+        (standard_gates.RZGate, 1, 1),
+        (standard_gates.RGate, 1, 2),
+        (standard_gates.HGate, 1, 0),
+        (standard_gates.PhaseGate, 1, 1),
+        (standard_gates.RXGate, 1, 1),
+        (standard_gates.RYGate, 1, 1),
+        (standard_gates.SGate, 1, 0),
+        (standard_gates.SdgGate, 1, 0),
+        (standard_gates.SXdgGate, 1, 0),
+        (standard_gates.TGate, 1, 0),
+        (standard_gates.TdgGate, 1, 0),
+        (standard_gates.UGate, 1, 3),
+        (standard_gates.U1Gate, 1, 1),
+        (standard_gates.U2Gate, 1, 2),
+        (standard_gates.U3Gate, 1, 3),
+        (standard_gates.YGate, 1, 0),
+        (standard_gates.ZGate, 1, 0),
+    ]
+    if reset:
+        gates_1q.append((Reset, 1, 0))
+    gates_2q = [
+        (standard_gates.CXGate, 2, 0),
+        (standard_gates.DCXGate, 2, 0),
+        (standard_gates.CHGate, 2, 0),
+        (standard_gates.CPhaseGate, 2, 1),
+        (standard_gates.CRXGate, 2, 1),
+        (standard_gates.CRYGate, 2, 1),
+        (standard_gates.CRZGate, 2, 1),
+        (standard_gates.CSXGate, 2, 0),
+        (standard_gates.CUGate, 2, 4),
+        (standard_gates.CU1Gate, 2, 1),
+        (standard_gates.CU3Gate, 2, 3),
+        (standard_gates.CYGate, 2, 0),
+        (standard_gates.CZGate, 2, 0),
+        (standard_gates.RXXGate, 2, 1),
+        (standard_gates.RYYGate, 2, 1),
+        (standard_gates.RZZGate, 2, 1),
+        (standard_gates.RZXGate, 2, 1),
+        (standard_gates.XXMinusYYGate, 2, 2),
+        (standard_gates.XXPlusYYGate, 2, 2),
+        (standard_gates.ECRGate, 2, 0),
+        (standard_gates.CSGate, 2, 0),
+        (standard_gates.CSdgGate, 2, 0),
+        (standard_gates.SwapGate, 2, 0),
+        (standard_gates.iSwapGate, 2, 0),
+    ]
+    gates_3q = [
+        (standard_gates.CCXGate, 3, 0),
+        (standard_gates.CSwapGate, 3, 0),
+        (standard_gates.CCZGate, 3, 0),
+        (standard_gates.RCCXGate, 3, 0),
+    ]
+    gates_4q = [
+        (standard_gates.C3SXGate, 4, 0),
+        (standard_gates.RC3XGate, 4, 0),
     ]
-    one_param = [U1Gate, RXGate, RYGate, RZGate, RZZGate, CU1Gate, CRZGate]
-    two_param = [U2Gate]
-    three_param = [U3Gate, CU3Gate]
-    two_q_ops = [CXGate, CYGate, CZGate, CHGate, CRZGate, CU1Gate, CU3Gate, SwapGate, RZZGate]
-    three_q_ops = [CCXGate, CSwapGate]
 
-    qr = QuantumRegister(num_qubits, "q")
+    gates = gates_1q.copy()
+    if max_operands >= 2:
+        gates.extend(gates_2q)
+    if max_operands >= 3:
+        gates.extend(gates_3q)
+    if max_operands >= 4:
+        gates.extend(gates_4q)
+    gates = np.array(
+        gates, dtype=[("class", object), ("num_qubits", np.int64), ("num_params", np.int64)]
+    )
+    gates_1q = np.array(gates_1q, dtype=gates.dtype)
+
     qc = QuantumCircuit(num_qubits)
 
     if measure or conditional:
         cr = ClassicalRegister(num_qubits, "c")
         qc.add_register(cr)
 
-    if reset:
-        one_q_ops += [Reset]
-
     if seed is None:
         seed = np.random.randint(0, np.iinfo(np.int32).max)
     rng = np.random.default_rng(seed)
 
-    # apply arbitrary random operations at every depth
+    qubits = np.array(qc.qubits, dtype=object, copy=True)
+
+    # Apply arbitrary random operations in layers across all qubits.
     for _ in range(depth):
-        # choose either 1, 2, or 3 qubits for the operation
-        remaining_qubits = list(range(num_qubits))
-        rng.shuffle(remaining_qubits)
-        while remaining_qubits:
-            max_possible_operands = min(len(remaining_qubits), max_operands)
-            num_operands = rng.choice(range(max_possible_operands)) + 1
-            operands = [remaining_qubits.pop() for _ in range(num_operands)]
-            if num_operands == 1:
-                operation = rng.choice(one_q_ops)
-            elif num_operands == 2:
-                operation = rng.choice(two_q_ops)
-            elif num_operands == 3:
-                operation = rng.choice(three_q_ops)
-            if operation in one_param:
-                num_angles = 1
-            elif operation in two_param:
-                num_angles = 2
-            elif operation in three_param:
-                num_angles = 3
-            else:
-                num_angles = 0
-            angles = [rng.uniform(0, 2 * np.pi) for x in range(num_angles)]
-            register_operands = [qr[i] for i in operands]
-            op = operation(*angles)
-
-            # with some low probability, condition on classical bit values
-            if conditional and rng.choice(range(10)) == 0:
-                value = rng.integers(0, np.power(2, num_qubits))
-                op.condition = (cr, value)
-
-            qc.append(op, register_operands)
+        # We generate all the randomness for the layer in one go, to avoid many separate calls to
+        # the randomisation routines, which can be fairly slow.
+
+        # This reliably draws too much randomness, but it's less expensive than looping over more
+        # calls to the rng. After, trim it down by finding the point when we've used all the qubits.
+        gate_specs = rng.choice(gates, size=len(qubits))
+        cumulative_qubits = np.cumsum(gate_specs["num_qubits"], dtype=np.int64)
+        # Efficiently find the point in the list where the total gates would use as many as
+        # possible of, but not more than, the number of qubits in the layer.  If there's slack, fill
+        # it with 1q gates.
+        max_index = np.searchsorted(cumulative_qubits, num_qubits, side="right")
+        gate_specs = gate_specs[:max_index]
+        slack = num_qubits - cumulative_qubits[max_index - 1]
+        if slack:
+            gate_specs = np.hstack((gate_specs, rng.choice(gates_1q, size=slack)))
+
+        # For efficiency in the Python loop, this uses Numpy vectorisation to pre-calculate the
+        # indices into the lists of qubits and parameters for every gate, and then suitably
+        # randomises those lists.
+        q_indices = np.empty(len(gate_specs) + 1, dtype=np.int64)
+        p_indices = np.empty(len(gate_specs) + 1, dtype=np.int64)
+        q_indices[0] = p_indices[0] = 0
+        np.cumsum(gate_specs["num_qubits"], out=q_indices[1:])
+        np.cumsum(gate_specs["num_params"], out=p_indices[1:])
+        parameters = rng.uniform(0, 2 * np.pi, size=p_indices[-1])
+        rng.shuffle(qubits)
+
+        # We've now generated everything we're going to need.  Now just to add everything.  The
+        # conditional check is outside the two loops to make the more common case of no conditionals
+        # faster, since in Python we don't have a compiler to do this for us.
+        if conditional:
+            is_conditional = rng.random(size=len(gate_specs)) < 0.1
+            condition_values = rng.integers(
+                0, 1 << min(num_qubits, 63), size=np.count_nonzero(is_conditional)
+            )
+            c_ptr = 0
+            for gate, q_start, q_end, p_start, p_end, is_cond in zip(
+                gate_specs["class"],
+                q_indices[:-1],
+                q_indices[1:],
+                p_indices[:-1],
+                p_indices[1:],
+                is_conditional,
+            ):
+                operation = gate(*parameters[p_start:p_end])
+                if is_cond:
+                    operation.condition = (cr, condition_values[c_ptr])
+                    c_ptr += 1
+                qc._append(CircuitInstruction(operation=operation, qubits=qubits[q_start:q_end]))
+        else:
+            for gate, q_start, q_end, p_start, p_end in zip(
+                gate_specs["class"], q_indices[:-1], q_indices[1:], p_indices[:-1], p_indices[1:]
+            ):
+                operation = gate(*parameters[p_start:p_end])
+                qc._append(CircuitInstruction(operation=operation, qubits=qubits[q_start:q_end]))
 
     if measure:
-        qc.measure(qr, cr)
+        qc.measure(qc.qubits, cr)
 
     return qc
diff --git a/releasenotes/notes/speedup-random-circuits-8d3b724cce1faaad.yaml b/releasenotes/notes/speedup-random-circuits-8d3b724cce1faaad.yaml
new file mode 100644
index 000000000..d1e7696b9
--- /dev/null
+++ b/releasenotes/notes/speedup-random-circuits-8d3b724cce1faaad.yaml
@@ -0,0 +1,14 @@
+---
+features:
+  - |
+    Random-circuit generation with ``qiskit.circuit.random.random_circuit`` is
+    now significantly faster for large circuits.
+  - Random-circuit generation with ``qiskit.circuit.random.random_circuit`` will
+    now output all "standard" gates in Qiskit's circuit library (:mod:`.circuit.library`).
+    This includes two 4-qubit gates :class:`.C3SXGate` and :class:`.RC3XGate`, and the
+    allowed values of ``max_operands`` have been expanded accordingly.
+upgrade:
+  - |
+    The exact circuit returned by ``qiskit.circuit.random.random_circuit`` for a
+    given seed has changed.  This is due to efficiency improvements in the
+    internal random-number generation for the function.


Test patch:
diff --git a/test/python/circuit/test_random_circuit.py b/test/python/circuit/test_random_circuit.py
index 93ac4f6de..db5acb851 100644
--- a/test/python/circuit/test_random_circuit.py
+++ b/test/python/circuit/test_random_circuit.py
@@ -52,3 +52,14 @@ class TestCircuitRandom(QiskitTestCase):
         circ = random_circuit(num_qubits, depth, conditional=True, reset=True, seed=5)
         self.assertEqual(circ.width(), 2 * num_qubits)
         self.assertIn("reset", circ.count_ops())
+
+    def test_large_conditional(self):
+        """Test that conditions do not fail with large conditionals.  Regression test of gh-6994."""
+        # The main test is that this call actually returns without raising an exception.
+        circ = random_circuit(64, 2, conditional=True, seed=0)
+        # Test that at least one instruction had a condition generated.  It's possible that this
+        # fails due to very bad luck with the random seed - if so, change the seed to ensure that a
+        # condition _is_ generated, because we need to test that generation doesn't error.
+        conditions = (getattr(instruction.operation, "condition", None) for instruction in circ)
+        conditions = [x for x in conditions if x is not None]
+        self.assertNotEqual(conditions, [])