wai-quantum 0.5.0

//! Native programs: a circuit in a processor's own gates, on its own qubits.
//!
//! A cloud processor runs only its basis gates, and a two-qubit gate only on a
//! pair its coupling map connects. [`compile`] takes a circuit there in four
//! exact steps:
//!
//! 1. [`super::device::lower`] to the standard gates;
//! 2. each standard gate to `rz`, `sx`, `x` and `cz` by textbook identities
//!    (`H = Rz(π/2)·SX·Rz(π/2)`, `CX = H·CZ·H`, the controlled phase through two
//!    `CX` and half-angle phases, the Toffoli through seven `T`s), every angle a
//!    dyadic multiple of `2π`;
//! 3. onto physical qubits: an initial layout, then for each `cz` on a pair the
//!    map does not connect, `SWAP`s along a shortest path;
//! 4. into the [`Basis`]: consecutive `rz` merged, or, for phased-X processors,
//!    every `rz` carried forward into the phase of the next rotation.
//!
//! Single-qubit gates are replaced up to a global phase, which no measurement
//! sees. A phased-X program leaves each qubit's carried `Z` rotation unapplied at
//! the end ([`NativeProgram::residual`]): a `Z` rotation just before a
//! computational-basis measurement changes no outcome.
//!
//! # Layout
//!
//! Each physical qubit in turn seeds a breadth-first layout of the circuit's
//! qubits over the coupling map; the layout whose routing needs the fewest `SWAP`s
//! wins, the lowest seed breaking ties. Deterministic: the same circuit and map
//! give the same program on every machine.
//!
//! # Verification
//!
//! [`NativeProgram::verify`] checks the program against its circuit before it is
//! sent anywhere: every `cz` on a connected pair, and the program's final state
//! (with the residual rotations applied) equal to the circuit's, with each logical
//! qubit on the physical qubit it ends on, by this crate's own simulator.

use super::device::{lower, DeviceError, StdGate, Turn};
use crate::quantum::{BaseGate, Circuit, Gate, DYADIC_MAX, MAX_QUBITS, ONE};

/// The angle `2π·m/2^k`, reduced: `m < 2^k`, and `m` odd unless the angle is
/// zero (`m = k = 0`). `k ≤ DYADIC_MAX`.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub struct Dyadic {
    m: u64,
    k: u8,
}

impl Dyadic {
    pub const ZERO: Dyadic = Dyadic { m: 0, k: 0 };

    /// `2π·m/2^k` reduced modulo `2π`, or `None` past `k = DYADIC_MAX`.
    pub fn new(m: u64, k: u8) -> Option<Self> {
        (usize::from(k) <= DYADIC_MAX).then(|| Dyadic::reduce(m, k))
    }

    fn reduce(m: u64, k: u8) -> Self {
        let (mut m, mut k) = (m & ((1u64 << k) - 1), k);
        if m == 0 {
            return Dyadic::ZERO;
        }
        while m & 1 == 0 {
            m >>= 1;
            k -= 1;
        }
        Dyadic { m, k }
    }

    /// The angle as `(m, k)`: `m/2^k` of a full turn, which a service taking
    /// turns reads exactly.
    pub fn turns(self) -> (u64, u8) {
        (self.m, self.k)
    }

    pub fn is_zero(self) -> bool {
        self.m == 0
    }

    /// Half this representative in `[0, 2π)`, or `None` past `DYADIC_MAX`.
    pub fn half(self) -> Option<Dyadic> {
        if self.m == 0 {
            return Some(self);
        }
        (usize::from(self.k) < DYADIC_MAX).then_some(Dyadic { m: self.m, k: self.k + 1 })
    }

    /// The angle in `(−π, π]` as OpenQASM writes it: `0`, `pi`, `pi/4`,
    /// `-3*pi/8`.
    pub fn qasm(self) -> String {
        match self.k {
            0 => "0".into(),
            1 => "pi".into(),
            k => {
                let full = 1u64 << k;
                let (sign, m) = if self.m > full / 2 { ("-", full - self.m) } else { ("", self.m) };
                let den = 1u64 << (k - 1);
                if m == 1 {
                    format!("{sign}pi/{den}")
                } else {
                    format!("{sign}{m}*pi/{den}")
                }
            }
        }
    }
}

impl core::ops::Add for Dyadic {
    type Output = Dyadic;
    fn add(self, o: Dyadic) -> Dyadic {
        let k = self.k.max(o.k);
        Dyadic::reduce((self.m << (k - self.k)) + (o.m << (k - o.k)), k)
    }
}

impl core::ops::Neg for Dyadic {
    type Output = Dyadic;
    fn neg(self) -> Dyadic {
        if self.m == 0 {
            self
        } else {
            Dyadic { m: (1u64 << self.k) - self.m, k: self.k }
        }
    }
}

impl From<Turn> for Dyadic {
    fn from(t: Turn) -> Dyadic {
        let d = Dyadic::reduce(1, t.k);
        if t.negative {
            -d
        } else {
            d
        }
    }
}

