Skip to content

Organize two_qubit_decompose code #15880

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ShellyGarion merged 4 commits into from
Apr 1, 2026
Merged

Organize two_qubit_decompose code #15880

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1e50dca
organize two_qubit_decompose code
ShellyGarion Mar 26, 2026
fb86c61
move back consts and functions from common. add comments
ShellyGarion Mar 29, 2026
0569b40
minor change order
ShellyGarion Mar 29, 2026
f4eec58
merge with main branch
ShellyGarion Apr 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 20 additions & 56 deletions crates/synthesis/src/two_qubit_decompose/basis_decomposer.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ use approx::{abs_diff_eq, relative_eq};
use num_complex::{Complex64, ComplexFloat};
use num_traits::Zero;
use smallvec::{SmallVec, smallvec};
use std::f64::consts::{FRAC_1_SQRT_2, PI};
use std::f64::consts::{FRAC_1_SQRT_2, FRAC_PI_2, FRAC_PI_4, PI, TAU}; // TAU=2*PI
Comment thread
ShellyGarion marked this conversation as resolved.
Outdated

use ndarray::linalg::kron;
use ndarray::prelude::*;
Expand All @@ -24,7 +24,10 @@ use numpy::PyReadonlyArray2;
use pyo3::exceptions::PyValueError;
use pyo3::prelude::*;

use super::common::{DEFAULT_FIDELITY, TraceToFidelity, rx_matrix, rz_matrix, transpose_conjugate};
use super::common::{
DEFAULT_FIDELITY, IPZ, K12L_ARR, K12R_ARR, TraceToFidelity, real_trace_transform, rx_matrix,
rz_matrix, transpose_conjugate, u4_to_su4,
};
use super::gate_sequence::{TwoQubitGateSequence, TwoQubitSequenceVec};
use super::weyl_decomposition::{__num_basis_gates, _num_basis_gates, TwoQubitWeylDecomposition};

Expand All @@ -44,12 +47,9 @@ use qiskit_circuit::gate_matrix::{CX_GATE, H_GATE, ONE_QUBIT_IDENTITY};
use qiskit_circuit::instruction::{Instruction, Parameters};
use qiskit_circuit::operations::{Operation, OperationRef, Param, StandardGate};
use qiskit_circuit::packed_instruction::PackedOperation;
use qiskit_circuit::util::{C_M_ONE, C_ONE, C_ZERO, GateArray1Q, IM, M_IM, c64};
use qiskit_circuit::util::{C_M_ONE, C_ONE, GateArray1Q, IM, M_IM, c64};
use qiskit_circuit::{NoBlocks, Qubit};

const PI2: f64 = PI / 2.;
const PI4: f64 = PI / 4.;
const TWO_PI: f64 = 2.0 * PI;
// Worst case length is 5x 1q gates for each 1q decomposition + 1x 2q gate
// We might overallocate a bit if the euler basis is different but
// the worst case is just 16 extra elements with just a String and 2 smallvecs
Expand All @@ -58,8 +58,6 @@ const TWO_PI: f64 = 2.0 * PI;
// Python space.
const TWO_QUBIT_SEQUENCE_DEFAULT_CAPACITY: usize = 21;

static IPZ: GateArray1Q = [[IM, C_ZERO], [C_ZERO, M_IM]];

