qudit_expr/library/
gates.rs

1use crate::UnitaryExpression;
2use crate::UnitarySystemExpression;
3use qudit_core::QuditSystem;
4use qudit_core::Radices;
5
6/// The identity or no-op gate.
7#[cfg_attr(feature = "python", pyo3::pyfunction)]
8#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
9pub fn IGate(radix: usize) -> UnitaryExpression {
10    let proto = format!("I<{}>()", radix);
11    let mut body = "".to_string();
12    body += "[";
13    for i in 0..radix {
14        body += "[";
15        for j in 0..radix {
16            if i == j {
17                body += "1,";
18            } else {
19                body += "0,";
20            }
21        }
22        body += "],";
23    }
24    body += "]";
25
26    UnitaryExpression::new(proto + "{" + &body + "}")
27}
28
29/// The one-qudit Hadamard gate. This is a Clifford/Weyl-Heisenberg gate.
30///
31/// The qubit (radix = 2) Hadamard gate is given by the following matrix:
32///
33/// $$
34/// \begin{pmatrix}
35///     \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\\\
36///     \frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} \\\\
37/// \end{pmatrix}
38/// $$
39///
40/// However, generally it is given by the following formula:
41///
42/// $$
43/// H = \frac{1}{\sqrt{d}} \sum_{ij} \omega^{ij} \ket{i}\bra{j}
44/// $$
45///
46/// where
47///
48/// $$
49/// \omega = \exp\Big(\frac{2\pi i}{d}\Big)
50/// $$
51///
52/// and $d$ is the number of levels (2 levels is a qubit, 3 levels is a qutrit,
53/// etc.)
54///
55/// References:
56/// - <https://www.frontiersin.org/articles/10.3389/fphy.2020.589504/full>
57/// - <https://pubs.aip.org/aip/jmp/article-abstract/56/3/032202/763827>
58/// - <https://arxiv.org/pdf/1701.07902.pdf>
59#[cfg_attr(feature = "python", pyo3::pyfunction)]
60#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
61pub fn HGate(radix: usize) -> UnitaryExpression {
62    let proto = format!("H<{}>()", radix);
63    let mut body = "".to_string();
64    if radix == 2 {
65        body += "[[1/sqrt(2), 1/sqrt(2)], [1/sqrt(2), ~1/sqrt(2)]]";
66        return UnitaryExpression::new(proto + "{" + &body + "}");
67    }
68    let omega = format!("e^(2*π*i/{})", radix);
69    let invsqrt = format!("1/sqrt({})", radix);
70    body += invsqrt.as_str();
71    body += " * ";
72    body += "[";
73    for i in 0..radix {
74        body += "[";
75        for j in 0..radix {
76            body += &format!("{}^({}*{}), ", omega, i, j);
77        }
78        body += "],";
79    }
80    body += "]";
81
82    UnitaryExpression::new(proto + "{" + &body + "}")
83}
84
85/// The qudit swap gate. This is a two-qudit Clifford/Weyl-Heisenberg gate
86/// that swaps the state of two qudits.
87///
88/// The qubit (radix = 2) version is given by the following matrix:
89///
90/// $$
91/// \begin{pmatrix}
92///     1 & 0 & 0 & 0 \\\\
93///     0 & 0 & 1 & 0 \\\\
94///     0 & 1 & 0 & 0 \\\\
95///     0 & 0 & 0 & 1 \\\\
96/// \end{pmatrix}
97/// $$
98///
99/// The qutrit (radix = 3) version is given by the following matrix:
100///
101/// $$
102/// \begin{pmatrix}
103///     1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
104///     0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\
105///     0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\\
106///     0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
107///     0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
108///     0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\\\
109///     0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
110///     0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\\\
111///     0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\\
112/// \end{pmatrix}
113/// $$
114///
115/// However, generally it is given by the following formula:
116///
117/// $$
118/// SWAP_d = \sum_{a, b} \ket{ab}\bra{ba}
119/// $$
120///
121/// where $d$ is the number of levels (2 levels is a qubit, 3 levels is a
122/// qutrit, etc.)
123///
124/// References:
125/// - <https://link.springer.com/article/10.1007/s11128-013-0621-x>
126/// - <https://arxiv.org/pdf/1105.5485.pdf>
127#[cfg_attr(feature = "python", pyo3::pyfunction)]
128#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
129pub fn SwapGate(radix: usize) -> UnitaryExpression {
130    let proto = format!("Swap<{}, {}>()", radix, radix);
131    let mut body = "".to_string();
132    body += "[";
133    for i in 0..radix {
134        body += "[";
135        let a_i = i / radix;
136        let b_i = i % radix;
137        for j in 0..radix {
138            let a_j = j / radix;
139            let b_j = j % radix;
140            if a_i == b_j && b_i == a_j {
141                body += "1,";
142            } else {
143                body += "0,";
144            }
145        }
146        body += "],";
147    }
148    body += "]";
149
150    UnitaryExpression::new(proto + "{" + &body + "}")
151}
152
153/// The one-qudit shift (X) gate. This is a Weyl-Heisenberg gate.
154///
155/// This gate shifts the state of a qudit up by one level modulo. For
156/// example, the shift gate on a qubit is the Pauli-X gate. The shift
157/// gate on a qutrit is the following matrix:
158///
159/// $$
160/// \begin{pmatrix}
161///     0 & 0 & 1 \\\\
162///     1 & 0 & 0 \\\\
163///     0 & 1 & 0 \\\\
164/// \end{pmatrix}
165/// $$
166///
167/// The shift gate is generally given by the following formula:
168///
169/// $$
170/// \begin{equation}
171///     X = \sum_a |a + 1 mod d ><a|
172/// \end{equation}
173/// $$
174///
175/// where d is the number of levels (2 levels is a qubit, 3 levels is
176/// a qutrit, etc.)
177///
178/// References:
179///     - <https://arxiv.org/pdf/2302.07966.pdf>
180#[cfg_attr(feature = "python", pyo3::pyfunction)]
181#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
182pub fn XGate(radix: usize) -> UnitaryExpression {
183    let proto = format!("X<{}>()", radix);
184    let mut body = "[".to_string();
185    for i in 0..radix {
186        body += "[";
187        for j in 0..radix {
188            if (j + 1) % radix == i {
189                body += "1, ";
190            } else {
191                body += "0, ";
192            }
193        }
194        body += "],";
195    }
196    body += "]";
197
198    UnitaryExpression::new(proto + "{" + &body + "}")
199}
200
201/// The one-qudit clock (Z) gate. This is a Weyl-Heisenberg gate.
202///
203/// This gate shifts the state of a qudit up by one level modulo. For
204/// example, the clock gate on a qubit is the Pauli-Z gate. The clock
205/// gate on a qutrit is the following matrix:
206///
207/// $$
208/// \begin{pmatrix}
209///     1 & 0 & 0 \\\\
210///     0 & e^{\frac{2\pi i}{3} & 0 \\\\
211///     0 & 0 & e^{\frac{4\pi i}{3} \\\\
212/// \end{pmatrix}
213/// $$
214///
215/// The clock gate is generally given by the following formula:
216///
217/// $$
218/// \begin{equation}
219///     X = \sum_a e^{\frac{2a\pi i}{d}|a><a|
220/// \end{equation}
221/// $$
222///
223/// where d is the number of levels (2 levels is a qubit, 3 levels is
224/// a qutrit, etc.)
225///
226/// References:
227///     - <https://arxiv.org/pdf/2302.07966.pdf>
228#[cfg_attr(feature = "python", pyo3::pyfunction)]
229#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
230pub fn ZGate(radix: usize) -> UnitaryExpression {
231    let proto = format!("Z<{}>()", radix);
232    let mut body = "[".to_string();
233    for i in 0..radix {
234        body += "[";
235        for j in 0..radix {
236            if i == j {
237                if i == 0 {
238                    body += "1, ";
239                } else {
240                    body += &format!("e^((2*{}*π*i)/{}), ", j, radix);
241                }
242            } else {
243                body += "0, ";
244            }
245        }
246        body += "],";
247    }
248    body += "]";
249
250    UnitaryExpression::new(proto + "{" + &body + "}")
251}
252
253/// The single-qudit phase gate.
254///
255/// The common qubit phase gate is given by the following matrix:
256///
257/// $$
258/// \begin{pmatrix}
259///     1 & 0 \\\\
260///     0 & \exp({i\theta}) \\\\
261/// \end{pmatrix}
262/// $$
263///
264/// The qutrit phase gate has two parameterized relative phases:
265///
266/// $$
267/// \begin{pmatrix}
268///     1 & 0 & 0 \\\\
269///     0 & \exp({i\theta_0}) & 0 \\\\
270///    0 & 0 & \exp({i\theta_1}) \\\\
271/// \end{pmatrix}
272/// $$
273///
274/// The d-level phase gate has d-1 parameterized relative phases. This
275/// gate is Clifford iff all of the relative phases are powers of roots
276/// of unity.
277///
278/// References:
279/// - <https://www.nature.com/articles/s41467-022-34851-z>
280/// - <https://arxiv.org/pdf/2204.13681.pdf>
281#[cfg_attr(feature = "python", pyo3::pyfunction)]
282#[cfg_attr(feature = "python", pyo3(signature = (radix = 2)))]
283pub fn PGate(radix: usize) -> UnitaryExpression {
284    let mut proto = format!("P<{}>(", radix);
285    for i in 0..radix - 1 {
286        proto += "θ";
287        proto += &i.to_string();
288        proto += ", ";
289    }
290    proto += ")";
291
292    let mut body = "".to_string();
293    body += "[";
294    for i in 0..radix {
295        body += "[";
296        for j in 0..radix {
297            if i == j {
298                if i == 0 {
299                    body += "1, ";
300                } else {
301                    body += &format!("e^(i*θ{}), ", i - 1);
302                }
303            } else {
304                body += "0, ";
305            }
306        }
307        body += "],";
308    }
309    body += "]";
310
311    UnitaryExpression::new(proto + "{" + &body + "}")
312}
313
314#[cfg_attr(feature = "python", pyo3::pyfunction)]
315pub fn U3Gate() -> UnitaryExpression {
316    let proto = "U(θ0, θ1, θ2)";
317    let body = "[
318            [cos(θ0/2), ~e^(i*θ2)*sin(θ0/2)],
319            [e^(i*θ1)*sin(θ0/2), e^(i*(θ1+θ2))*cos(θ0/2)]
320    ]";
321    UnitaryExpression::new(proto.to_owned() + "{" + body + "}")
322}
323
324fn generate_embedded_two_param_su2(dimension: usize, i: usize, j: usize) -> UnitaryExpression {
325    let proto = format!("U<{dimension}>(θ0, θ1)");
326    let mut body = String::from("[");
327    for row in 0..dimension {
328        body += "[";
329        for col in 0..dimension {
330            if row == i && col == i {
331                body += "e^(i*θ1)*cos(θ0/2)";
332            } else if row == i && col == j {
333                body += "~sin(θ0/2)";
334            } else if row == j && col == i {
335                body += "sin(θ0/2)";
336            } else if row == j && col == j {
337                body += "e^(i*(~θ1))*cos(θ0/2)";
338            } else if row == col {
339                body += "1";
340            } else {
341                body += "0";
342            }
343            body += ",";
344        }
345        body += "]";
346    }
347    body += "]";
348    UnitaryExpression::new(proto.to_owned() + "{" + &body + "}")
349}
350
351fn generate_embedded_su2(dimension: usize, i: usize, j: usize) -> UnitaryExpression {
352    let proto = format!("U<{dimension}>(θ0, θ1, θ2)");
353    let mut body = String::from("[");
354    for row in 0..dimension {
355        body += "[";
356        for col in 0..dimension {
357            if row == i && col == i {
358                body += "cos(θ0/2)";
359            } else if row == i && col == j {
360                body += "~e^(i*θ2)*sin(θ0/2)";
361            } else if row == j && col == i {
362                body += "e^(i*θ1)*sin(θ0/2)";
363            } else if row == j && col == j {
364                body += "e^(i*(θ1+θ2))*cos(θ0/2)";
365            } else if row == col {
366                body += "1";
367            } else {
368                body += "0";
369            }
370            body += ",";
371        }
372        body += "]";
373    }
374    body += "]";
375    UnitaryExpression::new(proto.to_owned() + "{" + &body + "}")
376}
377
378fn embed_one_larger(unitary: UnitaryExpression) -> UnitaryExpression {
379    let dimension = unitary.dimension() + 1;
380    let mut one_larger = UnitaryExpression::identity(unitary.name(), [dimension]);
381    one_larger.embed(unitary, 1, 1);
382    one_larger
383}
384
385/// Generates a fully parameterized unitary expression
386///
387/// References:
388/// - de Guise, Hubert, Olivia Di Matteo, and Luis L. Sánchez-Soto.
389///   "Simple factorization of unitary transformations."
390///   Physical Review A 97.2 (2018): 022328.
391#[cfg_attr(feature = "python", pyo3::pyfunction)]
392#[cfg_attr(feature = "python", pyo3(signature = (radices = Radices::new([2]))))]
393pub fn ParameterizedUnitary(radices: Radices) -> UnitaryExpression {
394    let dimension = radices.dimension();
395
396    if dimension < 2 {
397        panic!("Cannot generate parameterized unitary for dimension 1 or less");
398    }
399
400    if dimension == 2 {
401        return U3Gate();
402    }
403
404    let right = {
405        let one_smaller = ParameterizedUnitary(Radices::new([dimension - 1]));
406        embed_one_larger(one_smaller)
407    };
408
409    let left = {
410        let mut acm = generate_embedded_su2(dimension, dimension - 2, dimension - 1);
411        for i in (0..=(dimension - 3)).rev() {
412            let j = i + 1;
413            let two = generate_embedded_two_param_su2(dimension, i, j);
414            acm = acm.dot(two)
415        }
416        acm
417    };
418
419    let mut expression = left.dot(right);
420    expression.set_radices(radices);
421    expression
422}
423
424/// Invert an expression
425#[cfg_attr(feature = "python", pyo3::pyfunction)]
426pub fn Invert(mut expr: UnitaryExpression) -> UnitaryExpression {
427    expr.dagger();
428    expr
429}
430
431/// Calculates the cartesian product of the control levels.
432///
433/// # Examples
434///
435/// ```ignore
436/// let control_levels = vec![vec![0, 1], vec![0, 1]];
437/// let prod = cartesian_product(control_levels);
438/// assert_eq!(prod, vec![
439///    vec![0, 0],
440///    vec![1, 0],
441///    vec![0, 1],
442///    vec![1, 1],
443/// ]);
444/// ```
445fn cartesian_product(control_levels: Vec<Vec<usize>>) -> Vec<Vec<usize>> {
446    let mut prod = vec![];
447    for level in control_levels.into_iter() {
448        if prod.is_empty() {
449            for l in level.into_iter() {
450                prod.push(vec![l]);
451            }
452        } else {
453            let mut new_prod = vec![];
454            for l in level.into_iter() {
455                for v in prod.iter() {
456                    let mut v_new = v.clone();
457                    v_new.push(l);
458                    new_prod.push(v_new);
459                }
460            }
461            prod = new_prod;
462        }
463    }
464    prod
465}
466
467/// An arbitrary controlled gate.
468///
469/// Given any gate, ControlledGate can add control qudits.
470///
471/// A controlled gate adds arbitrarily controls, and is generalized
472/// for qudit or even mixed-qudit representation.
473///
474/// A controlled gate has a circuit structure as follows:
475///
476/// ```text
477///     controls ----/----■----
478///                       |
479///                      .-.
480///     targets  ----/---|G|---
481///                      '-'
482/// ```
483///
484/// Where $G$ is the gate being controlled.
485///
486/// To calculate the unitary for a controlled gate, given the unitary of
487/// the gate being controlled, we can use the following equation:
488///
489/// $$U_{control} = P_i \otimes I + P_c \otimes G$$
490///
491/// Where $P_i$ is the projection matrix for the states that don't
492/// activate the gate, $P_c$ is the projection matrix for the
493/// states that do activate the gate, $I$ is the identity matrix
494/// of dimension equal to the gate being controlled, and $G$ is
495/// the unitary matrix of the gate being controlled.
496///
497/// In the simple case of a normal qubit CNOT ($G = X$), $P_i$ and $P_c$
498/// are defined as follows:
499///
500/// $$
501///     P_i = \ket{0}\bra{0}
502///     P_c = \ket{1}\bra{1}
503/// $$
504///
505/// This is because the $\ket{0}$ state is the state that doesn't
506/// activate the gate, and the $\ket{1}$ state is the state that
507/// does activate the gate.
508///
509/// We can also decide to invert this, and have the $\ket{0}$
510/// state activate the gate, and the $\ket{1}$ state not activate
511/// the gate. This is equivalent to swapping $P_i$ and $P_c$,
512/// and usually drawn diagrammatically as follows:
513///
514/// ```text
515///     controls ----/----□----
516///                       |
517///                      .-.
518///     targets  ----/---|G|---
519///                      '-'
520/// ```
521///
522/// When we add more controls the projection matrices become more complex,
523/// but the basic idea stays the same: we have a projection matrix for
524/// the states that activate the gate, and a projection matrix for the
525/// states that don't activate the gate. As in the case of a toffoli gate,
526/// the projection matrices are defined as follows:
527///
528/// $$
529///     P_i = \ket{00}\bra{00} + \ket{01}\bra{01} + \ket{10}\bra{10}
530///     P_c = \ket{11}\bra{11}
531/// $$
532///
533/// This is because the $\ket{00}$, $\ket{01}$, and
534/// $\ket{10}$ states are the states that don't activate the
535/// gate, and the $\ket{11}$ state is the state that does
536/// activate the gate.
537///
538/// With qudits, we have more states and as such, more complex
539/// projection matrices; however, the basic idea is the same.
540/// For example, a qutrit controlled-not gate that is activated by
541/// the $\ket{2}$ state and not activated by the $\ket{0}$
542/// and $\ket{1}$ states is defined as follows:
543///
544/// $$
545///     P_i = \ket{0}\bra{0} + \ket{1}\bra{1}
546///     P_c = \ket{2}\bra{2}
547/// $$
548///
549/// One interesting concept with qudits is that we can have multiple
550/// active control levels. For example, a qutrit controlled-not gate that
551/// is activated by the $\ket{1}$ and $\ket{2}$ states
552/// and not activated by the $\ket{0}$ state is defined similarly
553/// as follows:
554///
555/// $$
556///     P_i = \ket{0}\bra{0}
557///     P_c = \ket{1}\bra{1} + \ket{2}\bra{2}
558/// $$
559///
560/// Note that we can always define $P_i$ simply from $P_c$:
561///
562/// $$P_i = I_p - P_c$$
563///
564/// Where $I_p$ is the identity matrix of dimension equal to the
565/// dimension of the control qudits. This leaves us with out final
566/// equation:
567///
568///
569/// $$U_{control} = (I_p - P_c) \otimes I + P_c \otimes G$$
570///
571/// If, G is a unitary-valued function of real parameters, then the
572/// gradient of the controlled gate simply discards the constant half
573/// of the equation:
574///
575/// $$
576///     \frac{\partial U_{control}}{\partial \theta} =
577///         P_c \otimes \frac{\partial G}{\partial \theta}
578/// $$
579///
580/// # Arguments
581///
582/// * `expr` - The gate to control.
583///
584/// * `control_radixes` - The number of levels for each control qudit.
585///
586/// * `control_levels` - The levels of the control qudits that activate the
587///   gate. If more than one level is selected, the subspace spanned by the
588///   levels acts as a control subspace. If all levels are selected for a
589///   given qudit, the operation is equivalent to the original gate without
590///   controls.
591/// # Panics
592///
593/// * If `control_radixes` and `control_levels` have different lengths.
594///
595/// * If `control_levels` contains an empty level.
596///
597/// * If any level in `control_levels` is greater than or equal to the
598///   corresponding radix in `control_radixes`.
599///
600/// * If any level in `control_levels` is not unique.
601#[cfg_attr(feature = "python", pyo3::pyfunction)]
602#[cfg_attr(feature = "python", pyo3(signature = (expr, control_radices = Radices::new([2]), control_levels = None)))]
603pub fn Controlled(
604    expr: UnitaryExpression,
605    control_radices: Radices,
606    control_levels: Option<Vec<Vec<usize>>>,
607) -> UnitaryExpression {
608    let control_levels = match control_levels {
609        Some(levels) => levels,
610        None => {
611            // Generate default control_levels: each control qudit activated by its highest level
612            control_radices
613                .iter()
614                .map(|&radix| vec![(usize::from(radix) - 1)])
615                .collect()
616        }
617    };
618
619    if control_radices.len() != control_levels.len() {
620        panic!("control_radices and control_levels must have the same length");
621    }
622
623    if control_levels.iter().any(|levels| levels.is_empty()) {
624        panic!("control_levels must not contain empty levels");
625    }
626
627    if control_levels
628        .iter()
629        .map(|levels| levels.iter().copied())
630        .zip(control_radices.iter())
631        .any(|(mut levels, radix)| levels.any(|level| level >= usize::from(*radix)))
632    {
633        panic!("Expected control levels to be less than the number of levels.");
634    }
635
636    // Check that all levels in control_levels are unique
637    let mut control_level_sets = control_levels.clone();
638    for level in control_level_sets.iter_mut() {
639        level.sort();
640        level.dedup();
641    }
642    if control_level_sets
643        .iter()
644        .zip(control_levels.iter())
645        .any(|(level_dedup, level)| level.len() != level_dedup.len())
646    {
647        panic!("Expected control levels to be unique.");
648    }
649
650    let gate_expr = expr;
651    let gate_dim = gate_expr.dimension();
652
653    // Build appropriately sized identity expression
654    let name = format!("Controlled({})", gate_expr.name());
655    let radices = control_radices.concat(&gate_expr.radices());
656    let mut expr = UnitaryExpression::identity(&name, radices);
657
658    // Embed gate expression into identity expression at correct spots
659    let diagonal_indices: Vec<usize> = cartesian_product(control_levels)
660        .into_iter()
661        .map(|block_idx_expansion| control_radices.compress(&block_idx_expansion))
662        .map(|block_diag_idx| block_diag_idx * gate_dim)
663        .collect();
664
665    for diagonal_idx in diagonal_indices.iter() {
666        expr.embed(gate_expr.clone(), *diagonal_idx, *diagonal_idx);
667    }
668
669    expr
670}
671
672#[cfg_attr(feature = "python", pyo3::pyfunction)]
673#[cfg_attr(feature = "python", pyo3(signature = (expr, control_radices = Radices::new([2]), control_levels = None)))]
674pub fn ClassicallyControlled(
675    expr: UnitaryExpression,
676    control_radices: Radices,
677    control_levels: Option<Vec<Vec<usize>>>,
678) -> UnitarySystemExpression {
679    let control_levels = match control_levels {
680        Some(levels) => levels,
681        None => {
682            // Generate default control_levels: each control qudit activated by its highest level
683            control_radices
684                .iter()
685                .map(|&radix| vec![(usize::from(radix) - 1)])
686                .collect()
687        }
688    };
689
690    if control_radices.len() != control_levels.len() {
691        panic!("control_radices and control_levels must have the same length");
692    }
693
694    if control_levels.iter().any(|levels| levels.is_empty()) {
695        panic!("control_levels must not contain empty levels");
696    }
697
698    if control_levels
699        .iter()
700        .map(|levels| levels.iter().copied())
701        .zip(control_radices.iter())
702        .any(|(mut levels, radix)| levels.any(|level| level >= usize::from(*radix)))
703    {
704        panic!("Expected control levels to be less than the number of levels.");
705    }
706
707    // Check that all levels in control_levels are unique
708    let mut control_level_sets = control_levels.clone();
709    for level in control_level_sets.iter_mut() {
710        level.sort();
711        level.dedup();
712    }
713    if control_level_sets
714        .iter()
715        .zip(control_levels.iter())
716        .any(|(level_dedup, level)| level.len() != level_dedup.len())
717    {
718        panic!("Expected control levels to be unique.");
719    }
720
721    let diagonal_indices: Vec<usize> = cartesian_product(control_levels)
722        .into_iter()
723        .map(|block_idx_expansion| control_radices.compress(&block_idx_expansion))
724        .collect();
725
726    let new_radices = control_radices
727        .iter()
728        .map(|&r| r.into())
729        .collect::<Vec<_>>();
730    expr.classically_control(&diagonal_indices, &new_radices)
731}
732
733#[cfg(feature = "python")]
734mod python {
735    use crate::python::PyExpressionRegistrar;
736    use pyo3::prelude::*;
737
738    /// Registers the gate library with the Python module.
739    fn register(parent_module: &Bound<'_, PyModule>) -> PyResult<()> {
740        parent_module.add_function(wrap_pyfunction!(super::IGate, parent_module)?)?;
741        parent_module.add_function(wrap_pyfunction!(super::HGate, parent_module)?)?;
742        parent_module.add_function(wrap_pyfunction!(super::SwapGate, parent_module)?)?;
743        parent_module.add_function(wrap_pyfunction!(super::XGate, parent_module)?)?;
744        parent_module.add_function(wrap_pyfunction!(super::ZGate, parent_module)?)?;
745        parent_module.add_function(wrap_pyfunction!(super::PGate, parent_module)?)?;
746        parent_module.add_function(wrap_pyfunction!(super::U3Gate, parent_module)?)?;
747        parent_module.add_function(wrap_pyfunction!(
748            super::ParameterizedUnitary,
749            parent_module
750        )?)?;
751        parent_module.add_function(wrap_pyfunction!(super::Invert, parent_module)?)?;
752        parent_module.add_function(wrap_pyfunction!(super::Controlled, parent_module)?)?;
753        parent_module.add_function(wrap_pyfunction!(
754            super::ClassicallyControlled,
755            parent_module
756        )?)?;
757        Ok(())
758    }
759    inventory::submit!(PyExpressionRegistrar { func: register });
760}