/// A processor's gate basis.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Basis {
    /// `rz(θ)`, `sx`, `x` and `cz`.
    RzSxXCz,
    /// `prx(θ, φ) = Rz(φ)·Rx(θ)·Rz(−φ)` and `cz`, the `Z` rotations carried in
    /// the phases.
    PrxCz,
    /// `gpi(φ)`, `gpi2(φ)` and `zz(θ)`, with angles in turns: `gpi2(φ)` is
    /// `prx(π/2, φ)`, `gpi(φ)` is `prx(π, φ)` up to a global phase, and
    /// `zz(θ) = exp(−iπθ·Z⊗Z)`. A `cz` is `zz(¼)` with a `Z` rotation of `−π/2`
    /// on each qubit, carried in the phases like every other.
    GpiZz,
}

/// A native gate on physical qubits.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum NativeGate {
    Rz(u16, Dyadic),
    Sx(u16),
    X(u16),
    /// `prx(θ, φ)`.
    Prx(u16, Dyadic, Dyadic),
    Cz(u16, u16),
    /// `gpi(φ)`, the phase a [`Dyadic`] read as turns.
    Gpi(u16, Dyadic),
    /// `gpi2(φ)`.
    Gpi2(u16, Dyadic),
    /// `zz(θ)` on a pair, the angle in turns.
    Zz(u16, u16, Dyadic),
}

/// A processor's connectivity: which physical pairs a two-qubit gate may act
/// on, in either order.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Coupling {
    adj: Vec<Vec<u16>>,
}

impl Coupling {
    /// `n` physical qubits and the pairs connected. A pair may be listed in
    /// either order or both; a qubit paired with itself or outside `0..n` is
    /// refused.
    pub fn new(n: u16, pairs: &[(u16, u16)]) -> Result<Self, NativeError> {
        let mut adj = vec![Vec::new(); usize::from(n)];
        for &(a, b) in pairs {
            if a == b || a >= n || b >= n {
                return Err(NativeError::Coupling(format!("pair ({a}, {b}) on {n} qubits")));
            }
            adj[usize::from(a)].push(b);
            adj[usize::from(b)].push(a);
        }
        for list in &mut adj {
            list.sort_unstable();
            list.dedup();
        }
        Ok(Coupling { adj })
    }

    /// `n` qubits in a line.
    pub fn line(n: u16) -> Self {
        let pairs: Vec<(u16, u16)> = (1..n).map(|i| (i - 1, i)).collect();
        Coupling::new(n, &pairs).expect("a line is a valid map")
    }

    /// A `rows × cols` grid, row-major.
    pub fn grid(rows: u16, cols: u16) -> Self {
        let mut pairs = Vec::new();
        for r in 0..rows {
            for c in 0..cols {
                let q = r * cols + c;
                if c + 1 < cols {
                    pairs.push((q, q + 1));
                }
                if r + 1 < rows {
                    pairs.push((q, q + cols));
                }
            }
        }
        Coupling::new(rows * cols, &pairs).expect("a grid is a valid map")
    }

    pub fn qubits(&self) -> u16 {
        self.adj.len() as u16
    }

    pub fn connected(&self, a: u16, b: u16) -> bool {
        self.adj.get(usize::from(a)).is_some_and(|l| l.binary_search(&b).is_ok())
    }

    /// Breadth-first order from `s`, neighbours in index order.
    fn bfs(&self, s: u16) -> Vec<u16> {
        let mut seen = vec![false; self.adj.len()];
        let mut order = vec![s];
        seen[usize::from(s)] = true;
        let mut i = 0;
        while i < order.len() {
            for &v in &self.adj[usize::from(order[i])] {
                if !seen[usize::from(v)] {
                    seen[usize::from(v)] = true;
                    order.push(v);
                }
            }
            i += 1;
        }
        order
    }

    /// For each target `b`, each qubit's next step towards `b` on a shortest
    /// path, and its distance; `u16::MAX` where `b` is unreachable.
    fn paths(&self) -> (Vec<Vec<u16>>, Vec<Vec<u16>>) {
        let n = self.adj.len();
        let mut next = vec![vec![u16::MAX; n]; n];
        let mut dist = vec![vec![u16::MAX; n]; n];
        for b in 0..n {
            dist[b][b] = 0;
            next[b][b] = b as u16;
            let mut queue = vec![b as u16];
            let mut i = 0;
            while i < queue.len() {
                let u = queue[i];
                for &v in &self.adj[usize::from(u)] {
                    if dist[b][usize::from(v)] == u16::MAX {
                        dist[b][usize::from(v)] = dist[b][usize::from(u)] + 1;
                        next[b][usize::from(v)] = u;
                        queue.push(v);
                    }
                }
                i += 1;
            }
        }
        (next, dist)
    }
}

