quantrs2_sim/

quantum_supremacy.rs

//! Quantum supremacy verification algorithms and benchmarks.
//!
//! This module implements various algorithms and statistical tests for verifying
//! quantum supremacy claims, including cross-entropy benchmarking, Porter-Thomas
//! distribution analysis, and linear cross-entropy benchmarking (Linear XEB).

use scirs2_core::ndarray::Array1;
use scirs2_core::random::prelude::*;
use scirs2_core::random::ChaCha8Rng;
use scirs2_core::random::{RngExt, SeedableRng};
use scirs2_core::Complex64;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::f64::consts::PI;

use crate::error::{Result, SimulatorError};
use crate::scirs2_integration::SciRS2Backend;

/// Cross-entropy benchmarking result
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct CrossEntropyResult {
    /// Linear cross-entropy benchmarking (XEB) fidelity
    pub linear_xeb_fidelity: f64,
    /// Cross-entropy difference
    pub cross_entropy_difference: f64,
    /// Number of samples used
    pub num_samples: usize,
    /// Statistical confidence (p-value)
    pub confidence: f64,
    /// Porter-Thomas test results
    pub porter_thomas: PorterThomasResult,
    /// Heavy output generation (HOG) score
    pub hog_score: f64,
    /// Computational cost comparison
    pub cost_comparison: CostComparison,
}

/// Porter-Thomas distribution analysis
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct PorterThomasResult {
    /// Chi-squared statistic
    pub chi_squared: f64,
    /// Degrees of freedom
    pub degrees_freedom: usize,
    /// P-value of the test
    pub p_value: f64,
    /// Whether distribution matches Porter-Thomas
    pub is_porter_thomas: bool,
    /// Kolmogorov-Smirnov test statistic
    pub ks_statistic: f64,
    /// Mean of the distribution
    pub mean: f64,
    /// Variance of the distribution
    pub variance: f64,
}

/// Heavy Output Generation analysis
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct HOGAnalysis {
    /// Fraction of heavy outputs
    pub heavy_fraction: f64,
    /// Expected heavy fraction for random circuit
    pub expected_heavy_fraction: f64,
    /// Threshold for heavy outputs
    pub threshold: f64,
    /// Total outputs analyzed
    pub total_outputs: usize,
    /// Heavy outputs count
    pub heavy_count: usize,
}

/// Computational cost comparison
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct CostComparison {
    /// Classical simulation time (seconds)
    pub classical_time: f64,
    /// Quantum execution time estimate (seconds)
    pub quantum_time: f64,
    /// Memory usage for classical simulation (bytes)
    pub classical_memory: usize,
    /// Speedup factor
    pub speedup_factor: f64,
    /// Number of operations
    pub operation_count: usize,
}

/// Quantum supremacy verifier
pub struct QuantumSupremacyVerifier {
    /// Number of qubits
    num_qubits: usize,
    /// Random number generator
    rng: ChaCha8Rng,
    /// `SciRS2` backend for optimization
    backend: Option<SciRS2Backend>,
    /// Verification parameters
    params: VerificationParams,
}

/// Parameters for quantum supremacy verification
#[derive(Debug, Clone)]
pub struct VerificationParams {
    /// Number of random circuits to test
    pub num_circuits: usize,
    /// Number of samples per circuit
    pub samples_per_circuit: usize,
    /// Circuit depth
    pub circuit_depth: usize,
    /// Gate set to use
    pub gate_set: GateSet,
    /// Significance level for statistical tests
    pub significance_level: f64,
    /// Porter-Thomas test bins
    pub pt_bins: usize,
    /// Heavy output threshold percentile
    pub heavy_threshold_percentile: f64,
}

/// Gate set for random circuit generation
#[derive(Debug, Clone)]
pub enum GateSet {
    /// Sycamore-like gate set (√X, √Y, √W, CZ)
    Sycamore,
    /// Google's quantum supremacy gate set
    GoogleSupremacy,
    /// Universal gate set (H, T, CNOT)
    Universal,
    /// IBM-like gate set (RZ, SX, CNOT)
    IBM,
    /// Custom gate set
    Custom(Vec<String>),
}

impl Default for VerificationParams {
    fn default() -> Self {
        Self {
            num_circuits: 100,
            samples_per_circuit: 1000,
            circuit_depth: 20,
            gate_set: GateSet::Sycamore,
            significance_level: 0.05,
            pt_bins: 100,
            heavy_threshold_percentile: 50.0,
        }
    }
}

impl QuantumSupremacyVerifier {
    /// Create new quantum supremacy verifier
    pub fn new(num_qubits: usize, params: VerificationParams) -> Result<Self> {
        Ok(Self {
            num_qubits,
            rng: ChaCha8Rng::from_rng(&mut thread_rng()),
            backend: None,
            params,
        })
    }

    /// Initialize with `SciRS2` backend
    pub fn with_scirs2_backend(mut self) -> Result<Self> {
        self.backend = Some(SciRS2Backend::new());
        Ok(self)
    }

    /// Set random seed for reproducibility
    #[must_use]
    pub fn with_seed(mut self, seed: u64) -> Self {
        self.rng = ChaCha8Rng::seed_from_u64(seed);
        self
    }

