|
| 1 | +// This code is part of Qiskit. |
| 2 | +// |
| 3 | +// (C) Copyright IBM 2026 |
| 4 | +// |
| 5 | +// This code is licensed under the Apache License, Version 2.0. You may |
| 6 | +// obtain a copy of this license in the LICENSE.txt file in the root directory |
| 7 | +// of this source tree or at https://www.apache.org/licenses/LICENSE-2.0. |
| 8 | +// |
| 9 | +// Any modifications or derivative works of this code must retain this |
| 10 | +// copyright notice, and modified files need to carry a notice indicating |
| 11 | +// that they have been altered from the originals. |
| 12 | + |
| 13 | +use pyo3::prelude::*; |
| 14 | +use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, PI}; |
| 15 | + |
| 16 | +use crate::QiskitError; |
| 17 | +use qiskit_circuit::dag_circuit::DAGCircuit; |
| 18 | +use qiskit_circuit::operations::{OperationRef, Param, StandardGate, add_param}; |
| 19 | +use qiskit_synthesis::ross_selinger::gridsynth_rz; |
| 20 | + |
| 21 | +const MINIMUM_EPSILON: f64 = 1e-12; // minimum epsilon for synthesis |
| 22 | +const DEFAULT_ACCURACY: f64 = 1e-10; // default synthesis accuracy |
| 23 | + |
| 24 | +/// Finds a canonical representation of an angle. |
| 25 | +/// |
| 26 | +/// Given `angle`, this returns `(interval, angle_normalized)` such that |
| 27 | +/// angle_normalized = angle (mod pi/2) |
| 28 | +/// (angle - angle_normalized - interval * pi/2) = 0 (mod 4pi) |
| 29 | +/// |
| 30 | +/// The canonical representation is limited by the f64 representation. |
| 31 | +/// Any angle that differs in decimal places beyond f64 will be non-unique. |
| 32 | +fn canonicalize_angle(angle: f64) -> (u8, f64) { |
| 33 | + let angle_normalized = angle.rem_euclid(FRAC_PI_2); |
| 34 | + let interval = ((angle - angle_normalized) / FRAC_PI_2) |
| 35 | + .round() |
| 36 | + .rem_euclid(8.) as u8; |
| 37 | + (interval, angle_normalized) |
| 38 | +} |
| 39 | + |
| 40 | +/// Lookup table for fixing the circuit based on interval computed during |
| 41 | +/// canonicalization. |
| 42 | +/// |
| 43 | +/// The table is based on the properties of `Rz(theta)` such as: |
| 44 | +/// * `Rz(theta + pi/2) = Rz(theta).S`, up to a global phase of `-pi/4`, |
| 45 | +/// * `Rz(theta + pi) = Rz(theta).Z`, up to a global phase of `-pi/2`, |
| 46 | +/// * `Rz(theta + 2*pi) = -Rz(theta)`, up to a global phase of `-pi`. |
| 47 | +static PHASE_GATE_LUT: [(f64, Option<StandardGate>); 8] = [ |
| 48 | + (0.0, None), |
| 49 | + (-FRAC_PI_4, Some(StandardGate::S)), |
| 50 | + (-FRAC_PI_2, Some(StandardGate::Z)), |
| 51 | + (-3. * FRAC_PI_4, Some(StandardGate::Sdg)), |
| 52 | + (PI, None), |
| 53 | + (-5. * FRAC_PI_4, Some(StandardGate::S)), |
| 54 | + (-6. * FRAC_PI_4, Some(StandardGate::Z)), |
| 55 | + (-7. * FRAC_PI_4, Some(StandardGate::Sdg)), |
| 56 | +]; |
| 57 | + |
| 58 | +/// Approximates RZ-rotation using gridsynth. |
| 59 | +/// |
| 60 | +/// Returns the sequence of gates in the synthesized circuit and |
| 61 | +/// an update to the global phase. |
| 62 | +fn synthesize_rz_gate_via_gridsynth( |
| 63 | + angle: f64, |
| 64 | + epsilon: f64, |
| 65 | +) -> PyResult<(Vec<StandardGate>, Param)> { |
| 66 | + let circ_data = gridsynth_rz(angle, epsilon)?; |
| 67 | + |
| 68 | + // obtain phase from circuit data |
| 69 | + let phase = circ_data.global_phase().clone(); |
| 70 | + |
| 71 | + // get sequence of standard gates |
| 72 | + let sequence: Vec<StandardGate> = circ_data |
| 73 | + .data() |
| 74 | + .iter() |
| 75 | + .map(|inst| { |
| 76 | + if let OperationRef::StandardGate(gate) = inst.op.view() { |
| 77 | + gate |
| 78 | + } else { |
| 79 | + unreachable!("gridsynth only produces standard gates"); |
| 80 | + } |
| 81 | + }) |
| 82 | + .collect(); |
| 83 | + Ok((sequence, phase)) |
| 84 | +} |
| 85 | + |
| 86 | +/// Synthesize RZ gates in the circuit, modifying the circuit in-place. |
| 87 | +/// |
| 88 | +/// # Arguments |
| 89 | +/// |
| 90 | +/// - `dag`: The DAG circuit in which the RZ gates will be synthesized. |
| 91 | +/// - `approximation_degree`: Controls the overall degree of approximation. |
| 92 | +/// - `synthesis_error`: Maximum allowed error for the approximate synthesis of |
| 93 | +/// :math:`RZ(\theta)`. |
| 94 | +/// - `cache_error`: Maximum allowed error when reusing a cached synthesis |
| 95 | +/// result for angles close to :math:`\theta`. |
| 96 | +/// |
| 97 | +/// If both `synthesis_error` and `cache_error` are provided, they specify the error budget |
| 98 | +/// due to approximate synthesis and due to caching respectively. If either value is not |
| 99 | +/// specified, the total allowed error is derived from `approximation_degree`, and |
| 100 | +/// suitable values for `synthesis_error` and `cache_error` are computed automatically. |
| 101 | +#[pyfunction] |
| 102 | +#[pyo3(name = "synthesize_rz_rotations")] |
| 103 | +#[pyo3(signature = (dag, approximation_degree=None, synthesis_error=None, cache_error=None))] |
| 104 | +pub fn py_run_synthesize_rz_rotations( |
| 105 | + dag: &mut DAGCircuit, |
| 106 | + approximation_degree: Option<f64>, |
| 107 | + synthesis_error: Option<f64>, |
| 108 | + cache_error: Option<f64>, |
| 109 | +) -> PyResult<()> { |
| 110 | + // Skip the pass if there are no RZ rotation gates. |
| 111 | + if dag.get_op_counts().keys().all(|k| k != "rz") { |
| 112 | + return Ok(()); |
| 113 | + } |
| 114 | + |
| 115 | + // Compute error budgets. When approximation degree is used, the total error is |
| 116 | + // computed as 1 - approximation_degree, and the error budget for synthesis and for |
| 117 | + // caching are distributed equally. |
| 118 | + let (synthesis_error, cache_error) = match (synthesis_error, cache_error) { |
| 119 | + (Some(synthesis_error), Some(cache_error)) => (synthesis_error, cache_error), |
| 120 | + _ => { |
| 121 | + let total_error = if let Some(approximation_degree) = approximation_degree { |
| 122 | + MINIMUM_EPSILON.max(1. - approximation_degree) |
| 123 | + } else { |
| 124 | + DEFAULT_ACCURACY |
| 125 | + }; |
| 126 | + (total_error / 2., total_error / 2.) |
| 127 | + } |
| 128 | + }; |
| 129 | + |
| 130 | + // By an explicit computation one can show that if the current angle is within |
| 131 | + // 4.0 * arcsin(cache_error / 2) from the previous angle, the error due to reusing the synthesis |
| 132 | + // result for the previous angle is precisely cache_error. Contact a Qiskit synthesis developer |
| 133 | + // for more details! |
| 134 | + let bin_width = 4. * (cache_error / 2.).asin(); |
| 135 | + |
| 136 | + // Iterate over nodes in the DAG and collect nodes that have RZ gates. |
| 137 | + // Canonicalize angles already at this stage, so that we can use them for sorting. |
| 138 | + let mut candidates: Vec<_> = dag |
| 139 | + .op_nodes(false) |
| 140 | + .filter_map(|(node_index, inst)| { |
| 141 | + if let OperationRef::StandardGate(StandardGate::RZ) = inst.op.view() { |
| 142 | + if let Param::Float(angle) = inst.params_view()[0] { |
| 143 | + let (interval_index, canonical_angle) = canonicalize_angle(angle); |
| 144 | + Some((node_index, canonical_angle, interval_index)) |
| 145 | + } else { |
| 146 | + None |
| 147 | + } |
| 148 | + } else { |
| 149 | + None |
| 150 | + } |
| 151 | + }) |
| 152 | + .collect(); |
| 153 | + |
| 154 | + // Sort candidates based on the canonicalized angles |
| 155 | + candidates.sort_unstable_by(|a, b| { |
| 156 | + a.1.partial_cmp(&b.1) |
| 157 | + .expect("Angles are never NaN here, so we can compare f64.") |
| 158 | + }); |
| 159 | + |
| 160 | + let mut prev_result: Option<(f64, (Vec<StandardGate>, Param))> = None; |
| 161 | + |
| 162 | + for (node_index, angle, interval_index) in candidates { |
| 163 | + // Get or compute the sequence and phase update. |
| 164 | + let should_recompute = prev_result |
| 165 | + .as_ref() |
| 166 | + .is_none_or(|(prev_angle, _)| *prev_angle + bin_width < angle); |
| 167 | + |
| 168 | + if should_recompute { |
| 169 | + let (sequence, phase_update) = synthesize_rz_gate_via_gridsynth(angle, synthesis_error) |
| 170 | + .map_err(|e| QiskitError::new_err(e.to_string()))?; |
| 171 | + |
| 172 | + prev_result = Some((angle, (sequence, phase_update))); |
| 173 | + } |
| 174 | + |
| 175 | + let (sequence, phase_update) = &prev_result |
| 176 | + .as_ref() |
| 177 | + .expect("is_none_or ensures prev_result is never None") |
| 178 | + .1; |
| 179 | + |
| 180 | + // Add the gates and phase update to DAG, remove old node |
| 181 | + for new_gate in sequence { |
| 182 | + dag.insert_1q_on_incoming_qubit((*new_gate, &[]), node_index); |
| 183 | + } |
| 184 | + if let Some(gate) = PHASE_GATE_LUT[interval_index as usize].1 { |
| 185 | + dag.insert_1q_on_incoming_qubit((gate, &[]), node_index); |
| 186 | + } |
| 187 | + dag.remove_1q_sequence(&[node_index]); |
| 188 | + |
| 189 | + let phase_update_with_shift = |
| 190 | + add_param(phase_update, PHASE_GATE_LUT[interval_index as usize].0); |
| 191 | + dag.add_global_phase(&phase_update_with_shift)?; |
| 192 | + } |
| 193 | + |
| 194 | + Ok(()) |
| 195 | +} |
| 196 | + |
| 197 | +pub fn synthesize_rz_rotations_mod(m: &Bound<PyModule>) -> PyResult<()> { |
| 198 | + m.add_wrapped(wrap_pyfunction!(py_run_synthesize_rz_rotations))?; |
| 199 | + Ok(()) |
| 200 | +} |
0 commit comments