/// A circuit compiled for a processor.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct NativeProgram {
    pub basis: Basis,
    /// The circuit's width.
    pub n_qubits: u8,
    /// Logical qubit `i` starts on physical qubit `layout[i]`.
    pub layout: Vec<u16>,
    /// Logical qubit `i` ends, and is measured into `c[i]`, on
    /// `final_layout[i]`.
    pub final_layout: Vec<u16>,
    pub gates: Vec<NativeGate>,
    /// The `Z` rotation still owed on each physical qubit when the gates end,
    /// as `(qubit, angle)`: none for [`Basis::RzSxXCz`]. A `Z` rotation just
    /// before a computational-basis measurement changes no outcome, so it is
    /// never sent.
    pub residual: Vec<(u16, Dyadic)>,
    /// The `SWAP`s routing inserted.
    pub swaps: usize,
}

/// Why a circuit cannot be compiled for a processor, or a program fails its
/// verification.
#[derive(Debug, PartialEq, Eq)]
pub enum NativeError {
    Device(DeviceError),
    Coupling(String),
    /// No connected part of the map has as many qubits as the circuit.
    TooFewQubits { need: u8, largest: usize },
    /// A controlled phase would need a turn finer than `2π/2^DYADIC_MAX`.
    TooFine,
    /// Gate `gate` is not in the program's basis, or not on a connected pair.
    NotNative { gate: usize },
    /// The program touches more qubits than this crate simulates.
    TooWideToVerify { qubits: usize },
    /// The program's state is not the circuit's: the fixed-point fidelity.
    Mismatch { fidelity: i64 },
}

impl core::fmt::Display for NativeError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            NativeError::Device(e) => write!(f, "{e}"),
            NativeError::Coupling(s) => write!(f, "coupling map: {s}"),
            NativeError::TooFewQubits { need, largest } => {
                write!(f, "the circuit needs {need} connected qubits; the largest connected part has {largest}")
            }
            NativeError::TooFine => write!(f, "a controlled phase would need a turn finer than 2π/2^{DYADIC_MAX}"),
            NativeError::NotNative { gate } => write!(f, "gate {gate} is not native to the program's processor"),
            NativeError::TooWideToVerify { qubits } => {
                write!(f, "the program touches {qubits} qubits; at most {MAX_QUBITS} are simulated")
            }
            NativeError::Mismatch { fidelity } => write!(f, "the program's state differs from the circuit's (fidelity {fidelity} of {ONE})"),
        }
    }
}

impl std::error::Error for NativeError {}

impl From<DeviceError> for NativeError {
    fn from(e: DeviceError) -> Self {
        NativeError::Device(e)
    }
}

/// An operation in the `rz`/`sx`/`x`/`cz` set on qubit indices: logical before
/// routing, physical after.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Op {
    Rz(u16, Dyadic),
    Sx(u16),
    X(u16),
    Cz(u16, u16),
}

const fn dy(m: u64, k: u8) -> Dyadic {
    Dyadic { m, k }
}
/// `π`, `π/2`, `−π/2`, `π/4`, `−π/4` as reduced turns.
const PI: Dyadic = dy(1, 1);
const HALF_PI: Dyadic = dy(1, 2);
const NEG_HALF_PI: Dyadic = dy(3, 2);
const QUARTER_PI: Dyadic = dy(1, 3);
const NEG_QUARTER_PI: Dyadic = dy(7, 3);

fn h(o: &mut Vec<Op>, q: u16) {
    o.extend([Op::Rz(q, HALF_PI), Op::Sx(q), Op::Rz(q, HALF_PI)]);
}

fn cx(o: &mut Vec<Op>, c: u16, t: u16) {
    h(o, t);
    o.push(Op::Cz(c, t));
    h(o, t);
}

/// `Ry(θ) = Rz(π/2)·Rx(θ)·Rz(−π/2)`, `Rx(θ) = H·Rz(θ)·H`.
fn ry(o: &mut Vec<Op>, q: u16, theta: Dyadic) {
    o.push(Op::Rz(q, NEG_HALF_PI));
    h(o, q);
    o.push(Op::Rz(q, theta));
    h(o, q);
    o.push(Op::Rz(q, HALF_PI));
}