    /// Perform comprehensive quantum supremacy verification
    pub fn verify_quantum_supremacy(&mut self) -> Result<CrossEntropyResult> {
        let start_time = std::time::Instant::now();

        // Generate random quantum circuits
        let circuits = self.generate_random_circuits()?;

        // Compute ideal amplitudes (classical simulation)
        let ideal_amplitudes = self.compute_ideal_amplitudes(&circuits)?;

        // Simulate quantum sampling using ideal state amplitudes as the probability distribution
        let quantum_samples = self.simulate_quantum_sampling(&ideal_amplitudes)?;

        // Perform cross-entropy benchmarking
        let linear_xeb = self.compute_linear_xeb(&ideal_amplitudes, &quantum_samples)?;

        // Analyze Porter-Thomas distribution
        let porter_thomas = self.analyze_porter_thomas(&ideal_amplitudes)?;

        // Compute Heavy Output Generation score
        let hog_score = self.compute_hog_score(&ideal_amplitudes, &quantum_samples)?;

        // Cross-entropy difference
        let cross_entropy_diff =
            self.compute_cross_entropy_difference(&ideal_amplitudes, &quantum_samples)?;

        // Statistical confidence
        let confidence =
            self.compute_statistical_confidence(&ideal_amplitudes, &quantum_samples)?;

        let classical_time = start_time.elapsed().as_secs_f64();

        // Cost comparison
        let cost_comparison = CostComparison {
            classical_time,
            quantum_time: self.estimate_quantum_time()?,
            classical_memory: self.estimate_classical_memory(),
            speedup_factor: classical_time / self.estimate_quantum_time()?,
            operation_count: circuits.len() * self.params.circuit_depth,
        };

        Ok(CrossEntropyResult {
            linear_xeb_fidelity: linear_xeb,
            cross_entropy_difference: cross_entropy_diff,
            num_samples: self.params.samples_per_circuit,
            confidence,
            porter_thomas,
            hog_score,
            cost_comparison,
        })
    }

    /// Generate random quantum circuits
    fn generate_random_circuits(&mut self) -> Result<Vec<RandomCircuit>> {
        let mut circuits = Vec::with_capacity(self.params.num_circuits);

        for _ in 0..self.params.num_circuits {
            let circuit = self.generate_single_random_circuit()?;
            circuits.push(circuit);
        }

        Ok(circuits)
    }

    /// Generate a single random circuit
    fn generate_single_random_circuit(&mut self) -> Result<RandomCircuit> {
        let mut layers = Vec::with_capacity(self.params.circuit_depth);

        for layer_idx in 0..self.params.circuit_depth {
            let layer = match &self.params.gate_set {
                GateSet::Sycamore => self.generate_sycamore_layer(layer_idx)?,
                GateSet::GoogleSupremacy => self.generate_google_supremacy_layer(layer_idx)?,
                GateSet::Universal => self.generate_universal_layer(layer_idx)?,
                GateSet::IBM => self.generate_ibm_layer(layer_idx)?,
                GateSet::Custom(gates) => {
                    let gates_clone = gates.clone();
                    self.generate_custom_layer(layer_idx, &gates_clone)?
                }
            };
            layers.push(layer);
        }

        Ok(RandomCircuit {
            num_qubits: self.num_qubits,
            layers,
            measurement_pattern: self.generate_measurement_pattern(),
        })
    }

    /// Generate Sycamore-style layer
    fn generate_sycamore_layer(&mut self, layer_idx: usize) -> Result<CircuitLayer> {
        let mut gates = Vec::new();

        // Single-qubit gates
        for qubit in 0..self.num_qubits {
            let gate_type = match self.rng.random_range(0..3) {
                0 => "SqrtX",
                1 => "SqrtY",
                _ => "SqrtW",
            };

            gates.push(QuantumGate {
                gate_type: gate_type.to_string(),
                qubits: vec![qubit],
                parameters: vec![],
            });
        }

        // Two-qubit gates (CZ gates with connectivity pattern)
        let offset = layer_idx % 2;
        for row in 0..((self.num_qubits as f64).sqrt() as usize) {
            for col in (offset..((self.num_qubits as f64).sqrt() as usize)).step_by(2) {
                let qubit1 = row * ((self.num_qubits as f64).sqrt() as usize) + col;
                let qubit2 = qubit1 + 1;

                if qubit2 < self.num_qubits {
                    gates.push(QuantumGate {
                        gate_type: "CZ".to_string(),
                        qubits: vec![qubit1, qubit2],
                        parameters: vec![],
                    });
                }
            }
        }

        Ok(CircuitLayer { gates })
    }

    /// Generate Google quantum supremacy layer
    fn generate_google_supremacy_layer(&mut self, layer_idx: usize) -> Result<CircuitLayer> {
        // Similar to Sycamore but with specific gate rotations
        self.generate_sycamore_layer(layer_idx)
    }

