wai-quantum 0.3.18

//! What quantum mechanics actually says about encoding and transporting
//! information — including, most usefully, what it says you cannot do.
//!
//! It is tempting to attach "quantum" to a codec. The theory does not support it,
//! and the theorem that refuses is worth more than the claim would be:
//!
//! - **Holevo**: `n` qubits carry at most `n` classical bits. Quantum does not
//!   compress classical media. Whatever a quantum channel is for, it is not a
//!   smaller pipe for the same bits.
//! - **Schumacher**: a *quantum* source compresses to `S(ρ)` qubits per symbol and
//!   no further — the quantum analogue of Shannon, with von Neumann entropy in
//!   place of Shannon's.
//! - Where quantum genuinely wins is **protocols, not codecs**: superdense coding
//!   sends two classical bits down one qubit *given a shared pair*, teleportation
//!   moves a state that cannot be measured, and key distribution secures the
//!   channel ([`crate::quantum_comm`]).
//!
//! So this module computes the limits rather than promising to beat them. Given an
//! ensemble of states it will tell you the most classical information anyone could
//! ever extract, and given a source it will tell you how far it compresses.
//!
//! # Determinism
//!
//! Eigenvalues come from a Jacobi sweep using only IEEE-754 `+ − × ÷ √`, the same
//! discipline as [`crate::quantum_mps`]. Entropy needs a logarithm, and `ln` is a
//! library function that differs between platforms, so this computes `log₂` itself
//! by repeated squaring — multiply and compare, nothing else. The results are
//! therefore bit-reproducible everywhere, not merely close.

use crate::quantum::ONE;
use crate::quantum_info::DensityMatrix;

/// `log₂(x)` for `x > 0`, computed rather than looked up.
///
/// The exponent comes from the bit pattern; the mantissa's log is refined one bit
/// at a time by squaring. Only multiplication, comparison and halving are used, so
/// every machine produces the same bits — unlike `ln`, whose last places are a
/// property of the platform's libm.
fn log2_reproducible(x: f64) -> f64 {
    if x <= 0.0 {
        return f64::NEG_INFINITY;
    }
    let bits = x.to_bits();
    let exp = (((bits >> 52) & 0x7FF) as i32) - 1023;
    let mut m = f64::from_bits((bits & 0x000F_FFFF_FFFF_FFFF) | 0x3FF0_0000_0000_0000);
    let mut frac = 0.0f64;
    let mut weight = 0.5f64;
    for _ in 0..52 {
        m *= m;
        if m >= 2.0 {
            m *= 0.5;
            frac += weight;
        }
        weight *= 0.5;
    }
    exp as f64 + frac
}

/// Eigenvalues of a density matrix, ascending.
///
/// A complex Hermitian `H = A + iB` has the same spectrum as the real symmetric
/// `[[A, −B], [B, A]]`, with every eigenvalue doubled. Working there means a plain
/// real Jacobi sweep suffices, and real Jacobi needs only the four operations and a
/// square root.
pub fn eigenvalues(rho: &DensityMatrix) -> Vec<f64> {
    let d = rho.d;
    let n = 2 * d;
    let scale = ONE as f64;
    let mut m = vec![vec![0.0f64; n]; n];
    for r in 0..d {
        for c in 0..d {
            let e = rho.get(r, c);
            let (a, b) = (e.re as f64 / scale, e.im as f64 / scale);
            m[r][c] = a;
            m[r + d][c + d] = a;
            m[r][c + d] = -b;
            m[r + d][c] = b;
        }
    }
    // cyclic Jacobi: rotate away the largest off-diagonal entries until quiet
    for _ in 0..100 {
        let mut off = 0.0f64;
        for r in 0..n {
            for c in 0..n {
                if r != c {
                    off += m[r][c] * m[r][c];
                }
            }
        }
        if off <= 1e-30 {
            break;
        }
        for p in 0..n - 1 {
            for q in p + 1..n {
                if m[p][q] == 0.0 {
                    continue;
                }
                let theta = (m[q][q] - m[p][p]) / (2.0 * m[p][q]);
                let t = if theta >= 0.0 {
                    1.0 / (theta + (1.0 + theta * theta).sqrt())
                } else {
                    -1.0 / (-theta + (1.0 + theta * theta).sqrt())
                };
                let c = 1.0 / (1.0 + t * t).sqrt();
                let s = t * c;
                for k in 0..n {
                    let (mkp, mkq) = (m[k][p], m[k][q]);
                    m[k][p] = c * mkp - s * mkq;
                    m[k][q] = s * mkp + c * mkq;
                }
                for k in 0..n {
                    let (mpk, mqk) = (m[p][k], m[q][k]);
                    m[p][k] = c * mpk - s * mqk;
                    m[q][k] = s * mpk + c * mqk;
                }
            }
        }
    }
    let mut ev: Vec<f64> = (0..n).map(|i| m[i][i]).collect();
    ev.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    // each eigenvalue of H appears twice in the doubled matrix
    ev.into_iter().step_by(2).collect()
}

/// The **von Neumann entropy** `S(ρ) = −Σ λ log₂ λ`, in bits.
///
/// Zero for a pure state — it is fully compressible, since you already know it.
/// One bit for a maximally mixed qubit — nothing to remove. By Schumacher this is
/// exactly the number of qubits per symbol a quantum source compresses to.
pub fn von_neumann_entropy(rho: &DensityMatrix) -> f64 {
    let mut s = 0.0f64;
    for l in eigenvalues(rho) {
        if l > 1e-12 {
            s -= l * log2_reproducible(l);
        }
    }
    if s < 0.0 { 0.0 } else { s }
}