#[derive(Clone, Debug)]
#[allow(non_snake_case)]
#[pyclass(
Expand Down Expand Up @@ -192,7 +190,8 @@ impl TwoQubitBasisDecomposer {
})
.collect();
let mut euler_matrix_q0 = rx_matrix(euler_q0[0][1]).dot(&rz_matrix(euler_q0[0][0]));
euler_matrix_q0 = rz_matrix(euler_q0[0][2] + euler_q0[1][0] + PI2).dot(&euler_matrix_q0);
euler_matrix_q0 =
rz_matrix(euler_q0[0][2] + euler_q0[1][0] + FRAC_PI_2).dot(&euler_matrix_q0);
self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0);
let mut euler_matrix_q1 = rz_matrix(euler_q1[0][1]).dot(&rx_matrix(euler_q1[0][0]));
euler_matrix_q1 = rx_matrix(euler_q1[0][2] + euler_q1[1][0]).dot(&euler_matrix_q1);
Expand All @@ -219,10 +218,10 @@ impl TwoQubitBasisDecomposer {
smallvec![euler_q1[1][1]],
smallvec![1],
));
global_phase += PI2;
global_phase += FRAC_PI_2;
gates.push((StandardGate::CX.into(), smallvec![], smallvec![0, 1]));
let mut euler_matrix_q0 =
rx_matrix(euler_q0[2][1]).dot(&rz_matrix(euler_q0[1][2] + euler_q0[2][0] + PI2));
rx_matrix(euler_q0[2][1]).dot(&rz_matrix(euler_q0[1][2] + euler_q0[2][0] + FRAC_PI_2));
euler_matrix_q0 = rz_matrix(euler_q0[2][2]).dot(&euler_matrix_q0);
self.append_1q_sequence(&mut gates, &mut global_phase, euler_matrix_q0.view(), 0);
let mut euler_matrix_q1 =
Expand Down Expand Up @@ -252,7 +251,7 @@ impl TwoQubitBasisDecomposer {
let mut gates = Vec::new();
let mut global_phase = target_decomposed.global_phase;
global_phase -= 3. * self.basis_decomposer.global_phase;
global_phase = global_phase.rem_euclid(TWO_PI);
global_phase = global_phase.rem_euclid(TAU);
let atol = 1e-10; // absolute tolerance for floats
// Decompose source unitaries to zxz
let euler_q0: Vec<[f64; 3]> = decomposition
Expand Down Expand Up @@ -286,7 +285,7 @@ impl TwoQubitBasisDecomposer {
x12_phase = PI * x12.cos();
}
let x02_add = x12 - euler_q0[1][0];
let x12_is_half_pi = abs_diff_eq!(x12, PI2, epsilon = atol);
let x12_is_half_pi = abs_diff_eq!(x12, FRAC_PI_2, epsilon = atol);

let mut euler_matrix_q0 = rx_matrix(euler_q0[0][1]).dot(&rz_matrix(euler_q0[0][0]));
if x12_is_non_zero && x12_is_pi_mult {
Expand Down Expand Up @@ -329,17 +328,17 @@ impl TwoQubitBasisDecomposer {
}
if x12_is_half_pi {
gates.push((StandardGate::SX.into(), smallvec![], smallvec![0]));
global_phase -= PI4;
global_phase -= FRAC_PI_4;
} else if x12_is_non_zero && !x12_is_pi_mult {
if self.pulse_optimize.is_none() {
self.append_1q_sequence(&mut gates, &mut global_phase, rx_matrix(x12).view(), 0);
} else {
return None;
}
}
if abs_diff_eq!(euler_q1[1][1], PI2, epsilon = atol) {
if abs_diff_eq!(euler_q1[1][1], FRAC_PI_2, epsilon = atol) {
gates.push((StandardGate::SX.into(), smallvec![], smallvec![1]));
global_phase -= PI4
global_phase -= FRAC_PI_4
} else if self.pulse_optimize.is_none() {
self.append_1q_sequence(
&mut gates,
Expand All @@ -361,9 +360,9 @@ impl TwoQubitBasisDecomposer {
smallvec![euler_q0[2][1]],
smallvec![0],
));
if abs_diff_eq!(euler_q1[2][1], PI2, epsilon = atol) {
if abs_diff_eq!(euler_q1[2][1], FRAC_PI_2, epsilon = atol) {
gates.push((StandardGate::SX.into(), smallvec![], smallvec![1]));
global_phase -= PI4;
global_phase -= FRAC_PI_4;
} else if self.pulse_optimize.is_none() {
self.append_1q_sequence(
&mut gates,
Expand Down Expand Up @@ -491,7 +490,7 @@ impl TwoQubitBasisDecomposer {
let ipz: ArrayView2<Complex64> = aview2(&IPZ);
let basis_decomposer =
TwoQubitWeylDecomposition::new_inner(gate_matrix, Some(DEFAULT_FIDELITY), None)?;
let super_controlled = relative_eq!(basis_decomposer.a, PI4, max_relative = 1e-09)
let super_controlled = relative_eq!(basis_decomposer.a, FRAC_PI_4, max_relative = 1e-09)
&& relative_eq!(basis_decomposer.c, 0.0, max_relative = 1e-09);

// Create some useful matrices U1, U2, U3 are equivalent to the basis,
Expand Down Expand Up @@ -715,16 +714,6 @@ impl TwoQubitBasisDecomposer {
}
}

static K12R_ARR: GateArray1Q = [
[c64(0., FRAC_1_SQRT_2), c64(FRAC_1_SQRT_2, 0.)],
[c64(-FRAC_1_SQRT_2, 0.), c64(0., -FRAC_1_SQRT_2)],
];

static K12L_ARR: GateArray1Q = [
[c64(0.5, 0.5), c64(0.5, 0.5)],
[c64(-0.5, 0.5), c64(0.5, -0.5)],
];

fn decomp0_inner(target: &TwoQubitWeylDecomposition) -> SmallVec<[Array2<Complex64>; 8]> {
smallvec![target.K1r.dot(&target.K2r), target.K1l.dot(&target.K2l),]
}
Expand Down Expand Up @@ -793,10 +782,10 @@ impl TwoQubitBasisDecomposer {
target.a.sin() * target.b.sin() * target.c.sin(),
),
4. * c64(
(PI4 - target.a).cos()
(FRAC_PI_4 - target.a).cos()
* (self.basis_decomposer.b - target.b).cos()
* target.c.cos(),
(PI4 - target.a).sin()
(FRAC_PI_4 - target.a).sin()
* (self.basis_decomposer.b - target.b).sin()
* target.c.sin(),
),
Expand Down Expand Up @@ -965,31 +954,6 @@ impl TwoQubitBasisDecomposer {
}
}

fn u4_to_su4(u4: ArrayView2<Complex64>) -> (Array2<Complex64>, f64) {
let det_u = ndarray_to_faer(u4).determinant();
let phase_factor = det_u.powf(-0.25).conj();
let su4 = u4.mapv(|x| x / phase_factor);
(su4, phase_factor.arg())
}

fn real_trace_transform(mat: ArrayView2<Complex64>) -> Array2<Complex64> {
let a1 = -mat[[1, 3]] * mat[[2, 0]] + mat[[1, 2]] * mat[[2, 1]] + mat[[1, 1]] * mat[[2, 2]]
- mat[[1, 0]] * mat[[2, 3]];
let a2 = mat[[0, 3]] * mat[[3, 0]] - mat[[0, 2]] * mat[[3, 1]] - mat[[0, 1]] * mat[[3, 2]]
+ mat[[0, 0]] * mat[[3, 3]];
let theta = 0.; // Arbitrary!
let phi = 0.; // This is extra arbitrary!
let psi = f64::atan2(a1.im + a2.im, a1.re - a2.re) - phi;
let im = Complex64::new(0., -1.);
let temp = [
(theta * im).exp(),
(phi * im).exp(),
(psi * im).exp(),
(-(theta + phi + psi) * im).exp(),
];
Array2::from_diag(&arr1(&temp))
}

#[pyfunction]
#[pyo3(name = "two_qubit_decompose_up_to_diagonal")]
pub fn py_two_qubit_decompose_up_to_diagonal(
Expand Down
58 changes: 56 additions & 2 deletions crates/synthesis/src/two_qubit_decompose/common.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,13 @@

use std::f64::consts::FRAC_1_SQRT_2;

use ndarray::{Array2, ArrayView2, array, aview2};
use ndarray::{Array2, ArrayView2, arr1, array, aview2};
use num_complex::{Complex64, ComplexFloat};
use numpy::{IntoPyArray, PyArray2, PyReadonlyArray1, PyReadonlyArray2};
use pyo3::prelude::*;

use crate::linalg::ndarray_to_faer;
use qiskit_circuit::util::{C_ZERO, GateArray2Q, IM, M_IM, c64};
use qiskit_circuit::util::{C_M_ONE, C_ONE, C_ZERO, GateArray1Q, GateArray2Q, IM, M_IM, c64};

pub(super) const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9;

Expand Down Expand Up @@ -97,6 +97,34 @@ static MAGIC_DAGGER: GateArray2Q = [
],
];

pub static K12R_ARR: GateArray1Q = [
[c64(0., FRAC_1_SQRT_2), c64(FRAC_1_SQRT_2, 0.)],
[c64(-FRAC_1_SQRT_2, 0.), c64(0., -FRAC_1_SQRT_2)],
];

pub static K12L_ARR: GateArray1Q = [
[c64(0.5, 0.5), c64(0.5, 0.5)],
[c64(-0.5, 0.5), c64(0.5, -0.5)],
];

pub static B_NON_NORMALIZED: GateArray2Q = [
[C_ONE, IM, C_ZERO, C_ZERO],
[C_ZERO, C_ZERO, IM, C_ONE],
[C_ZERO, C_ZERO, IM, C_M_ONE],
[C_ONE, M_IM, C_ZERO, C_ZERO],
];

pub static B_NON_NORMALIZED_DAGGER: GateArray2Q = [
[c64(0.5, 0.), C_ZERO, C_ZERO, c64(0.5, 0.)],
[c64(0., -0.5), C_ZERO, C_ZERO, c64(0., 0.5)],
[C_ZERO, c64(0., -0.5), c64(0., -0.5), C_ZERO],
[C_ZERO, c64(0.5, 0.), c64(-0.5, 0.), C_ZERO],
];

pub static IPZ: GateArray1Q = [[IM, C_ZERO], [C_ZERO, M_IM]];
pub static IPY: GateArray1Q = [[C_ZERO, C_ONE], [C_M_ONE, C_ZERO]];
pub static IPX: GateArray1Q = [[C_ZERO, IM], [IM, C_ZERO]];

#[inline(always)]
pub(crate) fn transpose_conjugate(mat: ArrayView2<Complex64>) -> Array2<Complex64> {
mat.t().mapv(|x| x.conj())
Expand Down Expand Up @@ -207,3 +235,29 @@ pub fn local_equivalence(weyl: PyReadonlyArray1<f64>) -> PyResult<[f64; 3]> {
- weyl.iter().map(|x| (4. * x).cos()).product::<f64>();
Ok([g0_equiv + 0., g1_equiv + 0., g2_equiv + 0.])
}

/// Convert a 4x4 unitary matrix into a unitary matrix with determinant 1
pub fn u4_to_su4(u4: ArrayView2<Complex64>) -> (Array2<Complex64>, f64) {
let det_u = ndarray_to_faer(u4).determinant();
let phase_factor = det_u.powf(-0.25).conj();
let su4 = u4.mapv(|x| x / phase_factor);
(su4, phase_factor.arg())
}

pub fn real_trace_transform(mat: ArrayView2<Complex64>) -> Array2<Complex64> {
let a1 = -mat[[1, 3]] * mat[[2, 0]] + mat[[1, 2]] * mat[[2, 1]] + mat[[1, 1]] * mat[[2, 2]]
- mat[[1, 0]] * mat[[2, 3]];
let a2 = mat[[0, 3]] * mat[[3, 0]] - mat[[0, 2]] * mat[[3, 1]] - mat[[0, 1]] * mat[[3, 2]]
+ mat[[0, 0]] * mat[[3, 3]];
let theta = 0.; // Arbitrary!
let phi = 0.; // This is extra arbitrary!
let psi = f64::atan2(a1.im + a2.im, a1.re - a2.re) - phi;
let im = Complex64::new(0., -1.);
let temp = [
(theta * im).exp(),
(phi * im).exp(),
(psi * im).exp(),
(-(theta + phi + psi) * im).exp(),
];
Array2::from_diag(&arr1(&temp))
}
6 changes: 2 additions & 4 deletions crates/synthesis/src/two_qubit_decompose/controlled_u_decomposer.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
use approx::abs_diff_eq;
use num_complex::Complex64;
use smallvec::{SmallVec, smallvec};
use std::f64::consts::PI;
use std::f64::consts::FRAC_PI_2;

use ndarray::prelude::*;
use numpy::PyReadonlyArray2;
Expand All @@ -38,8 +38,6 @@ use super::common::DEFAULT_FIDELITY;
use super::gate_sequence::TwoQubitGateSequence;
use super::weyl_decomposition::{Specialization, TwoQubitWeylDecomposition};

const PI2: f64 = PI / 2.;

/// invert 1q gate sequence
fn invert_1q_gate(
gate: (StandardGate, SmallVec<[f64; 3]>),
Expand Down Expand Up @@ -458,7 +456,7 @@ impl TwoQubitControlledUDecomposer {
/// Initialize the KAK decomposition.
pub fn new_inner(rxx_equivalent_gate: RXXEquivalent, euler_basis: &str) -> PyResult<Self> {
let atol = DEFAULT_ATOL;
let test_angles = [0.2, 0.3, PI2];
let test_angles = [0.2, 0.3, FRAC_PI_2];

let scales: PyResult<Vec<f64>> = test_angles
.into_iter()
Expand Down
Loading
⚡