wai-quantum 0.5.0

//! Device-ready OpenQASM: a circuit in the form a quantum processor's service
//! accepts, and the service's counts read back onto this crate's basis.
//!
//! [`super::to_qasm`] writes the crate's own gate set, which is wider than any
//! device's: a gate may carry any number of controls. A service reads only its
//! dialect's standard library (OpenQASM 2's `qelib1.inc`, OpenQASM 3's
//! `stdgates.inc`), and needs every qubit measured into a classical register.
//! [`device_qasm`] lowers the circuit to those standard gates and measures each
//! qubit `i` into bit `i` of a register `c`.
//!
//! # Exact or refused
//!
//! Every lowering is an identity, not an approximation, and every angle it
//! writes is a signed dyadic turn `±2π/2^k`, so the output reads back through
//! [`super::from_qasm`] without rounding:
//!
//! - a controlled `Z`, `S`, `T` or `P` is a controlled phase; with two or more
//!   controls it splits, by the identity `θ·ab = ½θ·(a + b − (a⊕b))`, into
//!   controlled phases of half the angle around a multi-controlled `X` on the last
//!   control;
//! - a multi-controlled `X` is a multi-controlled `Z` between two `H`s (two
//!   controls stay `ccx`);
//! - a controlled `Y` is a controlled `X` between `S†` and `S`;
//! - a multi-controlled `H` is a multi-controlled `Z` between `ry(−π/4)` and
//!   `ry(π/4)`, since `Ry(π/4)·Z·Ry(−π/4) = H`.
//!
//! Each halving costs one bit of angle, so a gate whose lowering would need a turn
//! finer than `2π/2^DYADIC_MAX` is refused, as `from_qasm` refuses such an angle.
//! The lowered length grows about threefold per control, so a gate with more than
//! [`MAX_LOWERED_CONTROLS`] controls is refused too.
//!
//! The angles written are the ideal ones whatever the circuit's gateset, as in
//! [`super::to_qasm`].
//!
//! # Bit order
//!
//! Qubit `i` is bit `i` of this crate's basis index, and is measured into
//! classical bit `c[i]`. A service that reports a sample as the integer whose bit
//! `i` is `c[i]` therefore reports this crate's index unchanged, and one that
//! prints the register with `c[0]` last prints it in binary.
//! [`counts_from_bitstrings`] and [`counts_from_integers`] read either back.

use crate::quantum::{BaseGate, Circuit, QuantumError, DYADIC_MAX};

/// The most controls a gate may carry and still be lowered: each control
/// roughly triples the lowered length (a gate with 8 controls lowers to some
/// thousands of standard gates).
pub const MAX_LOWERED_CONTROLS: usize = 8;

/// An OpenQASM dialect.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Dialect {
    /// OpenQASM 2.0 with `qelib1.inc`: phases are `u1` and `cu1`.
    Qasm2,
    /// OpenQASM 3.0 with `stdgates.inc`: phases are `p` and `cp`.
    Qasm3,
}

/// A signed dyadic turn `±2π/2^k`, `1 ≤ k ≤ DYADIC_MAX`.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub struct Turn {
    pub k: u8,
    pub negative: bool,
}

impl Turn {
    /// `+2π/2^k`.
    pub const fn of(k: u8) -> Self {
        Turn { k, negative: false }
    }

    /// The turn negated: `∓2π/2^k`.
    pub const fn neg(self) -> Self {
        Turn { k: self.k, negative: !self.negative }
    }

    /// Half the turn, or `None` past `2π/2^DYADIC_MAX`.
    pub fn half(self) -> Option<Self> {
        (usize::from(self.k) < DYADIC_MAX).then_some(Turn { k: self.k + 1, negative: self.negative })
    }

    /// The angle as OpenQASM writes it: `2π/2^k = π/2^(k−1)`, so `pi`, `pi/2`,
    /// `-pi/4`.
    pub fn qasm(self) -> String {
        let sign = if self.negative { "-" } else { "" };
        if self.k == 1 {
            format!("{sign}pi")
        } else {
            format!("{sign}pi/{}", 1u64 << (self.k - 1))
        }
    }
}