fn decompose(gates: &[StdGate]) -> Result<Vec<Op>, NativeError> {
    let mut o = Vec::new();
    for g in gates {
        match *g {
            StdGate::X(q) => o.push(Op::X(q.into())),
            // Y = i·X·Z.
            StdGate::Y(q) => o.extend([Op::Rz(q.into(), PI), Op::X(q.into())]),
            StdGate::Z(q) => o.push(Op::Rz(q.into(), PI)),
            StdGate::H(q) => h(&mut o, q.into()),
            StdGate::S(q) => o.push(Op::Rz(q.into(), HALF_PI)),
            StdGate::Sdg(q) => o.push(Op::Rz(q.into(), NEG_HALF_PI)),
            StdGate::T(q) => o.push(Op::Rz(q.into(), QUARTER_PI)),
            StdGate::Tdg(q) => o.push(Op::Rz(q.into(), NEG_QUARTER_PI)),
            StdGate::P(q, t) => o.push(Op::Rz(q.into(), t.into())),
            StdGate::Ry(q, t) => ry(&mut o, q.into(), t.into()),
            StdGate::Cx(c, t) => cx(&mut o, c.into(), t.into()),
            StdGate::Cy(c, t) => {
                o.push(Op::Rz(t.into(), NEG_HALF_PI));
                cx(&mut o, c.into(), t.into());
                o.push(Op::Rz(t.into(), HALF_PI));
            }
            StdGate::Cz(c, t) => o.push(Op::Cz(c.into(), t.into())),
            // H = Ry(π/4)·Z·Ry(−π/4).
            StdGate::Ch(c, t) => {
                ry(&mut o, t.into(), NEG_QUARTER_PI);
                o.push(Op::Cz(c.into(), t.into()));
                ry(&mut o, t.into(), QUARTER_PI);
            }
            // CP(θ) = P(θ/2)_c · CX · P(−θ/2)_t · CX · P(θ/2)_t.
            StdGate::Cp(c, t, a) => {
                let half = Dyadic::from(a).half().ok_or(NativeError::TooFine)?;
                let (c, t) = (u16::from(c), u16::from(t));
                o.push(Op::Rz(c, half));
                cx(&mut o, c, t);
                o.push(Op::Rz(t, -half));
                cx(&mut o, c, t);
                o.push(Op::Rz(t, half));
            }
            // The Toffoli through six CX and seven T gates.
            StdGate::Ccx(a, b, t) => {
                let (a, b, t) = (u16::from(a), u16::from(b), u16::from(t));
                h(&mut o, t);
                cx(&mut o, b, t);
                o.push(Op::Rz(t, NEG_QUARTER_PI));
                cx(&mut o, a, t);
                o.push(Op::Rz(t, QUARTER_PI));
                cx(&mut o, b, t);
                o.push(Op::Rz(t, NEG_QUARTER_PI));
                cx(&mut o, a, t);
                o.push(Op::Rz(b, QUARTER_PI));
                o.push(Op::Rz(t, QUARTER_PI));
                h(&mut o, t);
                cx(&mut o, a, b);
                o.push(Op::Rz(a, QUARTER_PI));
                o.push(Op::Rz(b, NEG_QUARTER_PI));
                cx(&mut o, a, b);
            }
        }
    }
    Ok(o)
}

/// Route `ops` (on logical qubits) from `layout`; the physical operations, the
/// final layout and the `SWAP` count.
fn route(ops: &[Op], layout: &[u16], next: &[Vec<u16>], dist: &[Vec<u16>], phys: usize) -> (Vec<Op>, Vec<u16>, usize) {
    let mut at: Vec<Option<u16>> = vec![None; phys];
    let mut on = layout.to_vec();
    for (l, &p) in layout.iter().enumerate() {
        at[usize::from(p)] = Some(l as u16);
    }
    let mut out = Vec::with_capacity(ops.len());
    let mut swaps = 0;
    for op in ops {
        match *op {
            Op::Rz(q, a) => out.push(Op::Rz(on[usize::from(q)], a)),
            Op::Sx(q) => out.push(Op::Sx(on[usize::from(q)])),
            Op::X(q) => out.push(Op::X(on[usize::from(q)])),
            Op::Cz(a, b) => {
                let pb = on[usize::from(b)];
                let mut pa = on[usize::from(a)];
                while dist[usize::from(pb)][usize::from(pa)] > 1 {
                    let step = next[usize::from(pb)][usize::from(pa)];
                    // SWAP = three CX.
                    cx(&mut out, pa, step);
                    cx(&mut out, step, pa);
                    cx(&mut out, pa, step);
                    swaps += 1;
                    let moved = at[usize::from(step)];
                    at[usize::from(step)] = Some(a);
                    at[usize::from(pa)] = moved;
                    if let Some(m) = moved {
                        on[usize::from(m)] = pa;
                    }
                    on[usize::from(a)] = step;
                    pa = step;
                }
                out.push(Op::Cz(pa, pb));
            }
        }
    }
    (out, on, swaps)
}

/// Compile `c` for a processor with basis `basis` and connectivity
/// `coupling`. See the module documentation.
pub fn compile(c: &Circuit, basis: Basis, coupling: &Coupling) -> Result<NativeProgram, NativeError> {
    let ops = decompose(&lower(c)?)?;
    let n = usize::from(c.n_qubits);
    let phys = usize::from(coupling.qubits());
    let (next, dist) = coupling.paths();
    // (SWAPs, layout, routed operations, final layout) of the best seed so far.
    type Routed = (usize, Vec<u16>, Vec<Op>, Vec<u16>);
    let mut best: Option<Routed> = None;
    let mut largest = 0;
    for seed in 0..coupling.qubits() {
        let order = coupling.bfs(seed);
        largest = largest.max(order.len());
        if order.len() < n {
            continue;
        }
        let layout = order[..n].to_vec();
        let (routed, fin, swaps) = route(&ops, &layout, &next, &dist, phys);
        if best.as_ref().is_none_or(|b| swaps < b.0) {
            best = Some((swaps, layout, routed, fin));
        }
    }
    let Some((swaps, layout, routed, final_layout)) = best else {
        return Err(NativeError::TooFewQubits { need: c.n_qubits, largest });
    };
    let (gates, residual) = to_basis(&routed, basis, phys);
    Ok(NativeProgram { basis, n_qubits: c.n_qubits, layout, final_layout, gates, residual, swaps })
}