    /// Generate universal gate set layer
    fn generate_universal_layer(&mut self, _layer_idx: usize) -> Result<CircuitLayer> {
        let mut gates = Vec::new();

        // Random single-qubit gates
        for qubit in 0..self.num_qubits {
            let gate_type = match self.rng.random_range(0..3) {
                0 => "H",
                1 => "T",
                _ => "S",
            };

            gates.push(QuantumGate {
                gate_type: gate_type.to_string(),
                qubits: vec![qubit],
                parameters: vec![],
            });
        }

        // Random CNOT gates
        let num_cnots = self.rng.random_range(1..=self.num_qubits / 2);
        for _ in 0..num_cnots {
            let control = self.rng.random_range(0..self.num_qubits);
            let target = self.rng.random_range(0..self.num_qubits);

            if control != target {
                gates.push(QuantumGate {
                    gate_type: "CNOT".to_string(),
                    qubits: vec![control, target],
                    parameters: vec![],
                });
            }
        }

        Ok(CircuitLayer { gates })
    }

    /// Generate IBM gate set layer
    fn generate_ibm_layer(&mut self, _layer_idx: usize) -> Result<CircuitLayer> {
        let mut gates = Vec::new();

        // RZ and SX gates
        for qubit in 0..self.num_qubits {
            // Random RZ rotation
            let angle = self.rng.random::<f64>() * 2.0 * PI;
            gates.push(QuantumGate {
                gate_type: "RZ".to_string(),
                qubits: vec![qubit],
                parameters: vec![angle],
            });

            // SX gate with probability 0.5
            if self.rng.random::<bool>() {
                gates.push(QuantumGate {
                    gate_type: "SX".to_string(),
                    qubits: vec![qubit],
                    parameters: vec![],
                });
            }
        }

        Ok(CircuitLayer { gates })
    }

    /// Generate custom gate set layer
    const fn generate_custom_layer(
        &self,
        _layer_idx: usize,
        _gates: &[String],
    ) -> Result<CircuitLayer> {
        // Implement custom gate generation
        Ok(CircuitLayer { gates: Vec::new() })
    }

    /// Generate measurement pattern
    fn generate_measurement_pattern(&self) -> Vec<usize> {
        // For now, measure all qubits in computational basis
        (0..self.num_qubits).collect()
    }

    /// Compute ideal amplitudes using classical state-vector simulation
    fn compute_ideal_amplitudes(
        &self,
        circuits: &[RandomCircuit],
    ) -> Result<Vec<Array1<Complex64>>> {
        let mut amplitudes = Vec::with_capacity(circuits.len());

        for circuit in circuits {
            let dim = 1usize << self.num_qubits;
            // Initialise to |0...0⟩
            let mut state = vec![Complex64::new(0.0, 0.0); dim];
            state[0] = Complex64::new(1.0, 0.0);

            // Apply circuit layers
            for layer in &circuit.layers {
                for gate in &layer.gates {
                    self.apply_gate_to_state(&mut state, gate)?;
                }
            }

            let final_state = Array1::from_vec(state);
            amplitudes.push(final_state);
        }

        Ok(amplitudes)
    }