/// A gate both standard libraries define, on qubit indices.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum StdGate {
    X(u8),
    Y(u8),
    Z(u8),
    H(u8),
    S(u8),
    Sdg(u8),
    T(u8),
    Tdg(u8),
    /// `p(θ)` (OpenQASM 2: `u1`).
    P(u8, Turn),
    /// `ry(θ)`.
    Ry(u8, Turn),
    /// Control, target.
    Cx(u8, u8),
    Cy(u8, u8),
    Cz(u8, u8),
    Ch(u8, u8),
    /// `cp(θ)` (OpenQASM 2: `cu1`): control, target.
    Cp(u8, u8, Turn),
    /// Two controls, target.
    Ccx(u8, u8, u8),
}

impl StdGate {
    /// The gate as one OpenQASM statement in `d`.
    pub fn qasm(&self, d: Dialect) -> String {
        let (p, cp) = match d {
            Dialect::Qasm2 => ("u1", "cu1"),
            Dialect::Qasm3 => ("p", "cp"),
        };
        match *self {
            StdGate::X(q) => format!("x q[{q}];"),
            StdGate::Y(q) => format!("y q[{q}];"),
            StdGate::Z(q) => format!("z q[{q}];"),
            StdGate::H(q) => format!("h q[{q}];"),
            StdGate::S(q) => format!("s q[{q}];"),
            StdGate::Sdg(q) => format!("sdg q[{q}];"),
            StdGate::T(q) => format!("t q[{q}];"),
            StdGate::Tdg(q) => format!("tdg q[{q}];"),
            StdGate::P(q, t) => format!("{p}({}) q[{q}];", t.qasm()),
            StdGate::Ry(q, t) => format!("ry({}) q[{q}];", t.qasm()),
            StdGate::Cx(c, t) => format!("cx q[{c}],q[{t}];"),
            StdGate::Cy(c, t) => format!("cy q[{c}],q[{t}];"),
            StdGate::Cz(c, t) => format!("cz q[{c}],q[{t}];"),
            StdGate::Ch(c, t) => format!("ch q[{c}],q[{t}];"),
            StdGate::Cp(c, t, a) => format!("{cp}({}) q[{c}],q[{t}];", a.qasm()),
            StdGate::Ccx(a, b, t) => format!("ccx q[{a}],q[{b}],q[{t}];"),
        }
    }
}

/// Why a circuit cannot be exported, or counts cannot be read back.
#[derive(Debug, PartialEq, Eq)]
pub enum DeviceError {
    /// The circuit fails [`Circuit::validate`].
    Invalid(QuantumError),
    /// Op `op` has more than [`MAX_LOWERED_CONTROLS`] controls.
    TooManyControls { op: usize, controls: usize },
    /// Lowering op `op` would need a turn finer than `2π/2^DYADIC_MAX`.
    TooFine { op: usize },
    /// A reported outcome is not a basis state of the register.
    Outcome(String),
}

impl core::fmt::Display for DeviceError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            DeviceError::Invalid(e) => write!(f, "invalid circuit: {e:?}"),
            DeviceError::TooManyControls { op, controls } => {
                write!(f, "op {op} has {controls} controls; at most {MAX_LOWERED_CONTROLS} are lowered")
            }
            DeviceError::TooFine { op } => {
                write!(f, "op {op} would need a turn finer than 2π/2^{DYADIC_MAX} to lower exactly")
            }
            DeviceError::Outcome(s) => write!(f, "outcome {s}"),
        }
    }
}

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