fn to_basis(ops: &[Op], basis: Basis, phys: usize) -> (Vec<NativeGate>, Vec<(u16, Dyadic)>) {
    // The Z rotation pending on each physical qubit: merged into the next `rz`
    // flush, or carried into the next rotation's phase.
    let mut frame = vec![Dyadic::ZERO; phys];
    let mut out = Vec::with_capacity(ops.len());
    let flush = |out: &mut Vec<NativeGate>, frame: &mut [Dyadic], q: u16| {
        let f = std::mem::replace(&mut frame[usize::from(q)], Dyadic::ZERO);
        if !f.is_zero() {
            out.push(NativeGate::Rz(q, f));
        }
    };
    for op in ops {
        match (basis, *op) {
            (_, Op::Rz(q, a)) => frame[usize::from(q)] = frame[usize::from(q)] + a,
            // CZ = Rz(−π/2)⊗Rz(−π/2)·zz(¼) up to a global phase; the rotations
            // join the frames.
            (Basis::GpiZz, Op::Cz(a, b)) => {
                out.push(NativeGate::Zz(a, b, HALF_PI));
                frame[usize::from(a)] = frame[usize::from(a)] + NEG_HALF_PI;
                frame[usize::from(b)] = frame[usize::from(b)] + NEG_HALF_PI;
            }
            // A Z rotation commutes with CZ: it waits.
            (_, Op::Cz(a, b)) => out.push(NativeGate::Cz(a, b)),
            (Basis::RzSxXCz, Op::Sx(q)) => {
                flush(&mut out, &mut frame, q);
                out.push(NativeGate::Sx(q));
            }
            (Basis::RzSxXCz, Op::X(q)) => {
                flush(&mut out, &mut frame, q);
                out.push(NativeGate::X(q));
            }
            // prx(θ, 0)·Rz(α) = Rz(α)·prx(θ, −α): the rotation moves on. SX and X
            // are Rx(π/2) and Rx(π) up to a global phase.
            (Basis::PrxCz, Op::Sx(q)) => out.push(NativeGate::Prx(q, HALF_PI, -frame[usize::from(q)])),
            (Basis::PrxCz, Op::X(q)) => out.push(NativeGate::Prx(q, PI, -frame[usize::from(q)])),
            (Basis::GpiZz, Op::Sx(q)) => out.push(NativeGate::Gpi2(q, -frame[usize::from(q)])),
            (Basis::GpiZz, Op::X(q)) => out.push(NativeGate::Gpi(q, -frame[usize::from(q)])),
        }
    }
    let mut residual = Vec::new();
    for q in 0..phys as u16 {
        match basis {
            Basis::RzSxXCz => flush(&mut out, &mut frame, q),
            Basis::PrxCz | Basis::GpiZz => {
                if !frame[usize::from(q)].is_zero() {
                    residual.push((q, frame[usize::from(q)]));
                }
            }
        }
    }
    (out, residual)
}

/// `prx(θ, φ) = Rz(φ)·Rx(θ)·Rz(−φ)` on `q[i]`, in the gates the reader takes.
fn phased_x(i: usize, theta: Dyadic, phi: Dyadic) -> String {
    format!("rz({}) q[{i}];\nrx({}) q[{i}];\nrz({}) q[{i}];\n", (-phi).qasm(), theta.qasm(), phi.qasm())
}

/// A [`Dyadic`] as a number of turns, `m/2^k`, written exactly in decimal.
pub fn turns(d: Dyadic) -> String {
    let (m, k) = d.turns();
    if k == 0 {
        return "0".into();
    }
    // m/2^k = m·5^k/10^k: k decimal digits, exactly.
    let digits = (u128::from(m) * 5u128.pow(u32::from(k))).to_string();
    let k = usize::from(k);
    let padded = format!("{digits:0>k$}");
    let (int, frac) = padded.split_at(padded.len() - k);
    format!("{}.{}", if int.is_empty() { "0" } else { int }, frac.trim_end_matches('0'))
}