    /// Apply a single gate to a raw state-vector slice.
    ///
    /// All gate matrices are given in the standard computational-basis order.
    fn apply_gate_to_state(&self, state: &mut Vec<Complex64>, gate: &QuantumGate) -> Result<()> {
        let n = self.num_qubits;
        let dim = state.len();

        /// Apply a 2×2 single-qubit gate matrix `m` (row-major, len 4) to `state`
        /// at `target_bit` (0 = least-significant qubit).
        fn apply_1q(state: &mut [Complex64], m: &[Complex64; 4], target_bit: usize, dim: usize) {
            for i in 0..dim {
                if (i >> target_bit) & 1 == 0 {
                    let j = i | (1 << target_bit);
                    let a = state[i];
                    let b = state[j];
                    state[i] = m[0] * a + m[1] * b;
                    state[j] = m[2] * a + m[3] * b;
                }
            }
        }

        /// Apply a 4×4 two-qubit gate matrix `m` (row-major, len 16) to `state`
        /// for (`ctrl_bit`, `tgt_bit`).  The matrix is ordered |00⟩,|01⟩,|10⟩,|11⟩
        /// where the first index is `ctrl_bit` and the second is `tgt_bit`.
        fn apply_2q(
            state: &mut [Complex64],
            m: &[Complex64; 16],
            ctrl_bit: usize,
            tgt_bit: usize,
            dim: usize,
        ) {
            for i in 0..dim {
                // Only process the "00" representative of each 4-element group
                if (i >> ctrl_bit) & 1 == 0 && (i >> tgt_bit) & 1 == 0 {
                    let i00 = i;
                    let i01 = i | (1 << tgt_bit);
                    let i10 = i | (1 << ctrl_bit);
                    let i11 = i | (1 << ctrl_bit) | (1 << tgt_bit);
                    let s00 = state[i00];
                    let s01 = state[i01];
                    let s10 = state[i10];
                    let s11 = state[i11];
                    state[i00] = m[0] * s00 + m[1] * s01 + m[2] * s10 + m[3] * s11;
                    state[i01] = m[4] * s00 + m[5] * s01 + m[6] * s10 + m[7] * s11;
                    state[i10] = m[8] * s00 + m[9] * s01 + m[10] * s10 + m[11] * s11;
                    state[i11] = m[12] * s00 + m[13] * s01 + m[14] * s10 + m[15] * s11;
                }
            }
        }

        match gate.gate_type.as_str() {
            // ─── single-qubit gates ───────────────────────────────────────
            "H" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("H gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let f = std::f64::consts::FRAC_1_SQRT_2;
                let c = Complex64::new(f, 0.0);
                let m = [c, c, c, -c];
                apply_1q(state, &m, q, dim);
            }
            "X" | "PauliX" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("X gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let m = [
                    Complex64::new(0.0, 0.0),
                    Complex64::new(1.0, 0.0),
                    Complex64::new(1.0, 0.0),
                    Complex64::new(0.0, 0.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            "Y" | "PauliY" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("Y gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let m = [
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, -1.0),
                    Complex64::new(0.0, 1.0),
                    Complex64::new(0.0, 0.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            "Z" | "PauliZ" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("Z gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let m = [
                    Complex64::new(1.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(-1.0, 0.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            "S" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("S gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let m = [
                    Complex64::new(1.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, 1.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            "T" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("T gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let f = std::f64::consts::FRAC_1_SQRT_2;
                let m = [
                    Complex64::new(1.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(f, f),
                ];
                apply_1q(state, &m, q, dim);
            }
            "SX" | "SqrtX" => {
                // √X = ½ [[1+i, 1-i], [1-i, 1+i]]
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("SX gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let pp = Complex64::new(0.5, 0.5);
                let pm = Complex64::new(0.5, -0.5);
                let m = [pp, pm, pm, pp];
                apply_1q(state, &m, q, dim);
            }
            "SqrtY" => {
                // √Y = ½ [[1+i, -1-i], [1+i,  1+i]]
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("SqrtY gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let pp = Complex64::new(0.5, 0.5);
                let nm = Complex64::new(-0.5, -0.5);
                let m = [pp, nm, pp, pp];
                apply_1q(state, &m, q, dim);
            }
            "SqrtW" => {
                // √W where W = (X+Y)/√2; matrix = [[cos(π/4), -i·sin(π/4)·e^(-iπ/4)], [i·sin(π/4)·e^(iπ/4), cos(π/4)]]
                // A common definition for √W used in Sycamore benchmarking:
                // √W = [[1, -√(i)/√2], [√(-i)/√2, 1]] / √2, simplified to a unitary below.
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("SqrtW gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                // Use: W = (X+Y)/√2, so W matrix = 1/√2 * [[0,1-i],[1+i,0]]
                // √W unitary (one principal square root):
                let f = 0.5_f64;
                let m = [
                    Complex64::new(f, f),
                    Complex64::new(f, -f),
                    Complex64::new(f, f),
                    Complex64::new(-f, f),
                ];
                apply_1q(state, &m, q, dim);
            }
            "RZ" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("RZ gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let angle = gate.parameters.first().copied().unwrap_or(0.0);
                // RZ(θ) = [[e^(-iθ/2), 0], [0, e^(iθ/2)]]
                let half = angle / 2.0;
                let m = [
                    Complex64::new((-half).cos(), (-half).sin()),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(0.0, 0.0),
                    Complex64::new(half.cos(), half.sin()),
                ];
                apply_1q(state, &m, q, dim);
            }
            "RX" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("RX gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let angle = gate.parameters.first().copied().unwrap_or(0.0);
                let half = angle / 2.0;
                let c = half.cos();
                let s = half.sin();
                let m = [
                    Complex64::new(c, 0.0),
                    Complex64::new(0.0, -s),
                    Complex64::new(0.0, -s),
                    Complex64::new(c, 0.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            "RY" => {
                let q = *gate.qubits.first().ok_or_else(|| {
                    SimulatorError::InvalidInput("RY gate requires a qubit argument".to_string())
                })?;
                if q >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index {q} out of range ({n} qubits)"
                    )));
                }
                let angle = gate.parameters.first().copied().unwrap_or(0.0);
                let half = angle / 2.0;
                let c = half.cos();
                let s = half.sin();
                let m = [
                    Complex64::new(c, 0.0),
                    Complex64::new(-s, 0.0),
                    Complex64::new(s, 0.0),
                    Complex64::new(c, 0.0),
                ];
                apply_1q(state, &m, q, dim);
            }
            // ─── two-qubit gates ──────────────────────────────────────────
            "CNOT" | "CX" => {
                if gate.qubits.len() < 2 {
                    return Err(SimulatorError::InvalidInput(
                        "CNOT gate requires 2 qubit arguments".to_string(),
                    ));
                }
                let ctrl = gate.qubits[0];
                let tgt = gate.qubits[1];
                if ctrl >= n || tgt >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index out of range ({n} qubits)"
                    )));
                }
                // CNOT matrix |ctrl,tgt⟩: [[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]
                let o = Complex64::new(0.0, 0.0);
                let l = Complex64::new(1.0, 0.0);
                let m: [Complex64; 16] = [l, o, o, o, o, l, o, o, o, o, o, l, o, o, l, o];
                apply_2q(state, &m, ctrl, tgt, dim);
            }
            "CZ" => {
                if gate.qubits.len() < 2 {
                    return Err(SimulatorError::InvalidInput(
                        "CZ gate requires 2 qubit arguments".to_string(),
                    ));
                }
                let ctrl = gate.qubits[0];
                let tgt = gate.qubits[1];
                if ctrl >= n || tgt >= n {
                    return Err(SimulatorError::InvalidInput(format!(
                        "Qubit index out of range ({n} qubits)"
                    )));
                }
                // CZ: flip phase of |11⟩ component
                let mask = (1 << ctrl) | (1 << tgt);
                for i in 0..dim {
                    if (i & mask) == mask {
                        state[i] = -state[i];
                    }
                }
            }
            // Unknown / custom gates: skip silently (no-op) to allow custom gate sets
            _ => {}
        }

        Ok(())
    }

    /// Simulate quantum sampling using precomputed ideal state amplitudes.
    ///
    /// For each circuit's amplitude vector we build a cumulative-probability
    /// table and draw `samples_per_circuit` bitstrings from it, faithfully
    /// mirroring the quantum probability distribution obtained by the classical
    /// state-vector simulation.
    fn simulate_quantum_sampling(
        &mut self,
        ideal_amplitudes: &[Array1<Complex64>],
    ) -> Result<Vec<Vec<Vec<u8>>>> {
        let mut all_samples = Vec::with_capacity(ideal_amplitudes.len());

        for amplitudes in ideal_amplitudes {
            let samples = self.sample_from_amplitudes(amplitudes)?;
            all_samples.push(samples);
        }

        Ok(all_samples)
    }

    /// Draw `samples_per_circuit` bitstrings from a state-amplitude vector
    /// using inverse-CDF (alias-free) sampling.
    fn sample_from_amplitudes(&mut self, amplitudes: &Array1<Complex64>) -> Result<Vec<Vec<u8>>> {
        let dim = amplitudes.len();
        let n = self.num_qubits;

        // Build probability distribution from |amplitude|²
        let mut probs: Vec<f64> = amplitudes.iter().map(|a| a.norm_sqr()).collect();

        // Normalise to guard against accumulated floating-point drift
        let total: f64 = probs.iter().sum();
        if total <= 0.0 {
            return Err(SimulatorError::InvalidInput(
                "State vector has zero norm — cannot sample".to_string(),
            ));
        }
        for p in &mut probs {
            *p /= total;
        }

        // Build cumulative distribution (CDF)
        let mut cdf = Vec::with_capacity(dim);
        let mut running = 0.0_f64;
        for p in &probs {
            running += p;
            cdf.push(running);
        }

        // Sample `samples_per_circuit` bitstrings via binary search on CDF
        let mut samples = Vec::with_capacity(self.params.samples_per_circuit);
        for _ in 0..self.params.samples_per_circuit {
            let u: f64 = self.rng.random();
            // Binary-search for the first CDF entry >= u
            let outcome = cdf.partition_point(|&c| c < u).min(dim - 1);

            // Convert outcome index to a per-qubit bitstring (MSB first)
            let mut bitstring = Vec::with_capacity(n);
            for q in (0..n).rev() {
                bitstring.push(((outcome >> q) & 1) as u8);
            }
            samples.push(bitstring);
        }

        Ok(samples)
    }

    /// Compute Linear Cross-Entropy Benchmarking (XEB) fidelity
    fn compute_linear_xeb(
        &self,
        ideal_amplitudes: &[Array1<Complex64>],
        quantum_samples: &[Vec<Vec<u8>>],
    ) -> Result<f64> {
        let mut total_fidelity = 0.0;

        for (amplitudes, samples) in ideal_amplitudes.iter().zip(quantum_samples.iter()) {
            let circuit_fidelity = self.compute_circuit_linear_xeb(amplitudes, samples)?;
            total_fidelity += circuit_fidelity;
        }

        Ok(total_fidelity / ideal_amplitudes.len() as f64)
    }

    /// Compute linear XEB for a single circuit
    fn compute_circuit_linear_xeb(
        &self,
        amplitudes: &Array1<Complex64>,
        samples: &[Vec<u8>],
    ) -> Result<f64> {
        let dim = amplitudes.len();
        let uniform_prob = 1.0 / dim as f64;

        let mut sum_probs = 0.0;

        for sample in samples {
            let bitstring_index = self.bitstring_to_index(sample);
            if bitstring_index < dim {
                let prob = amplitudes[bitstring_index].norm_sqr();
                sum_probs += prob;
            }
        }

        let mean_prob = sum_probs / samples.len() as f64;
        let linear_xeb = (mean_prob - uniform_prob) / uniform_prob;

        Ok(linear_xeb)
    }

    /// Convert bit string to state index
    fn bitstring_to_index(&self, bitstring: &[u8]) -> usize {
        let mut index = 0;
        for (i, &bit) in bitstring.iter().enumerate() {
            if bit == 1 {
                index |= 1 << (self.num_qubits - 1 - i);
            }
        }
        index
    }

    /// Analyze Porter-Thomas distribution
    fn analyze_porter_thomas(
        &self,
        ideal_amplitudes: &[Array1<Complex64>],
    ) -> Result<PorterThomasResult> {
        // Collect all probability amplitudes
        let mut all_probs = Vec::new();

        for amplitudes in ideal_amplitudes {
            for amplitude in amplitudes {
                all_probs.push(amplitude.norm_sqr());
            }
        }

        // Compute statistics
        let mean = all_probs.iter().sum::<f64>() / all_probs.len() as f64;
        let variance =
            all_probs.iter().map(|p| (p - mean).powi(2)).sum::<f64>() / all_probs.len() as f64;

        // Expected values for Porter-Thomas distribution
        let expected_mean = 1.0 / f64::from(1 << self.num_qubits);
        let expected_variance = expected_mean.powi(2);

        // Chi-squared test
        let bins = self.params.pt_bins;
        let mut observed = vec![0; bins];
        let mut expected = vec![0.0; bins];

        // Create histogram
        let max_prob = all_probs.iter().copied().fold(0.0f64, f64::max);
        for &prob in &all_probs {
            let bin = ((prob / max_prob) * (bins - 1) as f64) as usize;
            observed[bin.min(bins - 1)] += 1;
        }

        // Expected counts for Porter-Thomas
        let total_samples = all_probs.len() as f64;
        for i in 0..bins {
            let x = (i as f64 + 0.5) / bins as f64 * max_prob;
            // Porter-Thomas: p(P) = N * exp(-N*P) where N = 2^n
            let n = 1 << self.num_qubits;
            expected[i] = total_samples * f64::from(n) * (f64::from(-n) * x).exp() / bins as f64;
        }

        // Chi-squared statistic
        let mut chi_squared = 0.0;
        let mut degrees_freedom: usize = 0;

        for i in 0..bins {
            if expected[i] > 5.0 {
                // Only use bins with sufficient expected count
                chi_squared += (f64::from(observed[i]) - expected[i]).powi(2) / expected[i];
                degrees_freedom += 1;
            }
        }

        degrees_freedom = degrees_freedom.saturating_sub(1);

        // Approximate p-value (simplified)
        let p_value = self.chi_squared_p_value(chi_squared, degrees_freedom);

        // Kolmogorov-Smirnov test (simplified)
        let ks_statistic = self.kolmogorov_smirnov_test(&all_probs, expected_mean);

        Ok(PorterThomasResult {
            chi_squared,
            degrees_freedom,
            p_value,
            is_porter_thomas: p_value > self.params.significance_level,
            ks_statistic,
            mean,
            variance,
        })
    }

    /// Compute HOG (Heavy Output Generation) score
    fn compute_hog_score(
        &self,
        ideal_amplitudes: &[Array1<Complex64>],
        quantum_samples: &[Vec<Vec<u8>>],
    ) -> Result<f64> {
        let mut total_heavy_fraction = 0.0;

        for (amplitudes, samples) in ideal_amplitudes.iter().zip(quantum_samples.iter()) {
            let heavy_fraction = self.compute_circuit_hog_score(amplitudes, samples)?;
            total_heavy_fraction += heavy_fraction;
        }

        Ok(total_heavy_fraction / ideal_amplitudes.len() as f64)
    }

    /// Compute HOG score for single circuit
    fn compute_circuit_hog_score(
        &self,
        amplitudes: &Array1<Complex64>,
        samples: &[Vec<u8>],
    ) -> Result<f64> {
        // Find median probability
        let mut probs: Vec<f64> = amplitudes
            .iter()
            .map(scirs2_core::Complex::norm_sqr)
            .collect();
        probs.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let median_prob = probs[probs.len() / 2];

        // Count heavy outputs (above median)
        let mut heavy_count = 0;

        for sample in samples {
            let index = self.bitstring_to_index(sample);
            if index < amplitudes.len() {
                let prob = amplitudes[index].norm_sqr();
                if prob > median_prob {
                    heavy_count += 1;
                }
            }
        }

        Ok(f64::from(heavy_count) / samples.len() as f64)
    }

    /// Compute cross-entropy difference
    fn compute_cross_entropy_difference(
        &self,
        ideal_amplitudes: &[Array1<Complex64>],
        quantum_samples: &[Vec<Vec<u8>>],
    ) -> Result<f64> {
        let mut total_difference = 0.0;

        for (amplitudes, samples) in ideal_amplitudes.iter().zip(quantum_samples.iter()) {
            let difference = self.compute_circuit_cross_entropy(amplitudes, samples)?;
            total_difference += difference;
        }

        Ok(total_difference / ideal_amplitudes.len() as f64)
    }

    /// Compute cross-entropy for single circuit
    fn compute_circuit_cross_entropy(
        &self,
        amplitudes: &Array1<Complex64>,
        samples: &[Vec<u8>],
    ) -> Result<f64> {
        let dim = amplitudes.len();
        let uniform_entropy = (dim as f64).ln();

        let mut quantum_entropy = 0.0;
        let mut sample_counts = HashMap::new();

        // Count sample frequencies
        for sample in samples {
            *sample_counts.entry(sample.clone()).or_insert(0) += 1;
        }

        // Compute empirical entropy
        for count in sample_counts.values() {
            let prob = f64::from(*count) / samples.len() as f64;
            if prob > 0.0 {
                quantum_entropy -= prob * prob.ln();
            }
        }

        Ok(uniform_entropy - quantum_entropy)
    }

    /// Compute statistical confidence that the observed samples are
    /// consistent with the ideal quantum distribution.
    ///
    /// Performs a real chi-squared goodness-of-fit test per circuit,
    /// comparing each bitstring's observed sample frequency against its
    /// expected probability `|amplitude|^2` under the ideal state, then
    /// aggregates the chi-squared statistics (and their degrees of
    /// freedom) across all circuits and converts the combined statistic
    /// to a p-value via [`Self::chi_squared_p_value`].
    fn compute_statistical_confidence(
        &self,
        ideal_amplitudes: &[Array1<Complex64>],
        quantum_samples: &[Vec<Vec<u8>>],
    ) -> Result<f64> {
        let mut total_chi_squared = 0.0;
        let mut total_degrees_freedom = 0usize;

        for (amplitudes, samples) in ideal_amplitudes.iter().zip(quantum_samples.iter()) {
            if samples.is_empty() {
                continue;
            }
            let dim = amplitudes.len();
            let n = samples.len() as f64;

            let mut observed_counts: HashMap<usize, usize> = HashMap::new();
            for sample in samples {
                let index = self.bitstring_to_index(sample);
                *observed_counts.entry(index).or_insert(0) += 1;
            }

            let mut chi_squared = 0.0;
            let mut categories_with_support = 0usize;
            for index in 0..dim {
                let expected_prob = amplitudes[index].norm_sqr();
                let expected_count = expected_prob * n;
                let observed_count = *observed_counts.get(&index).unwrap_or(&0) as f64;
                if expected_count > 1e-12 {
                    chi_squared += (observed_count - expected_count).powi(2) / expected_count;
                    categories_with_support += 1;
                }
            }

            total_chi_squared += chi_squared;
            total_degrees_freedom += categories_with_support.saturating_sub(1);
        }

        if total_degrees_freedom == 0 {
            // No usable statistics (e.g. every circuit had zero samples):
            // honestly report zero confidence rather than a fabricated
            // constant.
            return Ok(0.0);
        }

        Ok(self.chi_squared_p_value(total_chi_squared, total_degrees_freedom))
    }

    /// Estimate quantum execution time
    fn estimate_quantum_time(&self) -> Result<f64> {
        // Estimate based on gate count and typical gate times
        let gates_per_circuit = self.params.circuit_depth * self.num_qubits;
        let total_gates = gates_per_circuit * self.params.num_circuits;
        let gate_time = 100e-9; // 100 ns per gate
        let readout_time = 1e-6; // 1 μs readout

        Ok((total_gates as f64).mul_add(gate_time, self.params.num_circuits as f64 * readout_time))
    }

    /// Estimate classical memory usage
    const fn estimate_classical_memory(&self) -> usize {
        let state_size = (1 << self.num_qubits) * std::mem::size_of::<Complex64>();
        state_size * self.params.num_circuits
    }

    /// Simplified chi-squared p-value computation
    fn chi_squared_p_value(&self, chi_squared: f64, degrees_freedom: usize) -> f64 {
        // Very simplified approximation
        if degrees_freedom == 0 {
            return 1.0;
        }

        let expected = degrees_freedom as f64;
        if chi_squared < expected {
            0.95 // High p-value
        } else if chi_squared < 2.0 * expected {
            0.5 // Medium p-value
        } else {
            0.01 // Low p-value
        }
    }

    /// Kolmogorov-Smirnov test statistic
    fn kolmogorov_smirnov_test(&self, data: &[f64], expected_mean: f64) -> f64 {
        // Simplified implementation
        let n = data.len() as f64;
        let mut sorted_data = data.to_vec();
        sorted_data.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));

        let mut max_diff: f64 = 0.0;

        for (i, &value) in sorted_data.iter().enumerate() {
            let empirical_cdf = (i + 1) as f64 / n;
            let theoretical_cdf = 1.0 - (-value / expected_mean).exp(); // Exponential CDF
            let diff = (empirical_cdf - theoretical_cdf).abs();
            max_diff = max_diff.max(diff);
        }

        max_diff
    }
}

/// Random quantum circuit representation
#[derive(Debug, Clone)]
pub struct RandomCircuit {
    pub num_qubits: usize,
    pub layers: Vec<CircuitLayer>,
    pub measurement_pattern: Vec<usize>,
}

/// Layer of gates in a circuit
#[derive(Debug, Clone)]
pub struct CircuitLayer {
    pub gates: Vec<QuantumGate>,
}

/// Quantum gate representation
#[derive(Debug, Clone)]
pub struct QuantumGate {
    pub gate_type: String,
    pub qubits: Vec<usize>,
    pub parameters: Vec<f64>,
}

/// Benchmark quantum supremacy verification
pub fn benchmark_quantum_supremacy(num_qubits: usize) -> Result<CrossEntropyResult> {
    let params = VerificationParams {
        num_circuits: 10,
        samples_per_circuit: 100,
        circuit_depth: 10,
        ..Default::default()
    };

    let mut verifier = QuantumSupremacyVerifier::new(num_qubits, params)?;
    verifier.verify_quantum_supremacy()
}

/// Verify a specific quantum supremacy claim against real experimental data.
///
/// Unlike a benchmark (which generates its own random circuits), verifying
/// an actual claim requires the exact circuit that produced
/// `experimental_data` -- there is no way to honestly cross-entropy-check
/// samples against a circuit we don't know. `circuit` must therefore be
/// supplied by the caller (e.g. reconstructed from the published circuit
/// description of the claim being verified); this function classically
/// simulates it to get the ideal amplitudes and then runs the real
/// Linear XEB / Porter-Thomas / HOG / chi-squared analysis against
/// `experimental_data`.
pub fn verify_supremacy_claim(
    num_qubits: usize,
    circuit_depth: usize,
    circuit: &RandomCircuit,
    experimental_data: &[Vec<u8>],
) -> Result<CrossEntropyResult> {
    if circuit.num_qubits != num_qubits {
        return Err(SimulatorError::InvalidConfiguration(format!(
            "supplied circuit acts on {} qubits, but verification was requested for {num_qubits} qubits",
            circuit.num_qubits
        )));
    }
    if experimental_data.is_empty() {
        return Err(SimulatorError::InvalidInput(
            "verify_supremacy_claim requires at least one experimental sample".to_string(),
        ));
    }

    let start_time = std::time::Instant::now();

    let params = VerificationParams {
        num_circuits: 1,
        samples_per_circuit: experimental_data.len(),
        circuit_depth,
        ..Default::default()
    };

    let verifier = QuantumSupremacyVerifier::new(num_qubits, params)?;

    let circuits = vec![circuit.clone()];
    let ideal_amplitudes = verifier.compute_ideal_amplitudes(&circuits)?;
    let quantum_samples = vec![experimental_data.to_vec()];

    let linear_xeb = verifier.compute_linear_xeb(&ideal_amplitudes, &quantum_samples)?;
    let porter_thomas = verifier.analyze_porter_thomas(&ideal_amplitudes)?;
    let hog_score = verifier.compute_hog_score(&ideal_amplitudes, &quantum_samples)?;
    let cross_entropy_diff =
        verifier.compute_cross_entropy_difference(&ideal_amplitudes, &quantum_samples)?;
    let confidence =
        verifier.compute_statistical_confidence(&ideal_amplitudes, &quantum_samples)?;

    let classical_time = start_time.elapsed().as_secs_f64();
    let quantum_time = verifier.estimate_quantum_time()?;

    let cost_comparison = CostComparison {
        classical_time,
        quantum_time,
        classical_memory: verifier.estimate_classical_memory(),
        speedup_factor: if quantum_time > 0.0 {
            classical_time / quantum_time
        } else {
            0.0
        },
        operation_count: circuit_depth * num_qubits,
    };

    Ok(CrossEntropyResult {
        linear_xeb_fidelity: linear_xeb,
        cross_entropy_difference: cross_entropy_diff,
        num_samples: experimental_data.len(),
        confidence,
        porter_thomas,
        hog_score,
        cost_comparison,
    })
}

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

    #[test]
    fn test_verifier_creation() {
        let params = VerificationParams::default();
        let verifier = QuantumSupremacyVerifier::new(10, params);
        assert!(verifier.is_ok());
    }

    #[test]
    fn test_linear_xeb_calculation() {
        let mut verifier = QuantumSupremacyVerifier::new(3, VerificationParams::default())
            .expect("Failed to create verifier");

        // Create dummy data
        let amplitudes = Array1::from_vec(vec![
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
            Complex64::new(0.0, 0.0),
            Complex64::new(0.0, 0.0),
            Complex64::new(0.0, 0.0),
            Complex64::new(0.0, 0.0),
        ]);

        let samples = vec![vec![0, 0, 0], vec![0, 0, 1], vec![0, 1, 0], vec![0, 1, 1]];

        let xeb = verifier
            .compute_circuit_linear_xeb(&amplitudes, &samples)
            .expect("Failed to compute linear XEB");
        assert!(xeb >= 0.0);
    }

    #[test]
    fn test_bitstring_conversion() {
        let verifier = QuantumSupremacyVerifier::new(3, VerificationParams::default())
            .expect("Failed to create verifier");

        assert_eq!(verifier.bitstring_to_index(&[0, 0, 0]), 0);
        assert_eq!(verifier.bitstring_to_index(&[0, 0, 1]), 1);
        assert_eq!(verifier.bitstring_to_index(&[1, 1, 1]), 7);
    }

    #[test]
    fn test_porter_thomas_analysis() {
        let verifier = QuantumSupremacyVerifier::new(2, VerificationParams::default())
            .expect("Failed to create verifier");

        // Create uniform random amplitudes
        let amplitudes = vec![Array1::from_vec(vec![
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
            Complex64::new(0.5, 0.0),
        ])];

        let result = verifier
            .analyze_porter_thomas(&amplitudes)
            .expect("Failed to analyze Porter-Thomas distribution");
        assert!(result.chi_squared >= 0.0);
        assert!(result.p_value >= 0.0 && result.p_value <= 1.0);
    }
}