/// Lower a circuit to gates both standard libraries define. Exact: see the
/// module documentation.
pub fn lower(c: &Circuit) -> Result<Vec<StdGate>, DeviceError> {
    c.validate().map_err(DeviceError::Invalid)?;
    let mut out = Vec::new();
    for (op, g) in c.ops.iter().enumerate() {
        if g.controls.len() > MAX_LOWERED_CONTROLS {
            return Err(DeviceError::TooManyControls { op, controls: g.controls.len() });
        }
        lower_gate(&mut out, g.base, g.param, &g.controls, g.target).ok_or(DeviceError::TooFine { op })?;
    }
    Ok(out)
}

/// The phase a diagonal base gate applies to `|1⟩`, as a turn.
fn phase_of(base: BaseGate, param: u16) -> Option<Turn> {
    match base {
        BaseGate::Z => Some(Turn::of(1)),
        BaseGate::S => Some(Turn::of(2)),
        BaseGate::Sdg => Some(Turn::of(2).neg()),
        BaseGate::T => Some(Turn::of(3)),
        BaseGate::Tdg => Some(Turn::of(3).neg()),
        // `validate` holds `param` to `1..=DYADIC_MAX`, which fits a `u8`.
        BaseGate::P => Some(Turn::of(param as u8)),
        _ => None,
    }
}

fn lower_gate(out: &mut Vec<StdGate>, base: BaseGate, param: u16, cs: &[u8], t: u8) -> Option<()> {
    match (base, cs) {
        (BaseGate::I, _) => {}
        (BaseGate::X, []) => out.push(StdGate::X(t)),
        (BaseGate::Y, []) => out.push(StdGate::Y(t)),
        (BaseGate::Z, []) => out.push(StdGate::Z(t)),
        (BaseGate::H, []) => out.push(StdGate::H(t)),
        (BaseGate::S, []) => out.push(StdGate::S(t)),
        (BaseGate::Sdg, []) => out.push(StdGate::Sdg(t)),
        (BaseGate::T, []) => out.push(StdGate::T(t)),
        (BaseGate::Tdg, []) => out.push(StdGate::Tdg(t)),
        (BaseGate::P, []) => out.push(StdGate::P(t, Turn::of(param as u8))),
        (BaseGate::X, _) => mcx(out, cs, t)?,
        (BaseGate::Y, _) => {
            out.push(StdGate::Sdg(t));
            mcx(out, cs, t)?;
            out.push(StdGate::S(t));
        }
        (BaseGate::Z, [c]) => out.push(StdGate::Cz(*c, t)),
        (BaseGate::Z, [_, _]) => {
            out.push(StdGate::H(t));
            mcx(out, cs, t)?;
            out.push(StdGate::H(t));
        }
        (BaseGate::H, [c]) => out.push(StdGate::Ch(*c, t)),
        (BaseGate::H, _) => {
            out.push(StdGate::Ry(t, Turn::of(3).neg()));
            mcphase(out, cs, t, Turn::of(1))?;
            out.push(StdGate::Ry(t, Turn::of(3)));
        }
        (base, _) => mcphase(out, cs, t, phase_of(base, param)?)?,
    }
    Some(())
}

/// `X` on `t` when every qubit in `cs` (at least one) is set.
fn mcx(out: &mut Vec<StdGate>, cs: &[u8], t: u8) -> Option<()> {
    match cs {
        [] => out.push(StdGate::X(t)),
        [c] => out.push(StdGate::Cx(*c, t)),
        [a, b] => out.push(StdGate::Ccx(*a, *b, t)),
        _ => {
            out.push(StdGate::H(t));
            mcphase(out, cs, t, Turn::of(1))?;
            out.push(StdGate::H(t));
        }
    }
    Some(())
}