impl NativeProgram {
    /// The program as OpenQASM 3 on physical qubits (`$p`), each logical
    /// qubit `i` measured into `c[i]` where it ends. A phased-X program defines
    /// `prx` from the standard gates.
    pub fn qasm3(&self) -> String {
        let mut s = String::from("OPENQASM 3.0;\ninclude \"stdgates.inc\";\n");
        match self.basis {
            Basis::RzSxXCz => {}
            Basis::PrxCz => s.push_str("gate prx(theta, phi) a { rz(-phi) a; rx(theta) a; rz(phi) a; }\n"),
            // The phases and the angle in turns.
            Basis::GpiZz => s.push_str(
                "gate gpi(phi) a { rz(-2*pi*phi) a; x a; rz(2*pi*phi) a; }\n\
                 gate gpi2(phi) a { rz(-2*pi*phi) a; rx(pi/2) a; rz(2*pi*phi) a; }\n\
                 gate zz(theta) a, b { cx a, b; rz(2*pi*theta) b; cx a, b; }\n",
            ),
        }
        s.push_str(&format!("bit[{}] c;\n", self.n_qubits));
        for g in &self.gates {
            s.push_str(&match *g {
                NativeGate::Rz(q, a) => format!("rz({}) ${q};\n", a.qasm()),
                NativeGate::Sx(q) => format!("sx ${q};\n"),
                NativeGate::X(q) => format!("x ${q};\n"),
                NativeGate::Prx(q, t, p) => format!("prx({}, {}) ${q};\n", t.qasm(), p.qasm()),
                NativeGate::Cz(a, b) => format!("cz ${a}, ${b};\n"),
                NativeGate::Gpi(q, p) => format!("gpi({}) ${q};\n", turns(p)),
                NativeGate::Gpi2(q, p) => format!("gpi2({}) ${q};\n", turns(p)),
                NativeGate::Zz(a, b, t) => format!("zz({}) ${a}, ${b};\n", turns(t)),
            });
        }
        for (i, p) in self.final_layout.iter().enumerate() {
            s.push_str(&format!("c[{i}] = measure ${p};\n"));
        }
        s
    }