/// Qubits per symbol for a quantum source — Schumacher's limit, which is the
/// entropy.
pub fn schumacher_limit(rho: &DensityMatrix) -> f64 {
    von_neumann_entropy(rho)
}

/// The **Holevo bound** `χ = S(Σ pᵢ ρᵢ) − Σ pᵢ S(ρᵢ)`, in bits.
///
/// The ceiling on classical information anyone can extract from an ensemble of
/// quantum states, however clever the measurement. Two consequences worth stating
/// plainly: encoding classical data in non-orthogonal states loses information
/// permanently, and `χ ≤ log₂ d`, so a qubit never carries more than one bit.
///
/// Weights are fixed point (`ONE` = 1) and should sum to `ONE`.
pub fn holevo_bound(ensemble: &[(i64, DensityMatrix)]) -> Option<f64> {
    let avg = DensityMatrix::mixture(ensemble)?;
    let mut mixed = von_neumann_entropy(&avg);
    for (w, r) in ensemble {
        mixed -= (*w as f64 / ONE as f64) * von_neumann_entropy(r);
    }
    Some(if mixed < 0.0 { 0.0 } else { mixed })
}

/// Purity restated as a one-line sanity companion to the entropy: `Tr ρ²` as a
/// float, one for pure and `1/d` for maximally mixed.
pub fn purity(rho: &DensityMatrix) -> f64 {
    rho.purity() as f64 / ONE as f64
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum::Circuit;
    use crate::quantum_info::DensityMatrix;

    fn pure(c: &mut Circuit) -> DensityMatrix {
        DensityMatrix::from_pure(&c.simulate().unwrap())
    }

    #[test]
    fn log2_matches_the_library_without_using_it() {
        for x in [0.5f64, 0.25, 1.0, 0.75, 1.0 / 3.0, 0.1, 0.999, 1e-6] {
            let got = log2_reproducible(x);
            assert!((got - x.log2()).abs() < 1e-12, "log2({x}) = {got}, want {}", x.log2());
        }
    }

    /// A pure state is fully compressible: you already know it, so there is
    /// nothing to send.
    #[test]
    fn a_pure_state_has_no_entropy() {
        let mut c = Circuit::new(2);
        c.h(0).cx(0, 1);
        let s = von_neumann_entropy(&pure(&mut c));
        assert!(s.abs() < 1e-6, "pure state entropy {s}, want 0");
    }

    /// Half a Bell pair is maximally mixed: exactly one bit, and no compression
    /// scheme can do better. This is the entropy the state vector could not report.
    #[test]
    fn half_a_bell_pair_costs_exactly_one_qubit() {
        let mut c = Circuit::new(2);
        c.h(0).cx(0, 1);
        let red = pure(&mut c).partial_trace(&[0]);
        let s = von_neumann_entropy(&red);
        assert!((s - 1.0).abs() < 1e-4, "maximally mixed entropy {s}, want 1 bit");
        assert!((schumacher_limit(&red) - 1.0).abs() < 1e-4);
        assert!((purity(&red) - 0.5).abs() < 1e-3);
    }

    /// Orthogonal states are perfectly distinguishable, so a two-state ensemble
    /// carries its full bit.
    #[test]
    fn orthogonal_states_carry_a_full_bit() {
        let mut z = Circuit::new(1);
        let mut o = Circuit::new(1);
        o.x(0);
        let chi = holevo_bound(&[(ONE / 2, pure(&mut z)), (ONE / 2, pure(&mut o))]).unwrap();
        assert!((chi - 1.0).abs() < 1e-4, "χ = {chi}, want 1 bit");
    }

    /// The anti-hype result. Encode a bit in |0> versus |+> and the information is
    /// not merely hard to read — a fraction of it is gone, for every measurement
    /// anyone will ever build.
    #[test]
    fn non_orthogonal_encoding_loses_information_permanently() {
        let mut z = Circuit::new(1);
        let mut p = Circuit::new(1);
        p.h(0);
        let chi = holevo_bound(&[(ONE / 2, pure(&mut z)), (ONE / 2, pure(&mut p))]).unwrap();
        assert!(chi < 0.7, "χ = {chi} must fall short of a bit");
        assert!(chi > 0.5, "χ = {chi} should still be substantial");
    }

    /// Holevo's ceiling: a qubit never carries more than one classical bit, so
    /// quantum is not a smaller pipe for classical media.
    #[test]
    fn a_qubit_never_carries_more_than_one_bit() {
        let mut z = Circuit::new(1);
        let mut o = Circuit::new(1);
        o.x(0);
        let mut p = Circuit::new(1);
        p.h(0);
        let third = ONE / 3;
        let chi = holevo_bound(&[
            (third, pure(&mut z)),
            (third, pure(&mut o)),
            (ONE - 2 * third, pure(&mut p)),
        ])
        .unwrap();
        assert!(chi <= 1.0 + 1e-6, "χ = {chi} would beat Holevo, which is impossible");
    }

    #[test]
    fn entropy_is_reproducible() {
        let mut c = Circuit::new(3);
        c.h(0).cx(0, 1).cx(1, 2);
        let r = pure(&mut c).partial_trace(&[0, 1]);
        assert_eq!(von_neumann_entropy(&r).to_bits(), von_neumann_entropy(&r).to_bits());
    }
}