/// The phase `turn` on `|1⟩` of `t` when every qubit in `cs` is set: for
/// `cs = rest ++ [last]`, `θ·r·l = ½θ·(l + r − (r⊕l))` with `r` the conjunction
/// of `rest`, the `r⊕l` written by a multi-controlled `X` onto `last`.
fn mcphase(out: &mut Vec<StdGate>, cs: &[u8], t: u8, turn: Turn) -> Option<()> {
    match cs {
        [] => out.push(StdGate::P(t, turn)),
        [c] => out.push(StdGate::Cp(*c, t, turn)),
        [rest @ .., last] => {
            let half = turn.half()?;
            out.push(StdGate::Cp(*last, t, half));
            mcx(out, rest, *last)?;
            out.push(StdGate::Cp(*last, t, half.neg()));
            mcx(out, rest, *last)?;
            mcphase(out, rest, t, half)?;
        }
    }
    Some(())
}

/// The circuit as OpenQASM a device service accepts: `d`'s header and standard
/// library, a quantum register `q` and a classical register `c` of the circuit's
/// width, the lowered gates, and qubit `i` measured into `c[i]`.
pub fn device_qasm(c: &Circuit, d: Dialect) -> Result<String, DeviceError> {
    let gates = lower(c)?;
    let n = c.n_qubits;
    let mut s = match d {
        Dialect::Qasm2 => format!("OPENQASM 2.0;\ninclude \"qelib1.inc\";\nqreg q[{n}];\ncreg c[{n}];\n"),
        Dialect::Qasm3 => format!("OPENQASM 3.0;\ninclude \"stdgates.inc\";\nqubit[{n}] q;\nbit[{n}] c;\n"),
    };
    for g in &gates {
        s.push_str(&g.qasm(d));
        s.push('\n');
    }
    for i in 0..n {
        s.push_str(&match d {
            Dialect::Qasm2 => format!("measure q[{i}] -> c[{i}];\n"),
            Dialect::Qasm3 => format!("c[{i}] = measure q[{i}];\n"),
        });
    }
    Ok(s)
}

/// How a service prints a measured register.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum BitOrder {
    /// `c[n−1] … c[0]`: the register in binary, `c[0]` last.
    Bit0Last,
    /// `c[0] … c[n−1]`: `c[0]` first.
    Bit0First,
}

/// Counts over this crate's basis index (`2^n` entries) from a service's counts
/// keyed by printed register, each exactly `n` characters of `0` and `1`.
/// A repeated outcome adds.
pub fn counts_from_bitstrings<'a>(
    n: u8,
    counts: impl IntoIterator<Item = (&'a str, u64)>,
    order: BitOrder,
) -> Result<Vec<u64>, DeviceError> {
    let width = register_width(n)?;
    let mut out = vec![0u64; 1 << width];
    for (s, k) in counts {
        if s.len() != width || !s.bytes().all(|b| b == b'0' || b == b'1') {
            return Err(DeviceError::Outcome(format!("{s:?} is not {width} binary digits")));
        }
        let mut x = 0usize;
        for (pos, b) in s.bytes().enumerate() {
            let bit = match order {
                BitOrder::Bit0Last => width - 1 - pos,
                BitOrder::Bit0First => pos,
            };
            x |= usize::from(b - b'0') << bit;
        }
        add(&mut out, x, k)?;
    }
    Ok(out)
}

/// Counts over this crate's basis index from a service's counts keyed by the
/// integer whose bit `i` is `c[i]`, which is this crate's index itself.
pub fn counts_from_integers(n: u8, counts: impl IntoIterator<Item = (u64, u64)>) -> Result<Vec<u64>, DeviceError> {
    let width = register_width(n)?;
    let mut out = vec![0u64; 1 << width];
    for (v, k) in counts {
        let x = usize::try_from(v)
            .ok()
            .filter(|&x| x < out.len())
            .ok_or_else(|| DeviceError::Outcome(format!("{v} is outside a {width}-bit register")))?;
        add(&mut out, x, k)?;
    }
    Ok(out)
}

fn register_width(n: u8) -> Result<usize, DeviceError> {
    let max = crate::quantum::MAX_QUBITS;
    if n == 0 || n > max {
        return Err(DeviceError::Outcome(format!("a {n}-qubit register; counts are read for 1..={max}")));
    }
    Ok(usize::from(n))
}