    /// Check the program against the circuit it was compiled from: every gate
    /// in the basis, every `cz` on a pair `coupling` connects, and the final
    /// state, residual rotations applied, equal to the circuit's with each
    /// logical qubit on its final physical qubit.
    pub fn verify(&self, c: &Circuit, coupling: &Coupling) -> Result<(), NativeError> {
        for (i, g) in self.gates.iter().enumerate() {
            let ok = match (self.basis, *g) {
                (Basis::RzSxXCz | Basis::PrxCz, NativeGate::Cz(a, b)) | (Basis::GpiZz, NativeGate::Zz(a, b, _)) => coupling.connected(a, b),
                (Basis::RzSxXCz, NativeGate::Rz(q, _) | NativeGate::Sx(q) | NativeGate::X(q)) => q < coupling.qubits(),
                (Basis::PrxCz, NativeGate::Prx(q, _, _)) => q < coupling.qubits(),
                (Basis::GpiZz, NativeGate::Gpi(q, _) | NativeGate::Gpi2(q, _)) => q < coupling.qubits(),
                _ => false,
            };
            if !ok {
                return Err(NativeError::NotNative { gate: i });
            }
        }
        // The physical qubits the program touches, renumbered from zero.
        let mut used: Vec<u16> = self.final_layout.clone();
        for g in &self.gates {
            match *g {
                NativeGate::Rz(q, _) | NativeGate::Sx(q) | NativeGate::X(q) | NativeGate::Prx(q, _, _) | NativeGate::Gpi(q, _) | NativeGate::Gpi2(q, _) => used.push(q),
                NativeGate::Cz(a, b) | NativeGate::Zz(a, b, _) => used.extend([a, b]),
            }
        }
        used.sort_unstable();
        used.dedup();
        if used.len() > usize::from(MAX_QUBITS) {
            return Err(NativeError::TooWideToVerify { qubits: used.len() });
        }
        let idx = |p: u16| used.binary_search(&p).expect("every touched qubit is listed");
        let m = used.len();
        let mut text = format!("OPENQASM 3.0;\nqubit[{m}] q;\n");
        for g in &self.gates {
            text.push_str(&match *g {
                NativeGate::Rz(q, a) => format!("rz({}) q[{}];\n", a.qasm(), idx(q)),
                NativeGate::Sx(q) => format!("sx q[{}];\n", idx(q)),
                NativeGate::X(q) => format!("x q[{}];\n", idx(q)),
                NativeGate::Prx(q, t, p) => phased_x(idx(q), t, p),
                NativeGate::Gpi(q, p) => phased_x(idx(q), PI, p),
                NativeGate::Gpi2(q, p) => phased_x(idx(q), HALF_PI, p),
                NativeGate::Cz(a, b) => format!("cz q[{}],q[{}];\n", idx(a), idx(b)),
                // exp(−iπθ·Z⊗Z) = CX·(I⊗Rz(2πθ))·CX.
                NativeGate::Zz(a, b, t) => {
                    let (a, b) = (idx(a), idx(b));
                    format!("cx q[{a}],q[{b}];\nrz({}) q[{b}];\ncx q[{a}],q[{b}];\n", t.qasm())
                }
            });
        }
        for &(q, a) in &self.residual {
            text.push_str(&format!("rz({}) q[{}];\n", a.qasm(), idx(q)));
        }
        let program = super::from_qasm(&text).map_err(|e| NativeError::Coupling(format!("re-reading the program: {e}")))?;
        let mut reference = Circuit::new(m as u8);
        let at = |l: u8| idx(self.final_layout[usize::from(l)]) as u8;
        for g in &c.ops {
            reference.ops.push(Gate {
                base: g.base,
                controls: g.controls.iter().map(|&q| at(q)).collect(),
                target: at(g.target),
                param: if g.base == BaseGate::P { g.param } else { 0 },
            });
        }
        let a = program.simulate().map_err(|e| NativeError::Device(DeviceError::Invalid(e)))?;
        let b = reference.simulate().map_err(|e| NativeError::Device(DeviceError::Invalid(e)))?;
        let fidelity = a.fidelity_fx(&b);
        // Fixed-point rounding over a long program: a few parts per million.
        if (fidelity - ONE).abs() > 1 << 14 {
            return Err(NativeError::Mismatch { fidelity });
        }
        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn gate(base: BaseGate, controls: &[u8], target: u8, param: u16) -> Gate {
        Gate { base, controls: controls.to_vec(), target, param }
    }

    fn prepared(n: u8) -> Circuit {
        let mut c = Circuit::new(n);
        for q in 0..n {
            c.ops.push(gate(BaseGate::H, &[], q, 0));
            c.ops.push(gate(BaseGate::P, &[], q, u16::from(q) + 2));
        }
        c
    }

    fn random_circuit(n: u8, gates: usize, seed: u64) -> Circuit {
        let mut s = seed;
        let mut next = || {
            s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            (s >> 33) as u32
        };
        let bases = [BaseGate::X, BaseGate::Y, BaseGate::Z, BaseGate::H, BaseGate::S, BaseGate::Sdg, BaseGate::T, BaseGate::P];
        let mut c = prepared(n);
        for _ in 0..gates {
            let target = (next() % u32::from(n)) as u8;
            let mut controls = Vec::new();
            for _ in 0..next() % 3 {
                let q = (next() % u32::from(n)) as u8;
                if q != target && !controls.contains(&q) {
                    controls.push(q);
                }
            }
            let base = bases[(next() as usize) % bases.len()];
            let param = if base == BaseGate::P { 1 + (next() % 10) as u16 } else { 0 };
            c.ops.push(gate(base, &controls, target, param));
        }
        c
    }

    #[test]
    fn dyadic_angles_reduce_add_and_print() {
        let d = |m, k| Dyadic::new(m, k).unwrap();
        assert_eq!(d(2, 2), d(1, 1));
        assert_eq!(d(4, 2), Dyadic::ZERO);
        assert_eq!(d(1, 2) + d(1, 2), d(1, 1));
        assert_eq!(d(1, 3) + d(7, 3), Dyadic::ZERO);
        assert_eq!(-d(3, 3), d(5, 3));
        assert_eq!(d(1, 1).half(), Some(d(1, 2)));
        assert_eq!(d(1, 32).half(), None);
        assert_eq!(Dyadic::new(1, 33), None);
        assert_eq!(d(0, 0).qasm(), "0");
        assert_eq!(d(1, 1).qasm(), "pi");
        assert_eq!(d(1, 3).qasm(), "pi/4");
        assert_eq!(d(3, 2).qasm(), "-pi/2");
        assert_eq!(d(13, 4).qasm(), "-3*pi/8");
        assert_eq!(d(3, 4).qasm(), "3*pi/8");
        assert_eq!(d(1, 3).turns(), (1, 3));
    }

    #[test]
    fn every_gate_compiles_exactly_in_both_bases() {
        let bases = [
            (BaseGate::X, 0),
            (BaseGate::Y, 0),
            (BaseGate::Z, 0),
            (BaseGate::H, 0),
            (BaseGate::S, 0),
            (BaseGate::Sdg, 0),
            (BaseGate::T, 0),
            (BaseGate::Tdg, 0),
            (BaseGate::P, 1),
            (BaseGate::P, 5),
        ];
        let map = Coupling::grid(2, 3);
        for (base, param) in bases {
            for k in 0..=3u8 {
                let controls: Vec<u8> = [3u8, 0, 4].iter().copied().take(usize::from(k)).collect();
                let mut c = prepared(5);
                c.ops.push(gate(base, &controls, 1, param));
                for basis in [Basis::RzSxXCz, Basis::PrxCz, Basis::GpiZz] {
                    let p = compile(&c, basis, &map).unwrap();
                    p.verify(&c, &map).unwrap_or_else(|e| panic!("{base:?} {controls:?} {basis:?}: {e}"));
                }
            }
        }
    }

    #[test]
    fn random_circuits_route_onto_a_line_and_a_grid() {
        for seed in 0..12u64 {
            let c = random_circuit(5, 16, seed);
            for map in [Coupling::line(7), Coupling::grid(3, 3)] {
                for basis in [Basis::RzSxXCz, Basis::PrxCz, Basis::GpiZz] {
                    let p = compile(&c, basis, &map).unwrap();
                    p.verify(&c, &map).unwrap_or_else(|e| panic!("seed {seed} {basis:?}: {e}"));
                    for g in &p.gates {
                        if let NativeGate::Cz(a, b) = *g {
                            assert!(map.connected(a, b));
                        }
                    }
                }
            }
        }
    }

    #[test]
    fn a_line_needs_swaps_and_the_program_records_where_qubits_end() {
        // CZ between the ends of a 4-qubit line chain: at least one SWAP.
        let mut c = Circuit::new(4);
        c.ops.push(gate(BaseGate::H, &[], 0, 0));
        c.ops.push(gate(BaseGate::X, &[0], 1, 0));
        c.ops.push(gate(BaseGate::X, &[1], 2, 0));
        c.ops.push(gate(BaseGate::X, &[2], 3, 0));
        c.ops.push(gate(BaseGate::Z, &[0], 3, 0));
        let map = Coupling::line(4);
        let p = compile(&c, Basis::RzSxXCz, &map).unwrap();
        assert!(p.swaps >= 1);
        assert_ne!(p.layout, p.final_layout);
        p.verify(&c, &map).unwrap();
    }

    #[test]
    fn compilation_is_deterministic() {
        let c = random_circuit(4, 20, 99);
        let map = Coupling::grid(3, 3);
        assert_eq!(compile(&c, Basis::PrxCz, &map), compile(&c, Basis::PrxCz, &map));
    }

    #[test]
    fn a_wrong_program_fails_verification() {
        let c = random_circuit(4, 12, 7);
        let map = Coupling::grid(2, 3);
        for basis in [Basis::RzSxXCz, Basis::PrxCz, Basis::GpiZz] {
            let good = compile(&c, basis, &map).unwrap();
            good.verify(&c, &map).unwrap();
            // A dropped gate.
            let mut p = good.clone();
            let i = p.gates.iter().position(|g| !matches!(g, NativeGate::Rz(..))).unwrap();
            p.gates.remove(i);
            assert!(matches!(p.verify(&c, &map), Err(NativeError::Mismatch { .. })), "{basis:?}");
            // A cz on a pair the map does not connect.
            let mut p = good.clone();
            p.gates.push(NativeGate::Cz(0, 5));
            assert_eq!(p.verify(&c, &map), Err(NativeError::NotNative { gate: p.gates.len() - 1 }));
            // A wrong final layout.
            let mut p = good.clone();
            p.final_layout.swap(0, 1);
            assert!(p.verify(&c, &map).is_err(), "{basis:?}");
        }
        // A gate from the other basis.
        let mut p = compile(&c, Basis::PrxCz, &map).unwrap();
        p.gates.push(NativeGate::Sx(0));
        assert!(matches!(p.verify(&c, &map), Err(NativeError::NotNative { .. })));
    }

    #[test]
    fn a_map_too_small_or_split_is_refused() {
        let c = prepared(4);
        let split = Coupling::new(6, &[(0, 1), (1, 2), (3, 4), (4, 5)]).unwrap();
        assert_eq!(compile(&c, Basis::RzSxXCz, &split), Err(NativeError::TooFewQubits { need: 4, largest: 3 }));
        assert!(Coupling::new(3, &[(0, 0)]).is_err());
        assert!(Coupling::new(3, &[(0, 3)]).is_err());
    }

    #[test]
    fn the_qasm_names_physical_qubits_and_measures_where_each_logical_qubit_ends() {
        let mut c = Circuit::new(3);
        c.ops.push(gate(BaseGate::H, &[], 0, 0));
        c.ops.push(gate(BaseGate::X, &[0], 2, 0));
        let map = Coupling::line(5);
        let p = compile(&c, Basis::RzSxXCz, &map).unwrap();
        let text = p.qasm3();
        assert!(text.starts_with("OPENQASM 3.0;\ninclude \"stdgates.inc\";\nbit[3] c;\n"));
        for (i, q) in p.final_layout.iter().enumerate() {
            assert!(text.contains(&format!("c[{i}] = measure ${q};")), "{text}");
        }
        for line in text.lines().skip(3) {
            let name: String = line.chars().take_while(|c| c.is_ascii_alphanumeric()).collect();
            assert!(["rz", "sx", "x", "cz", "c"].contains(&name.as_str()), "{line}");
        }
        let prx = compile(&c, Basis::PrxCz, &map).unwrap().qasm3();
        assert!(prx.contains("gate prx(theta, phi) a"));
        let gpi = compile(&c, Basis::GpiZz, &map).unwrap().qasm3();
        assert!(gpi.contains("gate zz(theta) a, b") && gpi.contains("zz(0.25) $"), "{gpi}");
    }

    #[test]
    fn turns_print_exactly_in_decimal() {
        let d = |m, k| Dyadic::new(m, k).unwrap();
        assert_eq!(turns(Dyadic::ZERO), "0");
        assert_eq!(turns(d(1, 1)), "0.5");
        assert_eq!(turns(d(1, 2)), "0.25");
        assert_eq!(turns(d(3, 2)), "0.75");
        assert_eq!(turns(d(5, 4)), "0.3125");
        assert_eq!(turns(d(1, 32)), "0.00000000023283064365386962890625");
        // Read back as an f64, each is exact: m/2^k has at most 53 significant bits.
        for (m, k) in [(1u64, 1u8), (3, 2), (5, 4), (1, 32), (4294967295, 32)] {
            assert_eq!(turns(d(m, k)).parse::<f64>().unwrap(), m as f64 / (1u64 << k) as f64);
        }
    }
}