fn add(out: &mut [u64], x: usize, k: u64) -> Result<(), DeviceError> {
    out[x] = out[x].checked_add(k).ok_or_else(|| DeviceError::Outcome(format!("the count of {x} overflows")))?;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum::{Gate, ONE};
    use crate::quantum_toolchain::from_qasm;

    /// A circuit whose state, read back from `text`, matches `c`'s up to global
    /// phase, at the fixed-point simulator's resolution.
    fn same_state(c: &Circuit, text: &str) {
        let back = from_qasm(text).unwrap_or_else(|e| panic!("reparse: {e}\n{text}"));
        let f = back.simulate().unwrap().fidelity_fx(&c.simulate().unwrap());
        assert!((f - ONE).abs() < 1 << 12, "fidelity {f}, want {ONE}\n{text}");
    }

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

    /// Superpose every qubit first, with distinct phases, so a wrong phase or a
    /// wrong control condition anywhere shows in the state.
    fn prepared(n: u8) -> Circuit {
        let mut c = Circuit::new(n);
        for q in 0..n {
            c.ops.push(one(BaseGate::H, &[], q, 0));
            c.ops.push(one(BaseGate::P, &[], q, u16::from(q) + 2));
        }
        c
    }

    #[test]
    fn every_base_gate_under_zero_to_four_controls_lowers_exactly() {
        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, 4),
            (BaseGate::P, 9),
        ];
        for (base, param) in bases {
            for k in 0..=4u8 {
                // Controls out of order and on both sides of the target.
                let controls: Vec<u8> = [4u8, 0, 3, 1].iter().copied().take(usize::from(k)).collect();
                let mut c = prepared(5);
                c.ops.push(one(base, &controls, 2, param));
                for d in [Dialect::Qasm2, Dialect::Qasm3] {
                    let text = device_qasm(&c, d).unwrap();
                    assert!(!text.contains("ccz") && !text.contains("unmapped"), "{text}");
                    same_state(&c, &text);
                }
            }
        }
    }

    #[test]
    fn random_circuits_round_trip_in_both_dialects() {
        let mut s = 0x5eed_u64;
        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::T, BaseGate::Tdg, BaseGate::P];
        for _ in 0..24 {
            let n = 6u8;
            let mut c = prepared(n);
            for _ in 0..14 {
                let target = (next() % u32::from(n)) as u8;
                let mut controls = Vec::new();
                for _ in 0..next() % 4 {
                    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() % 12) as u16 } else { 0 };
                c.ops.push(one(base, &controls, target, param));
            }
            for d in [Dialect::Qasm2, Dialect::Qasm3] {
                same_state(&c, &device_qasm(&c, d).unwrap());
            }
        }
    }

    #[test]
    fn the_output_names_only_its_dialects_standard_gates_and_measures_every_qubit() {
        let mut c = prepared(4);
        c.ops.push(one(BaseGate::Z, &[0, 1, 3], 2, 0));
        c.ops.push(one(BaseGate::H, &[0, 1], 3, 0));
        c.ops.push(one(BaseGate::P, &[2], 0, 5));
        let v2: &[&str] = &["x", "y", "z", "h", "s", "sdg", "t", "tdg", "u1", "ry", "cx", "cy", "cz", "ch", "cu1", "ccx"];
        let v3: &[&str] = &["x", "y", "z", "h", "s", "sdg", "t", "tdg", "p", "ry", "cx", "cy", "cz", "ch", "cp", "ccx"];
        for (d, names, header) in [(Dialect::Qasm2, v2, "OPENQASM 2.0;"), (Dialect::Qasm3, v3, "OPENQASM 3.0;")] {
            let text = device_qasm(&c, d).unwrap();
            let mut lines = text.lines();
            assert_eq!(lines.next(), Some(header));
            let mut measured = Vec::new();
            for line in lines.skip(3) {
                if line.contains("measure") {
                    measured.push(line.to_string());
                    continue;
                }
                let name: String = line.chars().take_while(|c| c.is_ascii_alphanumeric()).collect();
                assert!(names.contains(&name.as_str()), "{d:?}: {line}");
            }
            let want: Vec<String> = (0..4)
                .map(|i| match d {
                    Dialect::Qasm2 => format!("measure q[{i}] -> c[{i}];"),
                    Dialect::Qasm3 => format!("c[{i}] = measure q[{i}];"),
                })
                .collect();
            assert_eq!(measured, want);
        }
    }

    #[test]
    fn turns_print_as_multiples_of_pi() {
        assert_eq!(Turn::of(1).qasm(), "pi");
        assert_eq!(Turn::of(2).qasm(), "pi/2");
        assert_eq!(Turn::of(3).neg().qasm(), "-pi/4");
        assert_eq!(Turn::of(32).qasm(), "pi/2147483648");
        assert_eq!(Turn::of(32).half(), None);
    }

    #[test]
    fn what_cannot_be_lowered_exactly_is_refused() {
        // Two controls halve the turn once: P(32) has no half in the table.
        let mut c = Circuit::new(3);
        c.ops.push(one(BaseGate::P, &[0, 1], 2, 32));
        assert_eq!(lower(&c), Err(DeviceError::TooFine { op: 0 }));
        // One control needs no halving.
        let mut c = Circuit::new(3);
        c.ops.push(one(BaseGate::P, &[0], 2, 32));
        assert!(lower(&c).is_ok());
        let mut c = Circuit::new(10);
        c.ops.push(one(BaseGate::X, &[0, 1, 2, 3, 4, 5, 6, 7, 8], 9, 0));
        assert_eq!(lower(&c), Err(DeviceError::TooManyControls { op: 0, controls: 9 }));
        let mut c = Circuit::new(2);
        c.ops.push(one(BaseGate::P, &[], 0, 0));
        assert!(matches!(lower(&c), Err(DeviceError::Invalid(_))));
    }

    #[test]
    fn a_bell_pair_lowers_to_the_textbook_program() {
        let text = device_qasm(&crate::quantum_toolchain::bell(), Dialect::Qasm2).unwrap();
        assert_eq!(
            text,
            "OPENQASM 2.0;\ninclude \"qelib1.inc\";\nqreg q[2];\ncreg c[2];\nh q[0];\ncx q[0],q[1];\nmeasure q[0] -> c[0];\nmeasure q[1] -> c[1];\n"
        );
    }

    #[test]
    fn counts_read_back_onto_the_crates_index() {
        // X on qubit 0 only: index 1 in this crate; a service prints `c[0]` last as
        // "001" and reports the integer 1.
        let mut c = Circuit::new(3);
        c.ops.push(one(BaseGate::X, &[], 0, 0));
        let shots = c.simulate().unwrap().sample_shots(1, 13);
        assert_eq!(shots[1], 13);
        let a = counts_from_bitstrings(3, [("001", 10), ("110", 2), ("001", 1)], BitOrder::Bit0Last).unwrap();
        let b = counts_from_bitstrings(3, [("100", 10), ("011", 2), ("100", 1)], BitOrder::Bit0First).unwrap();
        let i = counts_from_integers(3, [(1, 11), (6, 2)]).unwrap();
        assert_eq!(a, vec![0, 11, 0, 0, 0, 0, 2, 0]);
        assert_eq!(a, b);
        assert_eq!(a, i);
        assert!(counts_from_bitstrings(3, [("01", 1)], BitOrder::Bit0Last).is_err());
        assert!(counts_from_bitstrings(3, [("0 01", 1)], BitOrder::Bit0Last).is_err());
        assert!(counts_from_bitstrings(2, [("12", 1)], BitOrder::Bit0Last).is_err());
        assert!(counts_from_integers(3, [(8, 1)]).is_err());
        assert!(counts_from_integers(3, [(1, u64::MAX), (1, 1)]).is_err());
        assert!(counts_from_integers(0, []).is_err());
    }
}