diff options
Diffstat (limited to 'libpsi-core/src')
| -rw-r--r-- | libpsi-core/src/core/circuit.rs | 481 | ||||
| -rw-r--r-- | libpsi-core/src/core/classical_components.rs | 63 | ||||
| -rw-r--r-- | libpsi-core/src/core/custom_gate.rs | 245 | ||||
| -rw-r--r-- | libpsi-core/src/core/gates.rs | 253 | ||||
| -rw-r--r-- | libpsi-core/src/core/kernel.rs | 622 | ||||
| -rw-r--r-- | libpsi-core/src/core/mod.rs | 17 | ||||
| -rw-r--r-- | libpsi-core/src/core/noise.rs | 560 | ||||
| -rw-r--r-- | libpsi-core/src/core/quantum_components.rs | 318 | ||||
| -rw-r--r-- | libpsi-core/src/core/runtime.rs | 585 | ||||
| -rw-r--r-- | libpsi-core/src/lib.rs | 18 | ||||
| -rw-r--r-- | libpsi-core/src/maths/complex.rs | 180 | ||||
| -rw-r--r-- | libpsi-core/src/maths/format.rs | 141 | ||||
| -rw-r--r-- | libpsi-core/src/maths/matrix.rs | 310 | ||||
| -rw-r--r-- | libpsi-core/src/maths/mod.rs | 14 | ||||
| -rw-r--r-- | libpsi-core/src/maths/numeric.rs | 101 | ||||
| -rw-r--r-- | libpsi-core/src/maths/simd.rs | 510 | ||||
| -rw-r--r-- | libpsi-core/src/maths/vector.rs | 258 | ||||
| -rw-r--r-- | libpsi-core/src/maths/vector_ops.rs | 107 |
18 files changed, 0 insertions, 4783 deletions
diff --git a/libpsi-core/src/core/circuit.rs b/libpsi-core/src/core/circuit.rs deleted file mode 100644 index d0e4450..0000000 --- a/libpsi-core/src/core/circuit.rs +++ /dev/null @@ -1,481 +0,0 @@ -use super::{CustomGate, QuantumState, Runtime, RuntimeConfig}; -use crate::{format_amplitude, format_probability, Vector}; -use core::fmt; -use std::sync::Arc; - -#[derive(Clone)] -pub enum GateOp { - H(usize), - X(usize), - Y(usize), - Z(usize), - S(usize), - T(usize), - Sdg(usize), - Tdg(usize), - Sx(usize), - Sxdg(usize), - Rx(usize, f64), - Ry(usize, f64), - Rz(usize, f64), - P(usize, f64), - U1(usize, f64), - U2(usize, f64, f64), - U3(usize, f64, f64, f64), - CNOT(usize, usize), - CZ(usize, usize), - SWAP(usize, usize), - CRx(usize, usize, f64), - CRy(usize, usize, f64), - CRz(usize, usize, f64), - CP(usize, usize, f64), - CCNOT(usize, usize, usize), - CSWAP(usize, usize, usize), - Measure(usize, usize), - Custom(Arc<CustomGate>, Vec<usize>), -} - -impl GateOp { - pub fn name(&self) -> &str { - match self { - GateOp::H(_) => "H", - GateOp::X(_) => "X", - GateOp::Y(_) => "Y", - GateOp::Z(_) => "Z", - GateOp::S(_) => "S", - GateOp::T(_) => "T", - GateOp::Sdg(_) => "S†", - GateOp::Tdg(_) => "T†", - GateOp::Sx(_) => "√X", - GateOp::Sxdg(_) => "√X†", - GateOp::Rx(_, _) => "Rx", - GateOp::Ry(_, _) => "Ry", - GateOp::Rz(_, _) => "Rz", - GateOp::P(_, _) => "P", - GateOp::U1(_, _) => "U1", - GateOp::U2(_, _, _) => "U2", - GateOp::U3(_, _, _, _) => "U3", - GateOp::CRx(_, _, _) => "CRx", - GateOp::CRy(_, _, _) => "CRy", - GateOp::CRz(_, _, _) => "CRz", - GateOp::CP(_, _, _) => "CP", - GateOp::CNOT(_, _) => "CNOT", - GateOp::CZ(_, _) => "CZ", - GateOp::SWAP(_, _) => "SWAP", - GateOp::CCNOT(_, _, _) => "CCNOT", - GateOp::CSWAP(_, _, _) => "CSWAP", - GateOp::Measure(_, _) => "M", - GateOp::Custom(gate, _) => &gate.name, - } - } - - pub fn quantum_targets(&self) -> Vec<usize> { - match self { - GateOp::H(t) - | GateOp::X(t) - | GateOp::Y(t) - | GateOp::Z(t) - | GateOp::S(t) - | GateOp::T(t) - | GateOp::Sdg(t) - | GateOp::Tdg(t) - | GateOp::Sx(t) - | GateOp::Sxdg(t) - | GateOp::Rx(t, _) - | GateOp::Ry(t, _) - | GateOp::Rz(t, _) - | GateOp::P(t, _) - | GateOp::U1(t, _) - | GateOp::U2(t, _, _) - | GateOp::U3(t, _, _, _) => vec![*t], - GateOp::CNOT(c, t) - | GateOp::CZ(c, t) - | GateOp::SWAP(c, t) - | GateOp::CRx(c, t, _) - | GateOp::CRy(c, t, _) - | GateOp::CRz(c, t, _) - | GateOp::CP(c, t, _) => vec![*c, *t], - GateOp::CCNOT(c1, c2, t) | GateOp::CSWAP(c1, c2, t) => vec![*c1, *c2, *t], - GateOp::Measure(q, _) => vec![*q], - GateOp::Custom(_, targets) => targets.clone(), - } - } - - pub fn classical_targets(&self) -> Vec<usize> { - match self { - GateOp::Measure(_, c) => vec![*c], - _ => vec![], - } - } - - pub fn is_measurement(&self) -> bool { - matches!(self, GateOp::Measure(_, _)) - } - - pub fn is_custom(&self) -> bool { - matches!(self, GateOp::Custom(_, _)) - } - - pub fn is_non_clifford(&self) -> bool { - matches!( - self, - GateOp::T(_) - | GateOp::Tdg(_) - | GateOp::Sx(_) - | GateOp::Sxdg(_) - | GateOp::Rx(_, _) - | GateOp::Ry(_, _) - | GateOp::Rz(_, _) - | GateOp::P(_, _) - | GateOp::U1(_, _) - | GateOp::U2(_, _, _) - | GateOp::U3(_, _, _, _) - | GateOp::CRx(_, _, _) - | GateOp::CRy(_, _, _) - | GateOp::CRz(_, _, _) - | GateOp::CP(_, _, _) - ) - } -} - -pub struct QuantumCircuit { - num_qubits: usize, - num_classical: usize, - operations: Vec<GateOp>, - computed_state: Option<QuantumState>, -} - -impl QuantumCircuit { - pub fn new(num_qubits: usize) -> QuantumCircuit { - QuantumCircuit { - num_qubits, - num_classical: 0, - operations: Vec::new(), - computed_state: None, - } - } - - pub fn with_classical(num_qubits: usize, num_classical: usize) -> QuantumCircuit { - QuantumCircuit { - num_qubits, - num_classical, - operations: Vec::new(), - computed_state: None, - } - } - - pub fn num_qubits(&self) -> usize { - self.num_qubits - } - - pub fn num_classical(&self) -> usize { - self.num_classical - } - - pub fn operations(&self) -> &[GateOp] { - &self.operations - } - - pub fn is_computed(&self) -> bool { - self.computed_state.is_some() - } - - pub fn compute(&mut self) -> &QuantumState { - self.compute_with(Runtime::default()) - } - - pub fn compute_with(&mut self, runtime: Runtime) -> &QuantumState { - if self.computed_state.is_some() { - return self.computed_state.as_ref().unwrap(); - } - - self.computed_state = Some(runtime.compute(self.num_qubits, &self.operations)); - self.computed_state.as_ref().unwrap() - } - - pub fn compute_with_config(&mut self, config: RuntimeConfig) -> &QuantumState { - if self.computed_state.is_some() { - return self.computed_state.as_ref().unwrap(); - } - - self.computed_state = Some(config.compute(self.num_qubits, &self.operations)); - self.computed_state.as_ref().unwrap() - } - - pub fn state(&mut self) -> &QuantumState { - self.compute() - } - - pub fn state_with(&mut self, runtime: Runtime) -> &QuantumState { - self.compute_with(runtime) - } - - pub fn state_with_config(&mut self, config: RuntimeConfig) -> &QuantumState { - self.compute_with_config(config) - } - - pub fn h(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::H(target)); - self.computed_state = None; - self - } - - pub fn x(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::X(target)); - self.computed_state = None; - self - } - - pub fn y(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Y(target)); - self.computed_state = None; - self - } - - pub fn z(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Z(target)); - self.computed_state = None; - self - } - - pub fn s(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::S(target)); - self.computed_state = None; - self - } - - pub fn t(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::T(target)); - self.computed_state = None; - self - } - - pub fn sdg(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Sdg(target)); - self.computed_state = None; - self - } - - pub fn tdg(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Tdg(target)); - self.computed_state = None; - self - } - - pub fn sx(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Sx(target)); - self.computed_state = None; - self - } - - pub fn sxdg(&mut self, target: usize) -> &mut Self { - self.operations.push(GateOp::Sxdg(target)); - self.computed_state = None; - self - } - - pub fn rx(&mut self, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::Rx(target, theta)); - self.computed_state = None; - self - } - - pub fn ry(&mut self, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::Ry(target, theta)); - self.computed_state = None; - self - } - - pub fn rz(&mut self, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::Rz(target, theta)); - self.computed_state = None; - self - } - - pub fn p(&mut self, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::P(target, theta)); - self.computed_state = None; - self - } - - pub fn u1(&mut self, target: usize, lambda: f64) -> &mut Self { - self.operations.push(GateOp::U1(target, lambda)); - self.computed_state = None; - self - } - - pub fn u2(&mut self, target: usize, phi: f64, lambda: f64) -> &mut Self { - self.operations.push(GateOp::U2(target, phi, lambda)); - self.computed_state = None; - self - } - - pub fn u3(&mut self, target: usize, theta: f64, phi: f64, lambda: f64) -> &mut Self { - self.operations.push(GateOp::U3(target, theta, phi, lambda)); - self.computed_state = None; - self - } - - pub fn crx(&mut self, control: usize, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::CRx(control, target, theta)); - self.computed_state = None; - self - } - - pub fn cry(&mut self, control: usize, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::CRy(control, target, theta)); - self.computed_state = None; - self - } - - pub fn crz(&mut self, control: usize, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::CRz(control, target, theta)); - self.computed_state = None; - self - } - - pub fn cp(&mut self, control: usize, target: usize, theta: f64) -> &mut Self { - self.operations.push(GateOp::CP(control, target, theta)); - self.computed_state = None; - self - } - - pub fn cnot(&mut self, control: usize, target: usize) -> &mut Self { - self.operations.push(GateOp::CNOT(control, target)); - self.computed_state = None; - self - } - - pub fn cx(&mut self, control: usize, target: usize) -> &mut Self { - self.cnot(control, target) - } - - pub fn cz(&mut self, control: usize, target: usize) -> &mut Self { - self.operations.push(GateOp::CZ(control, target)); - self.computed_state = None; - self - } - - pub fn swap(&mut self, qubit1: usize, qubit2: usize) -> &mut Self { - self.operations.push(GateOp::SWAP(qubit1, qubit2)); - self.computed_state = None; - self - } - - pub fn ccnot(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self { - self.operations - .push(GateOp::CCNOT(control1, control2, target)); - self.computed_state = None; - self - } - - pub fn toffoli(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self { - self.ccnot(control1, control2, target) - } - - pub fn cswap(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self { - self.operations - .push(GateOp::CSWAP(control, target1, target2)); - self.computed_state = None; - self - } - - pub fn fredkin(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self { - self.cswap(control, target1, target2) - } - - pub fn measure(&mut self, qubit: usize, classical: usize) -> &mut Self { - if classical >= self.num_classical { - self.num_classical = classical + 1; - } - self.operations.push(GateOp::Measure(qubit, classical)); - self - } - - pub fn measure_all(&mut self) -> &mut Self { - for i in 0..self.num_qubits { - self.measure(i, i); - } - self - } - - pub fn custom(&mut self, gate: &Arc<CustomGate>, targets: &[usize]) -> &mut Self { - self.operations - .push(GateOp::Custom(Arc::clone(gate), targets.to_vec())); - self.computed_state = None; - self - } - - pub fn apply_custom(&mut self, gate: CustomGate, targets: &[usize]) -> &mut Self { - self.operations - .push(GateOp::Custom(Arc::new(gate), targets.to_vec())); - self.computed_state = None; - self - } - - pub fn reset(&mut self) -> &mut Self { - self.operations.clear(); - self.computed_state = None; - self - } - - pub fn probability(&mut self, state_index: usize) -> f64 { - self.compute(); - let state = self.computed_state.as_ref().unwrap(); - let amp = state.get(state_index); - amp.norm2() - } - - pub fn probabilities(&mut self) -> Vec<f64> { - self.compute(); - let n = 1 << self.num_qubits; - let state = self.computed_state.as_ref().unwrap(); - (0..n).map(|i| state.get(i).norm2()).collect() - } - - pub fn print_probabilities(&mut self) { - let probs = self.probabilities(); - let n = self.num_qubits; - println!("Probabilities:"); - for (i, p) in probs.iter().enumerate() { - if *p > 1e-10 { - let basis: String = format!("{:0width$b}", i, width = n); - println!(" |{}⟩: {}", basis, format_probability(*p)); - } - } - } -} - -impl fmt::Display for QuantumCircuit { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - writeln!( - f, - "QuantumCircuit ({} qubits, {} classical)", - self.num_qubits, self.num_classical - )?; - writeln!(f, "Operations:")?; - for (i, op) in self.operations.iter().enumerate() { - match op { - GateOp::Measure(q, c) => writeln!(f, " {}: {} q{} → c{}", i, op.name(), q, c)?, - GateOp::Custom(gate, targets) => { - writeln!(f, " {}: [{}] on {:?}", i, gate.name, targets)? - } - _ => writeln!(f, " {}: {} on {:?}", i, op.name(), op.quantum_targets())?, - } - } - if let Some(state) = &self.computed_state { - writeln!(f, "State:")?; - let n = 1 << self.num_qubits; - for i in 0..n { - let amp = state.get(i); - if amp.real.abs() > 1e-10 || amp.imaginary.abs() > 1e-10 { - let basis: String = format!("{:0width$b}", i, width = self.num_qubits); - writeln!(f, " |{}⟩: {}", basis, format_amplitude(&))?; - } - } - } else { - writeln!(f, "State: (not computed)")?; - } - Ok(()) - } -} diff --git a/libpsi-core/src/core/classical_components.rs b/libpsi-core/src/core/classical_components.rs deleted file mode 100644 index f6565af..0000000 --- a/libpsi-core/src/core/classical_components.rs +++ /dev/null @@ -1,63 +0,0 @@ -use core::ops; - -#[derive(Clone, Copy)] -pub struct ClassicalBit<'a> { - state: bool, - name: &'a str, -} - -#[derive(Clone)] -pub struct ClassicalRegister<'a> { - bits: Vec<ClassicalBit<'a>>, - name: &'a str, -} - -impl<'a> ClassicalBit<'a> { - pub fn new(name: &'a str, state: bool) -> ClassicalBit<'a> { - ClassicalBit { name, state } - } - - pub fn get_name(&self) -> &'a str { - self.name - } - - pub fn get_state(&self) -> bool { - self.state - } -} - -impl<'a> ClassicalRegister<'a> { - pub fn new(name: &'a str, names: &'a [&'a str]) -> ClassicalRegister<'a> { - let mut bits: Vec<ClassicalBit<'a>> = Vec::new(); - for i in 0..names.len() { - bits.push(ClassicalBit::new(names[i], false)); - } - ClassicalRegister { name, bits } - } - - pub fn set_bits(&mut self, bits: Vec<ClassicalBit<'a>>) { - self.bits = bits; - } - - pub fn get_bits(&self) -> Vec<ClassicalBit<'a>> { - self.bits.clone() - } - - pub fn get_name(&self) -> &'a str { - self.name - } -} - -impl<'a> ops::Index<usize> for ClassicalRegister<'a> { - type Output = ClassicalBit<'a>; - - fn index(&self, index: usize) -> &Self::Output { - &self.bits[index] - } -} - -impl<'a> ops::IndexMut<usize> for ClassicalRegister<'a> { - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.bits[index] - } -} diff --git a/libpsi-core/src/core/custom_gate.rs b/libpsi-core/src/core/custom_gate.rs deleted file mode 100644 index 9e49ac9..0000000 --- a/libpsi-core/src/core/custom_gate.rs +++ /dev/null @@ -1,245 +0,0 @@ -use crate::{Complex, Matrix, QuantumGate}; - -#[derive(Clone)] -pub enum CustomGateDefinition { - Matrix(Matrix<Complex<f64>>), - Composite(Vec<(CompositeOp, Vec<usize>)>), -} - -#[derive(Clone, Copy)] -pub enum CompositeOp { - H, - X, - Y, - Z, - S, - T, - CNOT, - CZ, - SWAP, - CCNOT, - CSWAP, -} - -#[derive(Clone)] -pub struct CustomGate { - pub name: String, - pub num_qubits: usize, - pub definition: CustomGateDefinition, -} - -impl CustomGate { - pub fn from_matrix(name: &str, matrix: Matrix<Complex<f64>>) -> Self { - let dim = matrix.rows; - let num_qubits = (dim as f64).log2() as usize; - assert_eq!( - 1 << num_qubits, - dim, - "Matrix dimension must be a power of 2" - ); - assert_eq!(matrix.rows, matrix.cols, "Matrix must be square"); - - CustomGate { - name: String::from(name), - num_qubits, - definition: CustomGateDefinition::Matrix(matrix), - } - } - - pub fn from_composite( - name: &str, - num_qubits: usize, - ops: Vec<(CompositeOp, Vec<usize>)>, - ) -> Self { - CustomGate { - name: String::from(name), - num_qubits, - definition: CustomGateDefinition::Composite(ops), - } - } - - pub fn to_quantum_gate(&self) -> QuantumGate<'static> { - match &self.definition { - CustomGateDefinition::Matrix(matrix) => { - let name: &'static str = Box::leak(self.name.clone().into_boxed_str()); - QuantumGate { - name, - matrix: matrix.clone(), - num_qubits: self.num_qubits, - } - } - CustomGateDefinition::Composite(ops) => { - let matrix = self.compute_composite_matrix(ops); - let name: &'static str = Box::leak(self.name.clone().into_boxed_str()); - QuantumGate { - name, - matrix, - num_qubits: self.num_qubits, - } - } - } - } - - fn compute_composite_matrix(&self, ops: &[(CompositeOp, Vec<usize>)]) -> Matrix<Complex<f64>> { - use crate::gates::*; - use crate::Complex; - - let dim = 1 << self.num_qubits; - let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]); - for i in 0..dim { - result.data[i * dim + i] = Complex::new(1.0, 0.0); - } - - for (op, targets) in ops { - let gate: &QuantumGate = match op { - CompositeOp::H => &HADAMARD, - CompositeOp::X => &PAULI_X, - CompositeOp::Y => &PAULI_Y, - CompositeOp::Z => &PAULI_Z, - CompositeOp::S => &S_GATE, - CompositeOp::T => &T_GATE, - CompositeOp::CNOT => &CNOT, - CompositeOp::CZ => &CZ, - CompositeOp::SWAP => &SWAP, - CompositeOp::CCNOT => &TOFFOLI, - CompositeOp::CSWAP => &FREDKIN, - }; - - let full_gate = build_full_operator(&gate.matrix, targets, self.num_qubits); - result = matrix_multiply(&full_gate, &result); - } - - result - } -} - -fn build_full_operator( - gate_matrix: &Matrix<Complex<f64>>, - targets: &[usize], - total_qubits: usize, -) -> Matrix<Complex<f64>> { - let dim = 1 << total_qubits; - let gate_dim = gate_matrix.rows; - let num_gate_qubits = targets.len(); - - let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]); - - for i in 0..dim { - for j in 0..dim { - let mut gate_i = 0usize; - let mut gate_j = 0usize; - let mut match_non_targets = true; - - for q in 0..total_qubits { - let bit_i = (i >> (total_qubits - 1 - q)) & 1; - let bit_j = (j >> (total_qubits - 1 - q)) & 1; - - if let Some(pos) = targets.iter().position(|&t| t == q) { - gate_i |= bit_i << (num_gate_qubits - 1 - pos); - gate_j |= bit_j << (num_gate_qubits - 1 - pos); - } else if bit_i != bit_j { - match_non_targets = false; - break; - } - } - - if match_non_targets { - result.data[i * dim + j] = gate_matrix.data[gate_i * gate_dim + gate_j]; - } - } - } - - result -} - -fn matrix_multiply(a: &Matrix<Complex<f64>>, b: &Matrix<Complex<f64>>) -> Matrix<Complex<f64>> { - let n = a.rows; - let mut result = Matrix::new(n, n, vec![Complex::new(0.0, 0.0); n * n]); - - for i in 0..n { - for j in 0..n { - let mut sum = Complex::new(0.0, 0.0); - for k in 0..n { - sum = sum + a.data[i * n + k] * b.data[k * n + j]; - } - result.data[i * n + j] = sum; - } - } - - result -} - -pub struct CustomGateBuilder { - name: String, - num_qubits: usize, - ops: Vec<(CompositeOp, Vec<usize>)>, -} - -impl CustomGateBuilder { - pub fn new(name: &str, num_qubits: usize) -> Self { - CustomGateBuilder { - name: String::from(name), - num_qubits, - ops: Vec::new(), - } - } - - pub fn h(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::H, vec![target])); - self - } - - pub fn x(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::X, vec![target])); - self - } - - pub fn y(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::Y, vec![target])); - self - } - - pub fn z(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::Z, vec![target])); - self - } - - pub fn s(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::S, vec![target])); - self - } - - pub fn t(mut self, target: usize) -> Self { - self.ops.push((CompositeOp::T, vec![target])); - self - } - - pub fn cnot(mut self, control: usize, target: usize) -> Self { - self.ops.push((CompositeOp::CNOT, vec![control, target])); - self - } - - pub fn cz(mut self, control: usize, target: usize) -> Self { - self.ops.push((CompositeOp::CZ, vec![control, target])); - self - } - - pub fn swap(mut self, a: usize, b: usize) -> Self { - self.ops.push((CompositeOp::SWAP, vec![a, b])); - self - } - - pub fn ccnot(mut self, c1: usize, c2: usize, target: usize) -> Self { - self.ops.push((CompositeOp::CCNOT, vec![c1, c2, target])); - self - } - - pub fn cswap(mut self, control: usize, t1: usize, t2: usize) -> Self { - self.ops.push((CompositeOp::CSWAP, vec![control, t1, t2])); - self - } - - pub fn build(self) -> CustomGate { - CustomGate::from_composite(&self.name, self.num_qubits, self.ops) - } -} diff --git a/libpsi-core/src/core/gates.rs b/libpsi-core/src/core/gates.rs deleted file mode 100644 index cee4ebc..0000000 --- a/libpsi-core/src/core/gates.rs +++ /dev/null @@ -1,253 +0,0 @@ -use crate::{complex, matrix, Complex, Matrix, QuantumGate}; -use std::f64::consts::FRAC_1_SQRT_2; - -pub fn rx_matrix(theta: f64) -> Matrix<Complex<f64>> { - let cos = (theta / 2.0).cos(); - let sin = (theta / 2.0).sin(); - matrix!( - [complex!(cos, 0.0), complex!(0.0, -sin)]; - [complex!(0.0, -sin), complex!(cos, 0.0)] - ) -} - -pub fn ry_matrix(theta: f64) -> Matrix<Complex<f64>> { - let cos = (theta / 2.0).cos(); - let sin = (theta / 2.0).sin(); - matrix!( - [complex!(cos, 0.0), complex!(-sin, 0.0)]; - [complex!(sin, 0.0), complex!(cos, 0.0)] - ) -} - -pub fn rz_matrix(theta: f64) -> Matrix<Complex<f64>> { - let half = theta / 2.0; - matrix!( - [complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(half.cos(), half.sin())] - ) -} - -pub fn p_matrix(theta: f64) -> Matrix<Complex<f64>> { - matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())] - ) -} - -pub fn u1_matrix(lambda: f64) -> Matrix<Complex<f64>> { - p_matrix(lambda) -} - -pub fn u2_matrix(phi: f64, lambda: f64) -> Matrix<Complex<f64>> { - let inv_sqrt2 = FRAC_1_SQRT_2; - matrix!( - [complex!(inv_sqrt2, 0.0), complex!(-inv_sqrt2 * lambda.cos(), -inv_sqrt2 * lambda.sin())]; - [complex!(inv_sqrt2 * phi.cos(), inv_sqrt2 * phi.sin()), complex!((phi + lambda).cos() * inv_sqrt2, (phi + lambda).sin() * inv_sqrt2)] - ) -} - -pub fn u3_matrix(theta: f64, phi: f64, lambda: f64) -> Matrix<Complex<f64>> { - let cos = (theta / 2.0).cos(); - let sin = (theta / 2.0).sin(); - matrix!( - [complex!(cos, 0.0), complex!(-sin * lambda.cos(), -sin * lambda.sin())]; - [complex!(sin * phi.cos(), sin * phi.sin()), complex!(cos * (phi + lambda).cos(), cos * (phi + lambda).sin())] - ) -} - -pub fn crx_matrix(theta: f64) -> Matrix<Complex<f64>> { - let cos = (theta / 2.0).cos(); - let sin = (theta / 2.0).sin(); - matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(0.0, -sin)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, -sin), complex!(cos, 0.0)] - ) -} - -pub fn cry_matrix(theta: f64) -> Matrix<Complex<f64>> { - let cos = (theta / 2.0).cos(); - let sin = (theta / 2.0).sin(); - matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(-sin, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(sin, 0.0), complex!(cos, 0.0)] - ) -} - -pub fn crz_matrix(theta: f64) -> Matrix<Complex<f64>> { - let half = theta / 2.0; - matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), half.sin())] - ) -} - -pub fn cp_matrix(theta: f64) -> Matrix<Complex<f64>> { - matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())] - ) -} - -#[rustfmt::skip] -lazy_static::lazy_static! { - pub static ref HADAMARD: QuantumGate<'static> = QuantumGate { - name: "H", - matrix: matrix!([complex!(1.0, 0.0), complex!( 1.0, 0.0)]; - [complex!(1.0, 0.0), complex!(-1.0, 0.0)]) * - complex!(1.0/2.0_f64.sqrt(), 0.0), - num_qubits: 1, - }; - - pub static ref PAULI_X: QuantumGate<'static> = QuantumGate { - name: "X", - matrix: matrix!([complex!(0.0, 0.0), complex!(1.0, 0.0)]; - [complex!(1.0, 0.0), complex!(0.0, 0.0)]), - num_qubits: 1, - }; - - pub static ref PAULI_Y: QuantumGate<'static> = QuantumGate { - name: "Y", - matrix: matrix!([complex!(0.0, 0.0), complex!(0.0, -1.0)]; - [complex!(0.0, 1.0), complex!(0.0, 0.0)]), - num_qubits: 1, - }; - - pub static ref PAULI_Z: QuantumGate<'static> = QuantumGate { - name: "Z", - matrix: matrix!([complex!(1.0, 0.0), complex!( 0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(-1.0, 0.0)]), - num_qubits: 1, - }; - - pub static ref S_GATE: QuantumGate<'static> = QuantumGate { - name: "S", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 1.0)]), - num_qubits: 1, - }; - - pub static ref T_GATE: QuantumGate<'static> = QuantumGate { - name: "T", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, core::f64::consts::FRAC_1_SQRT_2)]), - num_qubits: 1, - }; - - pub static ref SDG_GATE: QuantumGate<'static> = QuantumGate { - name: "S†", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, -1.0)]), - num_qubits: 1, - }; - - pub static ref TDG_GATE: QuantumGate<'static> = QuantumGate { - name: "T†", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, -core::f64::consts::FRAC_1_SQRT_2)]), - num_qubits: 1, - }; - - pub static ref SX_GATE: QuantumGate<'static> = QuantumGate { - name: "√X", - matrix: matrix!([complex!(0.5, 0.5), complex!(0.5, -0.5)]; - [complex!(0.5, -0.5), complex!(0.5, 0.5)]), - num_qubits: 1, - }; - - pub static ref SXDG_GATE: QuantumGate<'static> = QuantumGate { - name: "√X†", - matrix: matrix!([complex!(0.5, -0.5), complex!(0.5, 0.5)]; - [complex!(0.5, 0.5), complex!(0.5, -0.5)]), - num_qubits: 1, - }; - - pub static ref IDENTITY: QuantumGate<'static> = QuantumGate { - name: "I", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0)]), - num_qubits: 1, - }; - - pub static ref CNOT: QuantumGate<'static> = QuantumGate { - name: "CNOT", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]), - num_qubits: 2, - }; - - pub static ref CZ: QuantumGate<'static> = QuantumGate { - name: "CZ", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!( 0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(-1.0, 0.0)]), - num_qubits: 2, - }; - - pub static ref SWAP: QuantumGate<'static> = QuantumGate { - name: "SWAP", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]), - num_qubits: 2, - }; - - pub static ref ISWAP: QuantumGate<'static> = QuantumGate { - name: "iSWAP", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]), - num_qubits: 2, - }; - - pub static ref SQRT_SWAP: QuantumGate<'static> = QuantumGate { - name: "√SWAP", - matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.5, 0.5), complex!(0.5, -0.5), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.5, -0.5), complex!(0.5, 0.5), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]), - num_qubits: 2, - }; - - pub static ref TOFFOLI: QuantumGate<'static> = QuantumGate { - name: "CCNOT", - matrix: matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)] - ), - num_qubits: 3, - }; - - pub static ref FREDKIN: QuantumGate<'static> = QuantumGate { - name: "CSWAP", - matrix: matrix!( - [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)]; - [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)] - ), - num_qubits: 3, - }; -} diff --git a/libpsi-core/src/core/kernel.rs b/libpsi-core/src/core/kernel.rs deleted file mode 100644 index 7b66eea..0000000 --- a/libpsi-core/src/core/kernel.rs +++ /dev/null @@ -1,622 +0,0 @@ -use crate::maths::simd::{ - apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel, SimdCapability, -}; -use crate::{complex, Complex, Matrix}; -use rayon::prelude::*; -use std::collections::HashSet; - -#[derive(Clone, Copy, Debug, PartialEq, Eq)] -pub enum GateType { - Diagonal, - NonDiagonal, - Controlled, -} - -#[derive(Clone)] -pub struct Kernel { - pub matrix: Matrix<Complex<f64>>, - pub targets: Vec<usize>, - pub name: String, - pub gate_type: GateType, -} - -impl Kernel { - pub fn new(name: &str, matrix: Matrix<Complex<f64>>, targets: Vec<usize>) -> Self { - let gate_type = Self::detect_gate_type(name, &matrix); - Self { - matrix, - targets, - name: name.to_string(), - gate_type, - } - } - - fn detect_gate_type(name: &str, matrix: &Matrix<Complex<f64>>) -> GateType { - let diagonal_gates = [ - "Z", "S", "T", "Sdg", "Tdg", "Rz", "P", "U1", "CZ", "CP", "CRz", - ]; - if diagonal_gates.iter().any(|&g| name.starts_with(g)) { - return GateType::Diagonal; - } - - let controlled_gates = [ - "CNOT", "CZ", "SWAP", "CRx", "CRy", "CRz", "CP", "CCNOT", "CSWAP", - ]; - if controlled_gates.iter().any(|&g| name.starts_with(g)) { - return GateType::Controlled; - } - - if matrix.rows == 2 && matrix.cols == 2 { - let is_diag = matrix.data[1].real.abs() < 1e-10 - && matrix.data[1].imaginary.abs() < 1e-10 - && matrix.data[2].real.abs() < 1e-10 - && matrix.data[2].imaginary.abs() < 1e-10; - if is_diag { - return GateType::Diagonal; - } - } - - GateType::NonDiagonal - } - - pub fn num_qubits(&self) -> usize { - self.targets.len() - } - - pub fn target_set(&self) -> HashSet<usize> { - self.targets.iter().cloned().collect() - } - - pub fn shares_qubits(&self, other: &Kernel) -> bool { - self.targets.iter().any(|t| other.targets.contains(t)) - } - - pub fn commutes_with(&self, other: &Kernel) -> bool { - if !self.shares_qubits(other) { - return true; - } - - if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal { - if self.targets == other.targets { - return true; - } - } - - false - } - - pub fn can_fuse_with(&self, other: &Kernel) -> bool { - if self.targets.len() != 1 || other.targets.len() != 1 { - return false; - } - self.targets[0] == other.targets[0] - } - - pub fn fuse(&self, other: &Kernel) -> Option<Kernel> { - if !self.can_fuse_with(other) { - return None; - } - let fused_matrix = other.matrix.dot(&self.matrix)?; - let new_type = - if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal { - GateType::Diagonal - } else { - GateType::NonDiagonal - }; - Some(Kernel { - matrix: fused_matrix, - targets: self.targets.clone(), - name: format!("{}+{}", self.name, other.name), - gate_type: new_type, - }) - } -} - -pub struct KernelBatch { - kernels: Vec<Kernel>, - num_qubits: usize, -} - -impl KernelBatch { - pub fn new(num_qubits: usize) -> Self { - Self { - kernels: Vec::new(), - num_qubits, - } - } - - pub fn add(&mut self, kernel: Kernel) { - self.kernels.push(kernel); - } - - pub fn len(&self) -> usize { - self.kernels.len() - } - - pub fn is_empty(&self) -> bool { - self.kernels.is_empty() - } - - pub fn kernels(&self) -> &[Kernel] { - &self.kernels - } - - pub fn optimize(&mut self) { - if self.kernels.len() < 2 { - return; - } - - let mut optimized: Vec<Kernel> = Vec::with_capacity(self.kernels.len()); - let mut i = 0; - - while i < self.kernels.len() { - let current = &self.kernels[i]; - - if i + 1 < self.kernels.len() { - let next = &self.kernels[i + 1]; - if let Some(fused) = current.fuse(next) { - optimized.push(fused); - i += 2; - continue; - } - } - - optimized.push(current.clone()); - i += 1; - } - - self.kernels = optimized; - } - - pub fn execute(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - *state = apply_kernel(state, kernel, self.num_qubits); - } - } - - pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - *state = apply_kernel_parallel(state, kernel, self.num_qubits); - } - } - - pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - if kernel.targets.len() == 1 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); - } else { - *state = apply_kernel(state, kernel, self.num_qubits); - } - } - } - - pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - if kernel.targets.len() == 1 && self.num_qubits >= 10 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd_parallel( - state, - &gate, - kernel.targets[0], - self.num_qubits, - ); - } else if kernel.targets.len() == 1 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); - } else { - *state = apply_kernel_parallel(state, kernel, self.num_qubits); - } - } - } - - pub fn simd_capability(&self) -> SimdCapability { - SimdCapability::detect() - } -} - -fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] { - [ - [matrix.data[0], matrix.data[1]], - [matrix.data[2], matrix.data[3]], - ] -} - -fn apply_kernel(state: &[Complex<f64>], kernel: &Kernel, num_qubits: usize) -> Vec<Complex<f64>> { - let dim = 1 << num_qubits; - let g = kernel.targets.len(); - let gate_dim = 1 << g; - - let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect(); - - let mut non_target_mask: usize = (1 << num_qubits) - 1; - for &pos in &target_bits { - non_target_mask &= !(1 << pos); - } - - let mut new_state = vec![complex!(0.0, 0.0); dim]; - - for i in 0..dim { - let mut target_idx = 0usize; - for (k, &pos) in target_bits.iter().enumerate() { - if (i >> pos) & 1 == 1 { - target_idx |= 1 << (g - 1 - k); - } - } - - let mut sum = complex!(0.0, 0.0); - - for j in 0..gate_dim { - let gate_elem = kernel.matrix.data[target_idx * gate_dim + j]; - - if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 { - continue; - } - - let mut source_idx = i & non_target_mask; - for (k, &pos) in target_bits.iter().enumerate() { - if (j >> (g - 1 - k)) & 1 == 1 { - source_idx |= 1 << pos; - } - } - - sum = sum + gate_elem * state[source_idx]; - } - - new_state[i] = sum; - } - - new_state -} - -fn apply_kernel_parallel( - state: &[Complex<f64>], - kernel: &Kernel, - num_qubits: usize, -) -> Vec<Complex<f64>> { - let dim = 1 << num_qubits; - let g = kernel.targets.len(); - let gate_dim = 1 << g; - - let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect(); - - let mut non_target_mask: usize = (1 << num_qubits) - 1; - for &pos in &target_bits { - non_target_mask &= !(1 << pos); - } - - (0..dim) - .into_par_iter() - .map(|i| { - let mut target_idx = 0usize; - for (k, &pos) in target_bits.iter().enumerate() { - if (i >> pos) & 1 == 1 { - target_idx |= 1 << (g - 1 - k); - } - } - - let mut sum = complex!(0.0, 0.0); - - for j in 0..gate_dim { - let gate_elem = kernel.matrix.data[target_idx * gate_dim + j]; - - if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 { - continue; - } - - let mut source_idx = i & non_target_mask; - for (k, &pos) in target_bits.iter().enumerate() { - if (j >> (g - 1 - k)) & 1 == 1 { - source_idx |= 1 << pos; - } - } - - sum = sum + gate_elem * state[source_idx]; - } - - sum - }) - .collect() -} - -pub struct KernelBuilder { - num_qubits: usize, -} - -impl KernelBuilder { - pub fn new(num_qubits: usize) -> Self { - Self { num_qubits } - } - - pub fn num_qubits(&self) -> usize { - self.num_qubits - } -} - -#[derive(Clone)] -pub struct ExecutionLayer { - pub kernels: Vec<Kernel>, -} - -impl ExecutionLayer { - pub fn new() -> Self { - Self { - kernels: Vec::new(), - } - } - - pub fn can_add(&self, kernel: &Kernel) -> bool { - !self.kernels.iter().any(|k| k.shares_qubits(kernel)) - } - - pub fn add(&mut self, kernel: Kernel) { - self.kernels.push(kernel); - } - - pub fn affected_qubits(&self) -> HashSet<usize> { - self.kernels - .iter() - .flat_map(|k| k.targets.iter().cloned()) - .collect() - } -} - -impl Default for ExecutionLayer { - fn default() -> Self { - Self::new() - } -} - -pub struct StructureAwareKernelBatch { - kernels: Vec<Kernel>, - layers: Vec<ExecutionLayer>, - num_qubits: usize, - optimised: bool, -} - -impl StructureAwareKernelBatch { - pub fn new(num_qubits: usize) -> Self { - Self { - kernels: Vec::new(), - layers: Vec::new(), - num_qubits, - optimised: false, - } - } - - pub fn add(&mut self, kernel: Kernel) { - self.kernels.push(kernel); - self.optimised = false; - } - - pub fn len(&self) -> usize { - self.kernels.len() - } - - pub fn is_empty(&self) -> bool { - self.kernels.is_empty() - } - - pub fn kernels(&self) -> &[Kernel] { - &self.kernels - } - - pub fn layers(&self) -> &[ExecutionLayer] { - &self.layers - } - - pub fn num_layers(&self) -> usize { - self.layers.len() - } - - pub fn optimise(&mut self) { - if self.optimised || self.kernels.len() < 2 { - return; - } - - self.reorder_commuting_gates(); - self.multi_pass_fusion(); - self.build_execution_layers(); - self.optimised = true; - } - - fn reorder_commuting_gates(&mut self) { - let mut changed = true; - let mut iterations = 0; - const MAX_ITERATIONS: usize = 100; - - while changed && iterations < MAX_ITERATIONS { - changed = false; - iterations += 1; - - for i in 0..self.kernels.len().saturating_sub(1) { - let current = &self.kernels[i]; - let next = &self.kernels[i + 1]; - - if current.targets.len() == 1 - && next.targets.len() == 1 - && current.targets[0] != next.targets[0] - && current.commutes_with(next) - { - for j in (i + 2)..self.kernels.len() { - let candidate = &self.kernels[j]; - - if candidate.targets.len() == 1 - && candidate.targets[0] == current.targets[0] - { - let can_move = (i + 1..j).all(|k| { - let between = &self.kernels[k]; - !between.shares_qubits(current) || current.commutes_with(between) - }); - - if can_move && current.can_fuse_with(candidate) { - let kernel_to_move = self.kernels.remove(j); - self.kernels.insert(i + 1, kernel_to_move); - changed = true; - break; - } - } - } - } - } - } - } - - fn multi_pass_fusion(&mut self) { - let mut changed = true; - let mut iterations = 0; - const MAX_ITERATIONS: usize = 50; - - while changed && iterations < MAX_ITERATIONS { - changed = false; - iterations += 1; - - let mut new_kernels: Vec<Kernel> = Vec::with_capacity(self.kernels.len()); - let mut i = 0; - - while i < self.kernels.len() { - if i + 1 < self.kernels.len() { - let current = &self.kernels[i]; - let next = &self.kernels[i + 1]; - - if let Some(fused) = current.fuse(next) { - new_kernels.push(fused); - i += 2; - changed = true; - continue; - } - } - - new_kernels.push(self.kernels[i].clone()); - i += 1; - } - - self.kernels = new_kernels; - } - } - - fn build_execution_layers(&mut self) { - self.layers.clear(); - - for kernel in &self.kernels { - let mut placed = false; - - for layer in &mut self.layers { - if layer.can_add(kernel) { - layer.add(kernel.clone()); - placed = true; - break; - } - } - - if !placed { - let mut new_layer = ExecutionLayer::new(); - new_layer.add(kernel.clone()); - self.layers.push(new_layer); - } - } - } - - pub fn execute(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - *state = apply_kernel(state, kernel, self.num_qubits); - } - } - - pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - *state = apply_kernel_parallel(state, kernel, self.num_qubits); - } - } - - pub fn execute_layered(&self, state: &mut Vec<Complex<f64>>) { - for layer in &self.layers { - for kernel in &layer.kernels { - *state = apply_kernel(state, kernel, self.num_qubits); - } - } - } - - pub fn execute_layered_parallel(&self, state: &mut Vec<Complex<f64>>) { - for layer in &self.layers { - for kernel in &layer.kernels { - *state = apply_kernel_parallel(state, kernel, self.num_qubits); - } - } - } - - pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - if kernel.targets.len() == 1 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); - } else { - *state = apply_kernel(state, kernel, self.num_qubits); - } - } - } - - pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) { - for kernel in &self.kernels { - if kernel.targets.len() == 1 && self.num_qubits >= 10 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd_parallel( - state, - &gate, - kernel.targets[0], - self.num_qubits, - ); - } else if kernel.targets.len() == 1 { - let gate = matrix_to_2x2(&kernel.matrix); - apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); - } else { - *state = apply_kernel_parallel(state, kernel, self.num_qubits); - } - } - } - - pub fn stats(&self) -> KernelStats { - let single_qubit = self.kernels.iter().filter(|k| k.targets.len() == 1).count(); - let two_qubit = self.kernels.iter().filter(|k| k.targets.len() == 2).count(); - let multi_qubit = self.kernels.iter().filter(|k| k.targets.len() > 2).count(); - let diagonal = self - .kernels - .iter() - .filter(|k| k.gate_type == GateType::Diagonal) - .count(); - - KernelStats { - total_kernels: self.kernels.len(), - single_qubit, - two_qubit, - multi_qubit, - diagonal, - execution_layers: self.layers.len(), - } - } -} - -#[derive(Debug, Clone)] -pub struct KernelStats { - pub total_kernels: usize, - pub single_qubit: usize, - pub two_qubit: usize, - pub multi_qubit: usize, - pub diagonal: usize, - pub execution_layers: usize, -} - -impl std::fmt::Display for KernelStats { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!( - f, - "Kernels: {} (1q: {}, 2q: {}, 3q+: {}, diag: {}), Layers: {}", - self.total_kernels, - self.single_qubit, - self.two_qubit, - self.multi_qubit, - self.diagonal, - self.execution_layers - ) - } -} diff --git a/libpsi-core/src/core/mod.rs b/libpsi-core/src/core/mod.rs deleted file mode 100644 index c38cf38..0000000 --- a/libpsi-core/src/core/mod.rs +++ /dev/null @@ -1,17 +0,0 @@ -pub mod circuit; -pub mod classical_components; -pub mod custom_gate; -pub mod gates; -pub mod kernel; -pub mod noise; -pub mod quantum_components; -pub mod runtime; - -pub use circuit::*; -pub use classical_components::*; -pub use custom_gate::*; -pub use gates::*; -pub use kernel::*; -pub use noise::*; -pub use quantum_components::*; -pub use runtime::*; diff --git a/libpsi-core/src/core/noise.rs b/libpsi-core/src/core/noise.rs deleted file mode 100644 index 6836c00..0000000 --- a/libpsi-core/src/core/noise.rs +++ /dev/null @@ -1,560 +0,0 @@ -use crate::{complex, Complex, Matrix}; - -#[derive(Clone, Debug)] -pub struct KrausOperator { - pub matrix: Matrix<Complex<f64>>, - pub name: String, -} - -impl KrausOperator { - pub fn new(name: &str, matrix: Matrix<Complex<f64>>) -> Self { - Self { - matrix, - name: name.to_string(), - } - } -} - -#[derive(Clone, Debug)] -pub struct NoiseChannel { - pub name: String, - pub operators: Vec<KrausOperator>, - pub num_qubits: usize, -} - -impl NoiseChannel { - pub fn new(name: &str, operators: Vec<KrausOperator>, num_qubits: usize) -> Self { - Self { - name: name.to_string(), - operators, - num_qubits, - } - } - - pub fn depolarising(p: f64) -> Self { - let sqrt_1_p = (1.0 - p).sqrt(); - let sqrt_p3 = (p / 3.0).sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_1_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(sqrt_p3, 0.0), - complex!(sqrt_p3, 0.0), - complex!(0.0, 0.0), - ], - ); - - let k2 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(0.0, -sqrt_p3), - complex!(0.0, sqrt_p3), - complex!(0.0, 0.0), - ], - ); - - let k3 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_p3, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(-sqrt_p3, 0.0), - ], - ); - - Self::new( - "Depolarising", - vec![ - KrausOperator::new("K0", k0), - KrausOperator::new("K1(X)", k1), - KrausOperator::new("K2(Y)", k2), - KrausOperator::new("K3(Z)", k3), - ], - 1, - ) - } - - pub fn amplitude_damping(gamma: f64) -> Self { - let sqrt_gamma = gamma.sqrt(); - let sqrt_1_gamma = (1.0 - gamma).sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(1.0, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_gamma, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(sqrt_gamma, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - ], - ); - - Self::new( - "AmplitudeDamping", - vec![ - KrausOperator::new("K0", k0), - KrausOperator::new("K1", k1), - ], - 1, - ) - } - - pub fn phase_damping(gamma: f64) -> Self { - let sqrt_gamma = gamma.sqrt(); - let sqrt_1_gamma = (1.0 - gamma).sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(1.0, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_gamma, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_gamma, 0.0), - ], - ); - - Self::new( - "PhaseDamping", - vec![ - KrausOperator::new("K0", k0), - KrausOperator::new("K1", k1), - ], - 1, - ) - } - - pub fn bit_flip(p: f64) -> Self { - let sqrt_1_p = (1.0 - p).sqrt(); - let sqrt_p = p.sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_1_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(sqrt_p, 0.0), - complex!(sqrt_p, 0.0), - complex!(0.0, 0.0), - ], - ); - - Self::new( - "BitFlip", - vec![ - KrausOperator::new("K0(I)", k0), - KrausOperator::new("K1(X)", k1), - ], - 1, - ) - } - - pub fn phase_flip(p: f64) -> Self { - let sqrt_1_p = (1.0 - p).sqrt(); - let sqrt_p = p.sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_1_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(-sqrt_p, 0.0), - ], - ); - - Self::new( - "PhaseFlip", - vec![ - KrausOperator::new("K0(I)", k0), - KrausOperator::new("K1(Z)", k1), - ], - 1, - ) - } - - pub fn bit_phase_flip(p: f64) -> Self { - let sqrt_1_p = (1.0 - p).sqrt(); - let sqrt_p = p.sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_1_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(0.0, -sqrt_p), - complex!(0.0, sqrt_p), - complex!(0.0, 0.0), - ], - ); - - Self::new( - "BitPhaseFlip", - vec![ - KrausOperator::new("K0(I)", k0), - KrausOperator::new("K1(Y)", k1), - ], - 1, - ) - } - - pub fn generalised_amplitude_damping(p: f64, gamma: f64) -> Self { - let sqrt_p = p.sqrt(); - let sqrt_1_p = (1.0 - p).sqrt(); - let sqrt_gamma = gamma.sqrt(); - let sqrt_1_gamma = (1.0 - gamma).sqrt(); - - let k0 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_p, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_p * sqrt_1_gamma, 0.0), - ], - ); - - let k1 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(sqrt_p * sqrt_gamma, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - ], - ); - - let k2 = Matrix::new( - 2, - 2, - vec![ - complex!(sqrt_1_p * sqrt_1_gamma, 0.0), - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p, 0.0), - ], - ); - - let k3 = Matrix::new( - 2, - 2, - vec![ - complex!(0.0, 0.0), - complex!(0.0, 0.0), - complex!(sqrt_1_p * sqrt_gamma, 0.0), - complex!(0.0, 0.0), - ], - ); - - Self::new( - "GeneralisedAmplitudeDamping", - vec![ - KrausOperator::new("K0", k0), - KrausOperator::new("K1", k1), - KrausOperator::new("K2", k2), - KrausOperator::new("K3", k3), - ], - 1, - ) - } -} - -#[derive(Clone)] -pub struct DensityMatrix { - pub data: Vec<Complex<f64>>, - pub dim: usize, - pub num_qubits: usize, -} - -impl DensityMatrix { - pub fn new(num_qubits: usize) -> Self { - let dim = 1 << num_qubits; - let mut data = vec![complex!(0.0, 0.0); dim * dim]; - data[0] = complex!(1.0, 0.0); - Self { - data, - dim, - num_qubits, - } - } - - pub fn from_state_vector(state: &[Complex<f64>]) -> Self { - let dim = state.len(); - let num_qubits = (dim as f64).log2() as usize; - let mut data = vec![complex!(0.0, 0.0); dim * dim]; - - for i in 0..dim { - for j in 0..dim { - data[i * dim + j] = state[i] * state[j].get_conjugate(); - } - } - - Self { - data, - dim, - num_qubits, - } - } - - pub fn get(&self, row: usize, col: usize) -> Complex<f64> { - self.data[row * self.dim + col] - } - - pub fn set(&mut self, row: usize, col: usize, value: Complex<f64>) { - self.data[row * self.dim + col] = value; - } - - pub fn trace(&self) -> Complex<f64> { - let mut sum = complex!(0.0, 0.0); - for i in 0..self.dim { - sum = sum + self.get(i, i); - } - sum - } - - pub fn purity(&self) -> f64 { - let mut sum = complex!(0.0, 0.0); - for i in 0..self.dim { - for j in 0..self.dim { - let rho_ij = self.get(i, j); - let rho_ji = self.get(j, i); - sum = sum + rho_ij * rho_ji; - } - } - sum.real - } - - pub fn is_pure(&self, tolerance: f64) -> bool { - (self.purity() - 1.0).abs() < tolerance - } - - pub fn probabilities(&self) -> Vec<f64> { - (0..self.dim).map(|i| self.get(i, i).real).collect() - } - - pub fn apply_unitary(&mut self, gate: &Matrix<Complex<f64>>, targets: &[usize]) { - let g = targets.len(); - let gate_dim = 1 << g; - - let target_bits: Vec<usize> = targets - .iter() - .map(|&t| self.num_qubits - 1 - t) - .collect(); - - let mut non_target_mask: usize = (1 << self.num_qubits) - 1; - for &pos in &target_bits { - non_target_mask &= !(1 << pos); - } - - let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim]; - - for i in 0..self.dim { - for j in 0..self.dim { - let mut sum = complex!(0.0, 0.0); - - for k in 0..gate_dim { - for l in 0..gate_dim { - let mut src_i = i & non_target_mask; - let mut src_j = j & non_target_mask; - - for (idx, &pos) in target_bits.iter().enumerate() { - if (k >> (g - 1 - idx)) & 1 == 1 { - src_i |= 1 << pos; - } - if (l >> (g - 1 - idx)) & 1 == 1 { - src_j |= 1 << pos; - } - } - - let mut tgt_i = 0usize; - let mut tgt_j = 0usize; - for (idx, &pos) in target_bits.iter().enumerate() { - if (i >> pos) & 1 == 1 { - tgt_i |= 1 << (g - 1 - idx); - } - if (j >> pos) & 1 == 1 { - tgt_j |= 1 << (g - 1 - idx); - } - } - - let u_ik = gate.data[tgt_i * gate_dim + k]; - let u_jl_dag = gate.data[tgt_j * gate_dim + l].get_conjugate(); - let rho_kl = self.get(src_i, src_j); - - sum = sum + u_ik * rho_kl * u_jl_dag; - } - } - - new_data[i * self.dim + j] = sum; - } - } - - self.data = new_data; - } - - pub fn apply_noise_channel(&mut self, channel: &NoiseChannel, target: usize) { - if channel.num_qubits != 1 { - panic!("Only single-qubit noise channels are currently supported"); - } - - let target_bit = self.num_qubits - 1 - target; - let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim]; - - for kraus in &channel.operators { - let k = &kraus.matrix; - - for i in 0..self.dim { - for j in 0..self.dim { - let i_target = (i >> target_bit) & 1; - let j_target = (j >> target_bit) & 1; - - for ki in 0..2 { - for kj in 0..2 { - let src_i = (i & !(1 << target_bit)) | (ki << target_bit); - let src_j = (j & !(1 << target_bit)) | (kj << target_bit); - - let k_elem = k.data[i_target * 2 + ki]; - let k_dag_elem = k.data[j_target * 2 + kj].get_conjugate(); - let rho_elem = self.get(src_i, src_j); - - new_data[i * self.dim + j] = - new_data[i * self.dim + j] + k_elem * rho_elem * k_dag_elem; - } - } - } - } - } - - self.data = new_data; - } - - pub fn measure_probability(&self, qubit: usize, outcome: usize) -> f64 { - let target_bit = self.num_qubits - 1 - qubit; - let mut prob = 0.0; - - for i in 0..self.dim { - if (i >> target_bit) & 1 == outcome { - prob += self.get(i, i).real; - } - } - - prob - } - - pub fn fidelity_with_pure_state(&self, state: &[Complex<f64>]) -> f64 { - let mut sum = complex!(0.0, 0.0); - - for i in 0..self.dim { - for j in 0..self.dim { - sum = sum + state[i].get_conjugate() * self.get(i, j) * state[j]; - } - } - - sum.real - } -} - -impl std::fmt::Display for DensityMatrix { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - writeln!(f, "DensityMatrix ({} qubits, {}×{}):", self.num_qubits, self.dim, self.dim)?; - writeln!(f, " Trace: {:.6}", self.trace().real)?; - writeln!(f, " Purity: {:.6}", self.purity())?; - writeln!(f, " Pure: {}", self.is_pure(1e-10))?; - writeln!(f, " Probabilities: {:?}", self.probabilities())?; - Ok(()) - } -} - -impl std::fmt::Debug for DensityMatrix { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - writeln!(f, "DensityMatrix {}×{}:", self.dim, self.dim)?; - for i in 0..self.dim { - write!(f, " [")?; - for j in 0..self.dim { - let val = self.get(i, j); - if j > 0 { - write!(f, ", ")?; - } - write!(f, "{:.4}+{:.4}i", val.real, val.imaginary)?; - } - writeln!(f, "]")?; - } - Ok(()) - } -} - diff --git a/libpsi-core/src/core/quantum_components.rs b/libpsi-core/src/core/quantum_components.rs deleted file mode 100644 index 75f6325..0000000 --- a/libpsi-core/src/core/quantum_components.rs +++ /dev/null @@ -1,318 +0,0 @@ -use crate::{column_vector, complex, ColumnVector, Complex, Float, Matrix, Vector, VectorMatrix}; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! count { - () => { 0 }; - ($head:expr $(,$tail:expr)*) => { 1 + count!($( $tail ),*) }; -} - -#[macro_export] -macro_rules! qubit { - ($(($re:expr, $im:expr)),*) => { - { - let mut vector = Vec::new(); - $( - vector.push(complex!($re, $im)); - )* - QuantumBit::new(vector) - } - }; -} - -#[macro_export] -macro_rules! quantum_register { - ($($bit:expr),*) => { - { - const N: usize = count!($($bit),*); - let mut bits: [QuantumBit; N] = [$($bit),*]; - QuantumRegister::from(&mut bits) - } - }; -} - -pub type QuantumState = ColumnVector<Complex<f64>>; -impl QuantumState { - pub fn state_0() -> QuantumState { - column_vector![complex!(1.0, 0.0), complex!(0.0, 0.0)] - } - - pub fn state_1() -> QuantumState { - column_vector![complex!(0.0, 0.0), complex!(1.0, 0.0)] - } -} - -fn identity_matrix<T: Float>(size: usize) -> Matrix<T> { - let mut data = vec![T::zero(); size * size]; - for i in 0..size { - data[i * size + i] = T::one(); - } - Matrix::new(size, size, data) -} - -#[derive(Clone)] -pub struct QuantumBit<'a> { - state: QuantumState, - name: &'a str, -} - -#[derive(Clone)] -pub struct QuantumRegister<'a> { - state_vector: QuantumState, - name: &'a str, - qubits: Vec<QuantumBit<'a>>, -} - -#[derive(Clone)] -pub struct QuantumGate<'a> { - pub name: &'a str, - pub matrix: Matrix<Complex<f64>>, - pub num_qubits: usize, -} - -impl<'a> QuantumGate<'a> { - pub fn new(name: &'a str, matrix: Matrix<Complex<f64>>, num_qubits: usize) -> Self { - let expected_dim = 1 << num_qubits; - assert_eq!( - matrix.rows, expected_dim, - "Gate matrix rows must be 2^num_qubits" - ); - assert_eq!( - matrix.cols, expected_dim, - "Gate matrix cols must be 2^num_qubits" - ); - QuantumGate { - name, - matrix, - num_qubits, - } - } - - pub fn from_matrix(name: &'a str, matrix: Matrix<Complex<f64>>) -> Self { - assert_eq!(matrix.rows, matrix.cols, "Gate matrix must be square"); - let dim = matrix.rows; - assert!( - dim > 0 && (dim & (dim - 1)) == 0, - "Matrix dimension must be a power of 2" - ); - let num_qubits = (dim as f64).log2() as usize; - QuantumGate { - name, - matrix, - num_qubits, - } - } -} - -impl<'a> QuantumBit<'a> { - pub fn new(name: &'a str, state: QuantumState) -> QuantumBit<'a> { - QuantumBit { name, state } - } - - pub fn get_state(&self) -> QuantumState { - self.state.clone() - } - - pub fn get_name(&self) -> &'a str { - self.name - } -} - -impl<'a> QuantumRegister<'a> { - pub fn new(name: &'a str, names: &[&'a str]) -> QuantumRegister<'a> { - let mut bits: Vec<QuantumBit<'a>> = Vec::new(); - for i in 0..names.len() { - bits.push(QuantumBit::new(names[i], QuantumState::state_0())) - } - - QuantumRegister::from(name, &mut bits) - } - - pub fn from(name: &'a str, bits: &mut [QuantumBit<'a>]) -> QuantumRegister<'a> { - let mut register = QuantumRegister { - name, - qubits: bits.to_vec(), - state_vector: ColumnVector::new(vec![]), - }; - - register.update(); - register - } - - fn update(&mut self) { - let matrices: Vec<Matrix<Complex<f64>>> = self - .qubits - .iter() - .map(|qubit| qubit.state.to_matrix()) - .collect(); - let mut new_result = matrices[0].clone(); - for matrix in &matrices[1..] { - new_result = new_result.kronecker(matrix); - } - - self.state_vector = ColumnVector::from_matrix(&new_result); - } - - pub fn get_bits(&self) -> Vec<QuantumBit<'_>> { - self.qubits.clone() - } - - pub fn get_state(&self) -> QuantumState { - self.state_vector.clone() - } - - pub fn get_name(&self) -> &'a str { - self.name - } - - pub fn num_qubits(&self) -> usize { - self.qubits.len() - } - - pub fn apply_gate(&mut self, gate: &QuantumGate, targets: &[usize]) { - let n = self.num_qubits(); - - assert_eq!( - gate.num_qubits, - targets.len(), - "Number of target qubits must match gate's qubit count" - ); - for &t in targets { - assert!( - t < n, - "Target qubit index {} out of range for {}-qubit register", - t, - n - ); - } - - let mut sorted_targets = targets.to_vec(); - sorted_targets.sort(); - for i in 1..sorted_targets.len() { - assert_ne!( - sorted_targets[i], - sorted_targets[i - 1], - "Duplicate target qubit indices are not allowed" - ); - } - - let full_operator = self.build_full_operator(gate, targets); - - self.state_vector = self - .state_vector - .mul_matrix(&full_operator) - .expect("Matrix multiplication failed during gate application"); - } - - fn build_full_operator(&self, gate: &QuantumGate, targets: &[usize]) -> Matrix<Complex<f64>> { - let n = self.num_qubits(); - let g = gate.num_qubits; - let dim = 1 << n; - - let mut contiguous = true; - for i in 1..targets.len() { - if targets[i] != targets[i - 1] + 1 { - contiguous = false; - break; - } - } - - if contiguous && g == n { - return gate.matrix.clone(); - } - - if contiguous { - return self.build_contiguous_operator(gate, targets[0]); - } - - let mut result = Matrix::new(dim, dim, vec![complex!(0.0, 0.0); dim * dim]); - - for col in 0..dim { - for row in 0..dim { - let mut target_row_bits = 0usize; - let mut target_col_bits = 0usize; - - for (i, &t) in targets.iter().enumerate() { - let qubit_pos = n - 1 - t; - if (row >> qubit_pos) & 1 == 1 { - target_row_bits |= 1 << (g - 1 - i); - } - if (col >> qubit_pos) & 1 == 1 { - target_col_bits |= 1 << (g - 1 - i); - } - } - - let mut non_target_match = true; - for q in 0..n { - if !targets.contains(&q) { - let qubit_pos = n - 1 - q; - if ((row >> qubit_pos) & 1) != ((col >> qubit_pos) & 1) { - non_target_match = false; - break; - } - } - } - - if non_target_match { - result.set(row, col, gate.matrix.get(target_row_bits, target_col_bits)); - } - } - } - - result - } - - fn build_contiguous_operator( - &self, - gate: &QuantumGate, - start_idx: usize, - ) -> Matrix<Complex<f64>> { - let n = self.num_qubits(); - let g = gate.num_qubits; - - let mut result: Option<Matrix<Complex<f64>>> = None; - - for i in 0..n { - let part: Matrix<Complex<f64>> = if i == start_idx { - gate.matrix.clone() - } else if i > start_idx && i < start_idx + g { - continue; - } else { - identity_matrix(2) - }; - - result = Some(match result { - None => part, - Some(r) => r.kronecker(&part), - }); - } - - result.unwrap_or_else(|| identity_matrix(1 << n)) - } - - pub fn apply_gates(&mut self, operations: &[(&QuantumGate, &[usize])]) { - for (gate, targets) in operations { - self.apply_gate(gate, targets); - } - } -} - -impl<'a> ops::Index<usize> for QuantumRegister<'a> { - type Output = QuantumBit<'a>; - - fn index(&self, index: usize) -> &Self::Output { - &self.qubits[index] - } -} - -impl<'a> ops::IndexMut<usize> for QuantumRegister<'a> { - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.qubits[index] - } -} - -impl<'a> fmt::Display for QuantumGate<'a> { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "{}", self.name) - } -} diff --git a/libpsi-core/src/core/runtime.rs b/libpsi-core/src/core/runtime.rs deleted file mode 100644 index 382fa2f..0000000 --- a/libpsi-core/src/core/runtime.rs +++ /dev/null @@ -1,585 +0,0 @@ -use super::{ - GateOp, Kernel, KernelBatch, QuantumGate, QuantumRegister, QuantumState, - StructureAwareKernelBatch, -}; -use crate::gates::{ - cp_matrix, crx_matrix, cry_matrix, crz_matrix, p_matrix, rx_matrix, ry_matrix, rz_matrix, - u1_matrix, u2_matrix, u3_matrix, CNOT, CZ, FREDKIN, HADAMARD, PAULI_X, PAULI_Y, PAULI_Z, - SDG_GATE, SWAP, SXDG_GATE, SX_GATE, S_GATE, TDG_GATE, TOFFOLI, T_GATE, -}; -use crate::maths::simd::{apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel}; -use crate::maths::vector::Vector; -use crate::{complex, Complex, Matrix}; -use rayon::prelude::*; - -const PARALLEL_THRESHOLD: usize = 8; - -#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] -pub struct RuntimeConfig { - pub parallel: bool, - pub simd: bool, - pub batched: bool, - pub structure_aware: bool, - pub parallel_threshold: usize, -} - -impl RuntimeConfig { - pub fn new() -> Self { - Self { - parallel: false, - simd: false, - batched: false, - structure_aware: false, - parallel_threshold: PARALLEL_THRESHOLD, - } - } - - pub fn parallel(mut self) -> Self { - self.parallel = true; - self - } - - pub fn simd(mut self) -> Self { - self.simd = true; - self - } - - pub fn batched(mut self) -> Self { - self.batched = true; - self - } - - pub fn structure_aware(mut self) -> Self { - self.structure_aware = true; - self - } - - pub fn with_threshold(mut self, threshold: usize) -> Self { - self.parallel_threshold = threshold; - self - } - - pub fn optimal() -> Self { - Self::new().structure_aware().simd().parallel() - } - - pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState { - let dim = 1 << num_qubits; - let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim]; - state[0] = complex!(1.0, 0.0); - - let use_parallel = self.parallel && num_qubits >= self.parallel_threshold; - - if self.structure_aware { - let mut batch = Runtime::build_structure_aware_batch(num_qubits, operations); - batch.optimise(); - self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel); - } else if self.batched { - let mut batch = Runtime::build_kernel_batch(num_qubits, operations); - batch.optimize(); - self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel); - } else { - let batch = Runtime::build_kernel_batch(num_qubits, operations); - self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel); - } - - QuantumState::new(state) - } - - fn execute_kernels( - &self, - state: &mut Vec<Complex<f64>>, - kernels: &[Kernel], - num_qubits: usize, - use_parallel: bool, - ) { - for kernel in kernels { - if self.simd && kernel.targets.len() == 1 { - let gate = matrix_to_2x2(&kernel.matrix); - if use_parallel { - apply_single_qubit_gate_simd_parallel( - state, - &gate, - kernel.targets[0], - num_qubits, - ); - } else { - apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], num_qubits); - } - } else if use_parallel { - *state = apply_gate_parallel(state, &kernel.matrix, &kernel.targets, num_qubits); - } else { - *state = apply_kernel_direct(state, kernel, num_qubits); - } - } - } -} - -impl std::fmt::Display for RuntimeConfig { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - let mut features = Vec::new(); - if self.structure_aware { - features.push("structure-aware"); - } - if self.batched && !self.structure_aware { - features.push("batched"); - } - if self.simd { - features.push("SIMD"); - } - if self.parallel { - features.push("parallel"); - } - if features.is_empty() { - features.push("basic"); - } - write!(f, "Runtime[{}]", features.join("+")) - } -} - -#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] -pub enum Runtime { - #[default] - BasicRT, - BasicRTMT, - BatchedRT, - BatchedRTMT, - SimdRT, - SimdRTMT, - StructureAwareRT, - StructureAwareMT, - WFEvolution, - WFEvolutionMT, - GPUAccelerated, - Custom(RuntimeConfig), -} - -impl Runtime { - pub fn custom() -> RuntimeConfig { - RuntimeConfig::new() - } - - pub fn optimal() -> RuntimeConfig { - RuntimeConfig::optimal() - } - - pub fn to_config(&self) -> RuntimeConfig { - match self { - Runtime::BasicRT => RuntimeConfig::new(), - Runtime::BasicRTMT => RuntimeConfig::new().parallel(), - Runtime::BatchedRT => RuntimeConfig::new().batched(), - Runtime::BatchedRTMT => RuntimeConfig::new().batched().parallel(), - Runtime::SimdRT => RuntimeConfig::new().batched().simd(), - Runtime::SimdRTMT => RuntimeConfig::new().batched().simd().parallel(), - Runtime::StructureAwareRT => RuntimeConfig::new().structure_aware().simd(), - Runtime::StructureAwareMT => RuntimeConfig::new().structure_aware().simd().parallel(), - Runtime::Custom(config) => *config, - _ => RuntimeConfig::new(), - } - } - - pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState { - match self { - Runtime::BasicRT => Self::compute_basic(num_qubits, operations), - Runtime::BasicRTMT => Self::compute_basic_mt(num_qubits, operations), - Runtime::Custom(config) => config.compute(num_qubits, operations), - Runtime::WFEvolution => { - unimplemented!("WFEvolution (Schrödinger equation) runtime not yet implemented") - } - Runtime::WFEvolutionMT => { - unimplemented!( - "WFEvolutionMT (multi-threaded Schrödinger) runtime not yet implemented" - ) - } - Runtime::GPUAccelerated => { - unimplemented!("GPUAccelerated runtime not yet implemented") - } - _ => self.to_config().compute(num_qubits, operations), - } - } - - pub fn build_kernel_batch(num_qubits: usize, operations: &[GateOp]) -> KernelBatch { - let mut batch = KernelBatch::new(num_qubits); - - for op in operations { - if let Some(kernel) = Self::op_to_kernel(op) { - batch.add(kernel); - } - } - - batch - } - - fn op_to_kernel(op: &GateOp) -> Option<Kernel> { - let (matrix, targets, name): (Matrix<Complex<f64>>, Vec<usize>, &str) = match op { - GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t], "H"), - GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t], "X"), - GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t], "Y"), - GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t], "Z"), - GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t], "S"), - GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t], "T"), - GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t], "Sdg"), - GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t], "Tdg"), - GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t], "Sx"), - GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t], "Sxdg"), - GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t], "Rx"), - GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t], "Ry"), - GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t], "Rz"), - GateOp::P(t, theta) => (p_matrix(*theta), vec![*t], "P"), - GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t], "U1"), - GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t], "U2"), - GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t], "U3"), - GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t], "CNOT"), - GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t], "CZ"), - GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b], "SWAP"), - GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t], "CRx"), - GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t], "CRy"), - GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t], "CRz"), - GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t], "CP"), - GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t], "CCNOT"), - GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2], "CSWAP"), - GateOp::Measure(_, _) => return None, - GateOp::Custom(gate, tgts) => { - let qg = gate.to_quantum_gate(); - (qg.matrix, tgts.clone(), "Custom") - } - }; - - Some(Kernel::new(name, matrix, targets)) - } - - pub fn build_structure_aware_batch( - num_qubits: usize, - operations: &[GateOp], - ) -> StructureAwareKernelBatch { - let mut batch = StructureAwareKernelBatch::new(num_qubits); - - for op in operations { - if let Some(kernel) = Self::op_to_kernel(op) { - batch.add(kernel); - } - } - - batch - } - - fn compute_basic(num_qubits: usize, operations: &[GateOp]) -> QuantumState { - let names: Vec<String> = (0..num_qubits).map(|i| format!("q{}", i)).collect(); - let leaked_names: &'static [String] = Box::leak(names.into_boxed_slice()); - let name_refs: Vec<&'static str> = leaked_names.iter().map(|s| s.as_str()).collect(); - - let mut register = QuantumRegister::new( - Box::leak(Box::new("circuit".to_string())).as_str(), - &name_refs, - ); - - for op in operations { - match op { - // Clifford gates - GateOp::H(t) => register.apply_gate(&HADAMARD, &[*t]), - GateOp::X(t) => register.apply_gate(&PAULI_X, &[*t]), - GateOp::Y(t) => register.apply_gate(&PAULI_Y, &[*t]), - GateOp::Z(t) => register.apply_gate(&PAULI_Z, &[*t]), - GateOp::S(t) => register.apply_gate(&S_GATE, &[*t]), - GateOp::CNOT(c, t) => register.apply_gate(&CNOT, &[*c, *t]), - GateOp::CZ(c, t) => register.apply_gate(&CZ, &[*c, *t]), - GateOp::SWAP(a, b) => register.apply_gate(&SWAP, &[*a, *b]), - GateOp::CCNOT(c1, c2, t) => register.apply_gate(&TOFFOLI, &[*c1, *c2, *t]), - GateOp::CSWAP(c, t1, t2) => register.apply_gate(&FREDKIN, &[*c, *t1, *t2]), - - // Non-Clifford fixed gates - GateOp::T(t) => register.apply_gate(&T_GATE, &[*t]), - GateOp::Sdg(t) => register.apply_gate(&SDG_GATE, &[*t]), - GateOp::Tdg(t) => register.apply_gate(&TDG_GATE, &[*t]), - GateOp::Sx(t) => register.apply_gate(&SX_GATE, &[*t]), - GateOp::Sxdg(t) => register.apply_gate(&SXDG_GATE, &[*t]), - - // Parametric single-qubit gates (non-Clifford for most angles) - GateOp::Rx(t, theta) => { - let gate = QuantumGate { - name: "Rx", - matrix: rx_matrix(*theta), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::Ry(t, theta) => { - let gate = QuantumGate { - name: "Ry", - matrix: ry_matrix(*theta), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::Rz(t, theta) => { - let gate = QuantumGate { - name: "Rz", - matrix: rz_matrix(*theta), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::P(t, theta) => { - let gate = QuantumGate { - name: "P", - matrix: p_matrix(*theta), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::U1(t, lambda) => { - let gate = QuantumGate { - name: "U1", - matrix: u1_matrix(*lambda), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::U2(t, phi, lambda) => { - let gate = QuantumGate { - name: "U2", - matrix: u2_matrix(*phi, *lambda), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - GateOp::U3(t, theta, phi, lambda) => { - let gate = QuantumGate { - name: "U3", - matrix: u3_matrix(*theta, *phi, *lambda), - num_qubits: 1, - }; - register.apply_gate(&gate, &[*t]); - } - - // Controlled parametric gates - GateOp::CRx(c, t, theta) => { - let gate = QuantumGate { - name: "CRx", - matrix: crx_matrix(*theta), - num_qubits: 2, - }; - register.apply_gate(&gate, &[*c, *t]); - } - GateOp::CRy(c, t, theta) => { - let gate = QuantumGate { - name: "CRy", - matrix: cry_matrix(*theta), - num_qubits: 2, - }; - register.apply_gate(&gate, &[*c, *t]); - } - GateOp::CRz(c, t, theta) => { - let gate = QuantumGate { - name: "CRz", - matrix: crz_matrix(*theta), - num_qubits: 2, - }; - register.apply_gate(&gate, &[*c, *t]); - } - GateOp::CP(c, t, theta) => { - let gate = QuantumGate { - name: "CP", - matrix: cp_matrix(*theta), - num_qubits: 2, - }; - register.apply_gate(&gate, &[*c, *t]); - } - - // Measurement and custom gates - GateOp::Measure(_, _) => {} - GateOp::Custom(gate, targets) => { - let quantum_gate = gate.to_quantum_gate(); - register.apply_gate(&quantum_gate, targets); - } - } - } - - register.get_state() - } - - fn compute_basic_mt(num_qubits: usize, operations: &[GateOp]) -> QuantumState { - // For small circuits, fall back to single-threaded (overhead not worth it) - if num_qubits < PARALLEL_THRESHOLD { - return Self::compute_basic(num_qubits, operations); - } - - let dim = 1 << num_qubits; - - // Initialize state to |0...0⟩ - let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim]; - state[0] = complex!(1.0, 0.0); - - for op in operations { - let (gate_matrix, targets): (Matrix<Complex<f64>>, Vec<usize>) = match op { - // Clifford gates - GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t]), - GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t]), - GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t]), - GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t]), - GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t]), - GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t]), - GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t]), - GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b]), - GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t]), - GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2]), - - // Non-Clifford fixed gates - GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t]), - GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t]), - GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t]), - GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t]), - GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t]), - - // Parametric single-qubit gates - GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t]), - GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t]), - GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t]), - GateOp::P(t, theta) => (p_matrix(*theta), vec![*t]), - GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t]), - GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t]), - GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t]), - - // Controlled parametric gates - GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t]), - GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t]), - GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t]), - GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t]), - - // Measurement (skip) and custom gates - GateOp::Measure(_, _) => continue, - GateOp::Custom(custom_gate, tgts) => { - let quantum_gate = custom_gate.to_quantum_gate(); - state = apply_gate_parallel(&state, &quantum_gate.matrix, tgts, num_qubits); - continue; - } - }; - - state = apply_gate_parallel(&state, &gate_matrix, &targets, num_qubits); - } - - QuantumState::new(state) - } -} - -/// Apply a gate to the state vector in parallel using sparse application -/// This is O(2^n * 2^g) instead of O(2^2n) for full matrix multiplication -fn apply_gate_parallel( - state: &[Complex<f64>], - gate_matrix: &Matrix<Complex<f64>>, - targets: &[usize], - num_qubits: usize, -) -> Vec<Complex<f64>> { - let dim = 1 << num_qubits; - let g = targets.len(); - let gate_dim = 1 << g; - - // Convert target qubit indices to bit positions (from MSB) - let target_bits: Vec<usize> = targets.iter().map(|&t| num_qubits - 1 - t).collect(); - - // Create a mask for non-target qubits - let mut non_target_mask: usize = (1 << num_qubits) - 1; - for &pos in &target_bits { - non_target_mask &= !(1 << pos); - } - - // Parallel computation of new state - let new_state: Vec<Complex<f64>> = (0..dim) - .into_par_iter() - .map(|i| { - // Extract the target qubit bits from index i - let mut target_idx = 0usize; - for (k, &pos) in target_bits.iter().enumerate() { - if (i >> pos) & 1 == 1 { - target_idx |= 1 << (g - 1 - k); - } - } - - // Compute the contribution to state[i] - let mut sum = complex!(0.0, 0.0); - - // For each possible input state that could contribute - for j in 0..gate_dim { - // Get the gate matrix element - let gate_elem = gate_matrix.data[target_idx * gate_dim + j]; - - // Skip if zero (sparse optimization) - if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 { - continue; - } - - // Compute the source index by replacing target bits in i with bits from j - let mut source_idx = i & non_target_mask; - for (k, &pos) in target_bits.iter().enumerate() { - if (j >> (g - 1 - k)) & 1 == 1 { - source_idx |= 1 << pos; - } - } - - sum = sum + gate_elem * state[source_idx]; - } - - sum - }) - .collect(); - - new_state -} - -fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] { - [ - [matrix.data[0], matrix.data[1]], - [matrix.data[2], matrix.data[3]], - ] -} - -fn apply_kernel_direct( - state: &[Complex<f64>], - kernel: &Kernel, - num_qubits: usize, -) -> Vec<Complex<f64>> { - let dim = 1 << num_qubits; - let g = kernel.targets.len(); - let gate_dim = 1 << g; - - let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect(); - - let mut non_target_mask: usize = (1 << num_qubits) - 1; - for &pos in &target_bits { - non_target_mask &= !(1 << pos); - } - - let mut new_state = vec![complex!(0.0, 0.0); dim]; - - for i in 0..dim { - let mut target_idx = 0usize; - for (k, &pos) in target_bits.iter().enumerate() { - if (i >> pos) & 1 == 1 { - target_idx |= 1 << (g - 1 - k); - } - } - - let mut sum = complex!(0.0, 0.0); - - for j in 0..gate_dim { - let gate_elem = kernel.matrix.data[target_idx * gate_dim + j]; - - if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 { - continue; - } - - let mut source_idx = i & non_target_mask; - for (k, &pos) in target_bits.iter().enumerate() { - if (j >> (g - 1 - k)) & 1 == 1 { - source_idx |= 1 << pos; - } - } - - sum = sum + gate_elem * state[source_idx]; - } - - new_state[i] = sum; - } - - new_state -} diff --git a/libpsi-core/src/lib.rs b/libpsi-core/src/lib.rs deleted file mode 100644 index 5876918..0000000 --- a/libpsi-core/src/lib.rs +++ /dev/null @@ -1,18 +0,0 @@ -pub mod core; -pub mod maths; - -pub use maths::complex::*; -pub use maths::format::*; -pub use maths::matrix::*; -pub use maths::numeric::*; -pub use maths::simd::*; -pub use maths::vector::*; - -pub use core::circuit::*; -pub use core::classical_components::*; -pub use core::custom_gate::*; -pub use core::gates; -pub use core::kernel::*; -pub use core::noise::*; -pub use core::quantum_components::*; -pub use core::runtime::*; diff --git a/libpsi-core/src/maths/complex.rs b/libpsi-core/src/maths/complex.rs deleted file mode 100644 index 31eae69..0000000 --- a/libpsi-core/src/maths/complex.rs +++ /dev/null @@ -1,180 +0,0 @@ -use crate::Float; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! complex { - ($real:expr, $imaginary:expr) => { - $crate::Complex::new($real, $imaginary) - }; -} - -macro_rules! impl_ops { - ($trait:ident, $method:ident, $op:tt) => { - impl<T: Float> ops::$trait for Complex<T> { - type Output = Complex<T>; - - fn $method(self, other: Complex<T>) -> Complex<T> { - Complex { - real: self.real $op other.real, - imaginary: self.imaginary $op other.imaginary, - } - } - } - }; - - ($trait:ident, $method:ident, $op:tt, real) => { - impl<T: Float> ops::$trait<T> for Complex<T> { - type Output = Complex<T>; - - fn $method(self, other: T) -> Complex<T> { - Complex { - real: self.real $op other, - imaginary: self.imaginary, - } - } - } - }; - - ($trait_assign:ident, $method_assign:ident, $op:tt, assign) => { - impl<T: Float> ops::$trait_assign for Complex<T> { - fn $method_assign(&mut self, other: Complex<T>) { - self.real = self.real $op other.real; - self.imaginary = self.imaginary $op other.imaginary; - } - } - }; - - ($trait_assign:ident, $method_assign:ident, $op:tt, assign_real) => { - impl<T: Float> ops::$trait_assign<T> for Complex<T> { - fn $method_assign(&mut self, other: T) { - self.real = self.real $op other; - } - } - }; -} - -#[derive(Copy, Clone, PartialOrd, PartialEq)] -pub struct Complex<T: Float> { - pub real: T, - pub imaginary: T, -} - -impl<T: Float + fmt::Debug> fmt::Debug for Complex<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "Complex {{ real: {:?}, imaginary: {:?} }}", - self.real, self.imaginary - ) - } -} - -impl<T: Float + fmt::Display> fmt::Display for Complex<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "{} + {}i", self.real, self.imaginary) - } -} - -impl<T: Float> ops::Neg for Complex<T> { - type Output = Complex<T>; - - fn neg(self) -> Complex<T> { - Complex { - real: -self.real, - imaginary: -self.imaginary, - } - } -} - -impl<T: Float> From<T> for Complex<T> { - fn from(real: T) -> Complex<T> { - Complex { - real, - imaginary: T::zero(), - } - } -} - -impl<T: Float> Complex<T> { - pub fn new(real: T, imaginary: T) -> Complex<T> { - Complex { real, imaginary } - } - - pub fn get_conjugate(&self) -> Complex<T> { - Complex { - real: self.real, - imaginary: -self.imaginary, - } - } - - pub fn conjugate(&mut self) { - self.imaginary = -self.imaginary; - } - - pub fn phase(&self) -> T { - T::atan2(self.imaginary, self.real) - } - - pub fn norm2(&self) -> T { - self.real * self.real + self.imaginary * self.imaginary - } - - pub fn abs(&self) -> T { - T::sqrt(self.norm2()) - } -} - -impl_ops!(Add, add, +); -impl_ops!(Sub, sub, -); - -impl<T: Float> ops::Mul for Complex<T> { - type Output = Complex<T>; - - fn mul(self, other: Complex<T>) -> Complex<T> { - // (a + bi) * (c + di) = (ac - bd) + (ad + bc)i - Complex { - real: self.real * other.real - self.imaginary * other.imaginary, - imaginary: self.real * other.imaginary + self.imaginary * other.real, - } - } -} - -impl<T: Float> ops::Div for Complex<T> { - type Output = Complex<T>; - - fn div(self, other: Complex<T>) -> Complex<T> { - // (a + bi) / (c + di) = ((ac + bd) + (bc - ad)i) / (c² + d²) - let denom = other.real * other.real + other.imaginary * other.imaginary; - Complex { - real: (self.real * other.real + self.imaginary * other.imaginary) / denom, - imaginary: (self.imaginary * other.real - self.real * other.imaginary) / denom, - } - } -} - -impl_ops!(AddAssign, add_assign, +, assign); -impl_ops!(SubAssign, sub_assign, -, assign); - -impl<T: Float> ops::MulAssign for Complex<T> { - fn mul_assign(&mut self, other: Complex<T>) { - let new_real = self.real * other.real - self.imaginary * other.imaginary; - let new_imag = self.real * other.imaginary + self.imaginary * other.real; - self.real = new_real; - self.imaginary = new_imag; - } -} - -impl<T: Float> ops::DivAssign for Complex<T> { - fn div_assign(&mut self, other: Complex<T>) { - let denom = other.real * other.real + other.imaginary * other.imaginary; - let new_real = (self.real * other.real + self.imaginary * other.imaginary) / denom; - let new_imag = (self.imaginary * other.real - self.real * other.imaginary) / denom; - self.real = new_real; - self.imaginary = new_imag; - } -} - -impl_ops!(Add, add, +, real); -impl_ops!(Sub, sub, -, real); -impl_ops!(Mul, mul, *, real); -impl_ops!(Div, div, /, real); diff --git a/libpsi-core/src/maths/format.rs b/libpsi-core/src/maths/format.rs deleted file mode 100644 index 957f1b6..0000000 --- a/libpsi-core/src/maths/format.rs +++ /dev/null @@ -1,141 +0,0 @@ -use crate::Complex; - -const EPSILON: f64 = 1e-10; -const SQRT_2: f64 = 1.4142135623730951; -const INV_SQRT_2: f64 = 0.7071067811865475; -const INV_SQRT_8: f64 = 0.3535533905932738; -const INV_SQRT_32: f64 = 0.1767766952966369; - -fn approx_eq(a: f64, b: f64) -> bool { - (a - b).abs() < EPSILON -} - -fn format_real_symbolic(v: f64) -> Option<String> { - let abs_v = v.abs(); - let sign = if v < 0.0 { "-" } else { "" }; - - if approx_eq(abs_v, 0.0) { - return Some("0".to_string()); - } - if approx_eq(abs_v, 1.0) { - return Some(format!("{}1", sign)); - } - if approx_eq(abs_v, 0.5) { - return Some(format!("{}½", sign)); - } - if approx_eq(abs_v, 0.25) { - return Some(format!("{}¼", sign)); - } - if approx_eq(abs_v, 0.75) { - return Some(format!("{}¾", sign)); - } - if approx_eq(abs_v, 0.125) { - return Some(format!("{}⅛", sign)); - } - if approx_eq(abs_v, SQRT_2) { - return Some(format!("{}√2", sign)); - } - if approx_eq(abs_v, INV_SQRT_2) { - return Some(format!("{}¹⁄√2", sign)); - } - if approx_eq(abs_v, INV_SQRT_8) { - return Some(format!("{}¹⁄√8", sign)); - } - if approx_eq(abs_v, INV_SQRT_32) { - return Some(format!("{}¹⁄√32", sign)); - } - if approx_eq(abs_v, 2.0) { - return Some(format!("{}2", sign)); - } - if approx_eq(abs_v, 1.0 / 3.0) { - return Some(format!("{}⅓", sign)); - } - if approx_eq(abs_v, 2.0 / 3.0) { - return Some(format!("{}⅔", sign)); - } - - None -} - -pub fn format_amplitude(c: &Complex<f64>) -> String { - let re = c.real; - let im = c.imaginary; - - let re_zero = approx_eq(re.abs(), 0.0); - let im_zero = approx_eq(im.abs(), 0.0); - - if re_zero && im_zero { - return "0".to_string(); - } - - if im_zero { - if let Some(s) = format_real_symbolic(re) { - return s; - } - return format!("{:.4}", re); - } - - if re_zero { - if approx_eq(im.abs(), 1.0) { - return if im > 0.0 { - "i".to_string() - } else { - "-i".to_string() - }; - } - if let Some(s) = format_real_symbolic(im) { - return format!("{}i", s); - } - return format!("{:.4}i", im); - } - - let re_str = format_real_symbolic(re).unwrap_or_else(|| format!("{:.4}", re)); - let im_str = if approx_eq(im.abs(), 1.0) { - if im > 0.0 { - "+i".to_string() - } else { - "-i".to_string() - } - } else { - let im_sym = format_real_symbolic(im.abs()); - let sign = if im > 0.0 { "+" } else { "-" }; - match im_sym { - Some(s) => format!("{}{}i", sign, s.trim_start_matches('-')), - None => format!("{}{:.4}i", sign, im.abs()), - } - }; - - format!("{}{}", re_str, im_str) -} - -pub fn format_probability(p: f64) -> String { - if approx_eq(p, 0.0) { - return "0".to_string(); - } - if approx_eq(p, 1.0) { - return "1".to_string(); - } - if approx_eq(p, 0.5) { - return "½".to_string(); - } - if approx_eq(p, 0.25) { - return "¼".to_string(); - } - if approx_eq(p, 0.75) { - return "¾".to_string(); - } - if approx_eq(p, 0.125) { - return "⅛".to_string(); - } - if approx_eq(p, 0.0625) { - return "¹⁄₁₆".to_string(); - } - if approx_eq(p, 1.0 / 3.0) { - return "⅓".to_string(); - } - if approx_eq(p, 2.0 / 3.0) { - return "⅔".to_string(); - } - - format!("{:.4}", p) -} diff --git a/libpsi-core/src/maths/matrix.rs b/libpsi-core/src/maths/matrix.rs deleted file mode 100644 index 58492f9..0000000 --- a/libpsi-core/src/maths/matrix.rs +++ /dev/null @@ -1,310 +0,0 @@ -use super::Float; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! matrix { - ( $( $( $x:expr ),* );* ) => {{ - let mut data = Vec::new(); - let mut rows = 0; - let mut cols = 0; - - $( - let row_data = $( $x )*; - if cols == 0 { - cols = row_data.len(); - } - assert_eq!(cols, row_data.len(), "All rows must have the same number of columns."); - data.extend(row_data); - rows += 1; - )* - - $crate::Matrix::new(rows, cols, data) - }}; -} - -macro_rules! impl_matrix_ops { - ($($trait:ident, $method:ident, $other:ty, $output:ty, $scale_fn:ident),* $(,)?) => { - $( - impl<T: Float> core::ops::$trait<$other> for Matrix<T> { - type Output = $output; - - fn $method(self, other: $other) -> Self::Output { - self.$scale_fn(other) - } - } - )* - }; - ($($trait:ident, $method:ident, $other:ty, $scale_fn:ident),* $(,)?) => { - $( - impl<T: Float> core::ops::$trait<$other> for Matrix<T> { - fn $method(&mut self, other: $other) { - *self = self.$scale_fn(other); - } - } - )* - }; -} - -#[derive(Clone)] -pub struct Matrix<T: Float> { - pub data: Vec<T>, - pub rows: usize, - pub cols: usize, -} - -impl<T: Float> Matrix<T> { - pub fn new(rows: usize, cols: usize, data: Vec<T>) -> Self { - Matrix { data, rows, cols } - } - - pub fn get(&self, row: usize, col: usize) -> T { - self.data[row * self.cols + col] - } - - pub fn set(&mut self, row: usize, col: usize, value: T) { - self.data[row * self.cols + col] = value; - } - - pub fn dot(&self, other: &Self) -> Option<Matrix<T>> { - if self.cols != other.rows { - return None; - } - - let mut result = Matrix::new( - self.rows, - other.cols, - vec![T::zero(); self.rows * other.cols], - ); - for i in 0..self.rows { - for j in 0..other.cols { - let mut sum = T::zero(); - for k in 0..self.cols { - sum = sum + (self.get(i, k) * other.get(k, j)); - } - result.set(i, j, sum); - } - } - Some(result) - } - - pub fn kronecker(&self, other: &Self) -> Matrix<T> { - let new_rows = self.rows * other.rows; - let new_cols = self.cols * other.cols; - - let mut result = Matrix::new(new_rows, new_cols, vec![T::zero(); new_rows * new_cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let self_val = self.get(i, j); - for k in 0..other.rows { - for l in 0..other.cols { - let result_row = i * other.rows + k; - let result_col = j * other.cols + l; - result.set(result_row, result_col, self_val.clone() * other.get(k, l)); - } - } - } - } - - result - } - - pub fn transpose(&self) -> Matrix<T> { - let mut result = Matrix::new(self.cols, self.rows, vec![T::zero(); self.cols * self.rows]); - - for i in 0..self.rows { - for j in 0..self.cols { - let value = self.get(i, j); - result.set(j, i, value); - } - } - - result - } - - pub fn add_to(&self, other: &Self) -> Option<Matrix<T>> { - if self.rows != other.rows || self.cols != other.cols { - return None; - } - - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let sum = self.get(i, j) + other.get(i, j); - result.set(i, j, sum); - } - } - Some(result) - } - - pub fn subtract(&self, other: &Self) -> Option<Matrix<T>> { - if self.rows != other.rows || self.cols != other.cols { - return None; - } - - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let diff = self.get(i, j) - other.get(i, j); - result.set(i, j, diff); - } - } - Some(result) - } - - pub fn scale(&self, scalar: T) -> Matrix<T> { - let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]); - - for i in 0..self.rows { - for j in 0..self.cols { - let scaled_value = self.get(i, j) * scalar; - result.set(i, j, scaled_value); - } - } - result - } -} - -impl<T: Float> ops::Index<(usize, usize)> for Matrix<T> { - type Output = T; - - fn index(&self, index: (usize, usize)) -> &Self::Output { - &self.data[index.0 * self.cols + index.1] - } -} - -impl<T: Float> ops::IndexMut<(usize, usize)> for Matrix<T> { - fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output { - &mut self.data[index.0 * self.cols + index.1] - } -} - -impl<T: Float> ops::AddAssign<&Matrix<T>> for Matrix<T> { - fn add_assign(&mut self, other: &Matrix<T>) { - if let Some(result) = self.add_to(other) { - *self = result; - } - } -} - -impl<T: Float> ops::SubAssign<&Matrix<T>> for Matrix<T> { - fn sub_assign(&mut self, other: &Matrix<T>) { - if let Some(result) = self.subtract(other) { - *self = result; - } - } -} - -impl_matrix_ops! { - Add, add, &Matrix<T>, Option<Matrix<T>>, add_to, - Sub, sub, &Matrix<T>, Option<Matrix<T>>, subtract, - Mul, mul, T, Matrix<T>, scale, - Div, div, T, Matrix<T>, scale, -} - -impl_matrix_ops! { - MulAssign, mul_assign, T, scale, - DivAssign, div_assign, T, scale, -} - -impl<T: Float + fmt::Debug> fmt::Debug for Matrix<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - for i in 0..self.rows { - for j in 0..self.cols { - write!(f, "{:?} ", self.get(i, j))?; - } - writeln!(f)?; - } - Ok(()) - } -} - -impl<T: Float + fmt::Display> fmt::Display for Matrix<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let elements: Vec<String> = self.data.iter().map(ToString::to_string).collect(); - let is_complex = elements.iter().any(|element| element.contains("i")); - - let normalized: Vec<(f64, f64)> = self - .data - .iter() - .map(|element| { - let element_string = element.to_string(); - - if is_complex { - let element_string = element_string.trim_end_matches('i').trim(); - let element_split: Vec<&str> = element_string.split_whitespace().collect(); - let real = element_split[0].parse::<f64>().unwrap(); - let imaginary = element_split - .get(2) - .map_or(0.0, |&s| s.parse::<f64>().unwrap()); - (real, imaginary) - } else { - (element_string.parse::<f64>().unwrap(), 0.0) - } - }) - .collect(); - - let max_widths = normalized - .iter() - .fold((0, 0), |(max_0, max_1), &(real, imag)| { - let new_max_0 = max_0.max(format!("{:.2}", real).len()); - let new_max_1 = if is_complex { - max_1.max(format!("{:.2}", imag.abs()).len()) - } else { - max_1 - }; - (new_max_0, new_max_1) - }); - - let aligned: Vec<String> = normalized - .iter() - .map(|&(real, imag)| { - if is_complex { - format!( - "{:>rewidth$.2} {} {:>imwidth$.2}i", - real, - if imag > 0.0 { "+" } else { "-" }, - imag.abs(), - rewidth = max_widths.0, - imwidth = max_widths.1, - ) - } else { - format!("{:>width$.2}", real, width = max_widths.0) - } - }) - .collect(); - - for i in 0..self.rows { - if i == 0 { - write!(f, "┌")?; - } else if i == self.rows - 1 { - write!(f, "└")?; - } else { - write!(f, "│")?; - } - - for j in 0..self.cols { - write!(f, "{}", aligned[i + j * self.rows])?; - if j != self.cols - 1 { - write!(f, ", ")?; - } - } - - if i == 0 { - write!(f, "┐")?; - } else if i == self.rows - 1 { - write!(f, "┘")?; - } else { - write!(f, "│")?; - } - - if i != self.rows - 1 { - write!(f, "\n")?; - } - } - - Ok(()) - } -} diff --git a/libpsi-core/src/maths/mod.rs b/libpsi-core/src/maths/mod.rs deleted file mode 100644 index 85f4872..0000000 --- a/libpsi-core/src/maths/mod.rs +++ /dev/null @@ -1,14 +0,0 @@ -pub mod complex; -pub mod format; -pub mod matrix; -pub mod numeric; -pub mod simd; -pub mod vector; -pub mod vector_ops; - -pub use complex::*; -pub use format::*; -pub use matrix::*; -pub use numeric::*; -pub use simd::*; -pub use vector::*; diff --git a/libpsi-core/src/maths/numeric.rs b/libpsi-core/src/maths/numeric.rs deleted file mode 100644 index f5a649e..0000000 --- a/libpsi-core/src/maths/numeric.rs +++ /dev/null @@ -1,101 +0,0 @@ -use crate::Complex; -use core::ops; - -macro_rules! impl_numeric { - ($($t:ty),*) => { - $( - impl Numeric for $t { - fn zero() -> Self { - 0 as $t - } - - fn one() -> Self { - 1 as $t - } - } - )* - }; -} - -macro_rules! impl_cnumeric { - ($($t:ty),*) => { - $(impl Numeric for Complex<$t> { - fn zero() -> Self { Complex::new(0.0, 0.0) } - fn one() -> Self { Complex::new(1.0, 0.0) } - })* - }; -} - -macro_rules! impl_float { - ($($t:ty, $sqrt_fn:path, $atan2_fn:path),*) => { - $( - impl Float for $t { - fn sqrt(self) -> Self { - $sqrt_fn(self) - } - - fn atan2(y: Self, x: Self) -> Self { - $atan2_fn(y, x) - } - } - )* - }; -} - -macro_rules! impl_cfloat { - ($($t:ty, $sqrt_fn:path, $atan2_fn:path, $cos_fn:path, $sin_fn:path),*) => { - $( - impl Float for Complex<$t> { - fn sqrt(self) -> Self { - let r = self.abs(); - let theta = self.phase(); - - let sqrt_r = $sqrt_fn(r); - let sqrt_theta = theta / 2.0; - - Complex::new( - sqrt_r * $cos_fn(sqrt_theta), - sqrt_r * $sin_fn(sqrt_theta), - ) - } - - fn atan2(y: Self, x: Self) -> Self { - Complex::new( - $atan2_fn(y.real, x.real), - $atan2_fn(y.imaginary, x.imaginary), - ) - } - } - )* - }; -} - -pub trait Numeric: - Copy - + PartialOrd - + ops::Add<Output = Self> - + ops::Mul<Output = Self> - + ops::Sub<Output = Self> - + ops::Div<Output = Self> - + ops::Neg<Output = Self> - + ops::AddAssign - + ops::SubAssign - + ops::MulAssign - + ops::DivAssign -{ - fn zero() -> Self; - fn one() -> Self; -} - -impl_numeric!(i32, i64, f32, f64); -impl_cnumeric!(f32, f64); -impl_float!(f32, libm::sqrtf, libm::atan2f); -impl_float!(f64, libm::sqrt, libm::atan2); -impl_cfloat!(f32, libm::sqrtf, libm::atan2f, libm::cosf, libm::sinf); -impl_cfloat!(f64, libm::sqrt, libm::atan2, libm::cos, libm::sin); - -pub trait Integer: Numeric {} -pub trait Float: Numeric { - fn sqrt(self) -> Self; - fn atan2(y: Self, x: Self) -> Self; -} diff --git a/libpsi-core/src/maths/simd.rs b/libpsi-core/src/maths/simd.rs deleted file mode 100644 index de0370c..0000000 --- a/libpsi-core/src/maths/simd.rs +++ /dev/null @@ -1,510 +0,0 @@ -use crate::{complex, Complex}; - -#[cfg(target_arch = "x86_64")] -use std::arch::x86_64::*; - -#[cfg(target_arch = "aarch64")] -use std::arch::aarch64::*; - -#[derive(Debug, Clone, Copy, PartialEq, Eq)] -pub enum SimdCapability { - None, - #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] - Avx2, - #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] - Avx512, - #[cfg(target_arch = "aarch64")] - Neon, -} - -impl SimdCapability { - pub fn detect() -> Self { - #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] - { - if is_x86_feature_detected!("avx512f") && is_x86_feature_detected!("avx512dq") { - return SimdCapability::Avx512; - } - if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") { - return SimdCapability::Avx2; - } - } - - #[cfg(target_arch = "aarch64")] - { - return SimdCapability::Neon; - } - - #[allow(unreachable_code)] - SimdCapability::None - } - - pub fn name(&self) -> &'static str { - match self { - SimdCapability::None => "Scalar", - #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] - SimdCapability::Avx2 => "AVX2+FMA", - #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] - SimdCapability::Avx512 => "AVX-512", - #[cfg(target_arch = "aarch64")] - SimdCapability::Neon => "NEON", - } - } -} - -pub fn apply_single_qubit_gate_simd( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - let capability = SimdCapability::detect(); - - match capability { - #[cfg(target_arch = "x86_64")] - SimdCapability::Avx2 => unsafe { - apply_single_qubit_avx2(state, gate, target, num_qubits); - }, - #[cfg(target_arch = "x86_64")] - SimdCapability::Avx512 => unsafe { - apply_single_qubit_avx512(state, gate, target, num_qubits); - }, - #[cfg(target_arch = "aarch64")] - SimdCapability::Neon => unsafe { - apply_single_qubit_neon(state, gate, target, num_qubits); - }, - _ => { - apply_single_qubit_scalar(state, gate, target, num_qubits); - } - } -} - -#[cfg(target_arch = "x86_64")] -#[target_feature(enable = "avx2", enable = "fma")] -unsafe fn apply_single_qubit_avx2( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - let target_bit = num_qubits - 1 - target; - let step = 1 << target_bit; - let dim = 1 << num_qubits; - - let g00 = gate[0][0]; - let g01 = gate[0][1]; - let g10 = gate[1][0]; - let g11 = gate[1][1]; - - let pairs: Vec<(usize, usize)> = (0..dim) - .filter(|&i| (i >> target_bit) & 1 == 0) - .map(|i| (i, i | step)) - .collect(); - - let chunks = pairs.len() / 2; - - for chunk_idx in 0..chunks { - let (i0, j0) = pairs[chunk_idx * 2]; - let (i1, j1) = pairs[chunk_idx * 2 + 1]; - - let s0_re = _mm256_set_pd( - state[j1].real, - state[i1].real, - state[j0].real, - state[i0].real, - ); - let s0_im = _mm256_set_pd( - state[j1].imaginary, - state[i1].imaginary, - state[j0].imaginary, - state[i0].imaginary, - ); - - let g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real); - let g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); - let g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real); - let g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); - - let prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0)); - let prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0)); - - let prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1)); - let prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1)); - - let mut res0_re = [0.0f64; 4]; - let mut res0_im = [0.0f64; 4]; - let mut res1_re = [0.0f64; 4]; - let mut res1_im = [0.0f64; 4]; - - _mm256_storeu_pd(res0_re.as_mut_ptr(), prod0_re); - _mm256_storeu_pd(res0_im.as_mut_ptr(), prod0_im); - _mm256_storeu_pd(res1_re.as_mut_ptr(), prod1_re); - _mm256_storeu_pd(res1_im.as_mut_ptr(), prod1_im); - - state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); - state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); - state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); - state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); - } - - for &(i, j) in pairs.iter().skip(chunks * 2) { - let s0 = state[i]; - let s1 = state[j]; - - let new0 = complex!( - s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real - - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary - + s0.imaginary * g00.real - + s1.real * g01.imaginary - + s1.imaginary * g01.real - ); - - let new1 = complex!( - s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real - - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary - + s0.imaginary * g10.real - + s1.real * g11.imaginary - + s1.imaginary * g11.real - ); - - state[i] = new0; - state[j] = new1; - } -} - -#[cfg(target_arch = "x86_64")] -#[target_feature(enable = "avx512f", enable = "avx512dq")] -unsafe fn apply_single_qubit_avx512( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - let target_bit = num_qubits - 1 - target; - let step = 1 << target_bit; - let dim = 1 << num_qubits; - - let g00 = gate[0][0]; - let g01 = gate[0][1]; - let g10 = gate[1][0]; - let g11 = gate[1][1]; - - let pairs: Vec<(usize, usize)> = (0..dim) - .filter(|&i| (i >> target_bit) & 1 == 0) - .map(|i| (i, i | step)) - .collect(); - - let chunks = pairs.len() / 4; - - for chunk_idx in 0..chunks { - let base = chunk_idx * 4; - let (i0, j0) = pairs[base]; - let (i1, j1) = pairs[base + 1]; - let (i2, j2) = pairs[base + 2]; - let (i3, j3) = pairs[base + 3]; - - let s0_re = _mm512_set_pd( - state[j3].real, - state[i3].real, - state[j2].real, - state[i2].real, - state[j1].real, - state[i1].real, - state[j0].real, - state[i0].real, - ); - let s0_im = _mm512_set_pd( - state[j3].imaginary, - state[i3].imaginary, - state[j2].imaginary, - state[i2].imaginary, - state[j1].imaginary, - state[i1].imaginary, - state[j0].imaginary, - state[i0].imaginary, - ); - - let g_re_0 = _mm512_set_pd( - g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, - ); - let g_im_0 = _mm512_set_pd( - g01.imaginary, - g00.imaginary, - g01.imaginary, - g00.imaginary, - g01.imaginary, - g00.imaginary, - g01.imaginary, - g00.imaginary, - ); - let g_re_1 = _mm512_set_pd( - g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, - ); - let g_im_1 = _mm512_set_pd( - g11.imaginary, - g10.imaginary, - g11.imaginary, - g10.imaginary, - g11.imaginary, - g10.imaginary, - g11.imaginary, - g10.imaginary, - ); - - let prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0)); - let prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0)); - let prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1)); - let prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1)); - - let mut res0_re = [0.0f64; 8]; - let mut res0_im = [0.0f64; 8]; - let mut res1_re = [0.0f64; 8]; - let mut res1_im = [0.0f64; 8]; - - _mm512_storeu_pd(res0_re.as_mut_ptr(), prod0_re); - _mm512_storeu_pd(res0_im.as_mut_ptr(), prod0_im); - _mm512_storeu_pd(res1_re.as_mut_ptr(), prod1_re); - _mm512_storeu_pd(res1_im.as_mut_ptr(), prod1_im); - - state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); - state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); - state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); - state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); - state[i2] = complex!(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]); - state[j2] = complex!(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]); - state[i3] = complex!(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]); - state[j3] = complex!(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]); - } - - for &(i, j) in pairs.iter().skip(chunks * 4) { - let s0 = state[i]; - let s1 = state[j]; - - let new0 = complex!( - s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real - - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary - + s0.imaginary * g00.real - + s1.real * g01.imaginary - + s1.imaginary * g01.real - ); - - let new1 = complex!( - s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real - - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary - + s0.imaginary * g10.real - + s1.real * g11.imaginary - + s1.imaginary * g11.real - ); - - state[i] = new0; - state[j] = new1; - } -} - -#[cfg(target_arch = "aarch64")] -unsafe fn apply_single_qubit_neon( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - let target_bit = num_qubits - 1 - target; - let step = 1 << target_bit; - let dim = 1 << num_qubits; - - let g00 = gate[0][0]; - let g01 = gate[0][1]; - let g10 = gate[1][0]; - let g11 = gate[1][1]; - - let pairs: Vec<(usize, usize)> = (0..dim) - .filter(|&i| (i >> target_bit) & 1 == 0) - .map(|i| (i, i | step)) - .collect(); - - let chunks = pairs.len() / 2; - - for chunk_idx in 0..chunks { - let (i0, j0) = pairs[chunk_idx * 2]; - let (i1, j1) = pairs[chunk_idx * 2 + 1]; - - let s0_0 = state[i0]; - let s1_0 = state[j0]; - let s0_1 = state[i1]; - let s1_1 = state[j1]; - - let s0_re = vld1q_f64([s0_0.real, s0_1.real].as_ptr()); - let s0_im = vld1q_f64([s0_0.imaginary, s0_1.imaginary].as_ptr()); - let s1_re = vld1q_f64([s1_0.real, s1_1.real].as_ptr()); - let s1_im = vld1q_f64([s1_0.imaginary, s1_1.imaginary].as_ptr()); - - let g00_re = vdupq_n_f64(g00.real); - let g00_im = vdupq_n_f64(g00.imaginary); - let g01_re = vdupq_n_f64(g01.real); - let g01_im = vdupq_n_f64(g01.imaginary); - let g10_re = vdupq_n_f64(g10.real); - let g10_im = vdupq_n_f64(g10.imaginary); - let g11_re = vdupq_n_f64(g11.real); - let g11_im = vdupq_n_f64(g11.imaginary); - - let new0_re = vaddq_f64( - vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im), - vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im), - ); - let new0_im = vaddq_f64( - vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re), - vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re), - ); - - let new1_re = vaddq_f64( - vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im), - vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im), - ); - let new1_im = vaddq_f64( - vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re), - vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re), - ); - - state[i0] = complex!(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0)); - state[j0] = complex!(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0)); - state[i1] = complex!(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1)); - state[j1] = complex!(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1)); - } - - for &(i, j) in pairs.iter().skip(chunks * 2) { - let s0 = state[i]; - let s1 = state[j]; - - let new0 = complex!( - s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real - - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary - + s0.imaginary * g00.real - + s1.real * g01.imaginary - + s1.imaginary * g01.real - ); - - let new1 = complex!( - s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real - - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary - + s0.imaginary * g10.real - + s1.real * g11.imaginary - + s1.imaginary * g11.real - ); - - state[i] = new0; - state[j] = new1; - } -} - -fn apply_single_qubit_scalar( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - let target_bit = num_qubits - 1 - target; - let step = 1 << target_bit; - let dim = 1 << num_qubits; - - let g00 = gate[0][0]; - let g01 = gate[0][1]; - let g10 = gate[1][0]; - let g11 = gate[1][1]; - - for i in 0..dim { - if (i >> target_bit) & 1 == 1 { - continue; - } - - let j = i | step; - let s0 = state[i]; - let s1 = state[j]; - - let new0 = complex!( - s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real - - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary - + s0.imaginary * g00.real - + s1.real * g01.imaginary - + s1.imaginary * g01.real - ); - - let new1 = complex!( - s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real - - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary - + s0.imaginary * g10.real - + s1.real * g11.imaginary - + s1.imaginary * g11.real - ); - - state[i] = new0; - state[j] = new1; - } -} - -pub fn apply_single_qubit_gate_simd_parallel( - state: &mut [Complex<f64>], - gate: &[[Complex<f64>; 2]; 2], - target: usize, - num_qubits: usize, -) { - use rayon::prelude::*; - - let target_bit = num_qubits - 1 - target; - let step = 1 << target_bit; - let dim = 1 << num_qubits; - - let g00 = gate[0][0]; - let g01 = gate[0][1]; - let g10 = gate[1][0]; - let g11 = gate[1][1]; - - let pairs: Vec<(usize, usize)> = (0..dim) - .filter(|&i| (i >> target_bit) & 1 == 0) - .map(|i| (i, i | step)) - .collect(); - - let results: Vec<(usize, usize, Complex<f64>, Complex<f64>)> = pairs - .par_iter() - .map(|&(i, j)| { - let s0 = state[i]; - let s1 = state[j]; - - let new0 = complex!( - s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real - - s1.imaginary * g01.imaginary, - s0.real * g00.imaginary - + s0.imaginary * g00.real - + s1.real * g01.imaginary - + s1.imaginary * g01.real - ); - - let new1 = complex!( - s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real - - s1.imaginary * g11.imaginary, - s0.real * g10.imaginary - + s0.imaginary * g10.real - + s1.real * g11.imaginary - + s1.imaginary * g11.real - ); - - (i, j, new0, new1) - }) - .collect(); - - for (i, j, new0, new1) in results { - state[i] = new0; - state[j] = new1; - } -} - -pub fn get_simd_info() -> String { - let cap = SimdCapability::detect(); - format!("SIMD: {}", cap.name()) -} diff --git a/libpsi-core/src/maths/vector.rs b/libpsi-core/src/maths/vector.rs deleted file mode 100644 index 87b5353..0000000 --- a/libpsi-core/src/maths/vector.rs +++ /dev/null @@ -1,258 +0,0 @@ -use super::{Float, Matrix}; -use core::{fmt, ops}; - -#[macro_export] -macro_rules! row_vector { - ($($x:expr),*) => { - RowVector::new(vec![$($x),*]) - }; - ($($x:expr,)*) => { - RowVector::new(vec![$($x),*]) - }; -} - -#[macro_export] -macro_rules! column_vector { - ($($x:expr),*) => { - ColumnVector::new(vec![$($x),*]) - }; - ($($x:expr,)*) => { - ColumnVector::new(vec![$($x),*]) - }; -} - -pub trait Vector<T: Float> { - fn new(data: Vec<T>) -> Self; - fn get(&self, index: usize) -> T; - fn set(&mut self, index: usize, value: T); - fn size(&self) -> usize; - - fn dot(&self, other: &Self) -> T; - fn norm(&self) -> T; - - fn max(&self) -> T; - fn min(&self) -> T; - fn sum(&self) -> T; - - fn from_matrix(matrix: &Matrix<T>) -> Self; -} - -pub trait VectorMatrix<T: Float> { - fn to_matrix(&self) -> Matrix<T>; -} - -#[derive(Clone)] -pub struct VectorImpl<T: Float, const ROWS: usize, const COLS: usize>(Vec<T>); -pub type RowVector<T> = VectorImpl<T, 1, 0>; -pub type ColumnVector<T> = VectorImpl<T, 0, 1>; - -impl<T: Float> ColumnVector<T> { - pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<ColumnVector<T>> { - if matrix.cols != self.size() { - return None; - } - - let mut result = ColumnVector::new(vec![T::zero(); matrix.rows]); - - for i in 0..matrix.rows { - let mut sum = T::zero(); - for j in 0..matrix.cols { - sum = sum + (matrix.get(i, j) * self.get(j)); - } - result.set(i, sum); - } - - Some(result) - } - - pub fn transpose(&self) -> RowVector<T> { - RowVector::new(self.0.clone()) - } -} - -impl<T: Float> RowVector<T> { - pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<RowVector<T>> { - if self.size() != matrix.rows { - return None; - } - - let mut result = RowVector::new(vec![T::zero(); matrix.cols]); - - for j in 0..matrix.cols { - let mut sum = T::zero(); - for i in 0..matrix.rows { - sum = sum + (self.get(i) * matrix.get(i, j)); - } - result.set(j, sum); - } - - Some(result) - } - - pub fn transpose(&self) -> ColumnVector<T> { - ColumnVector::new(self.0.clone()) - } -} - -impl<T: Float> VectorMatrix<T> for RowVector<T> { - fn to_matrix(&self) -> Matrix<T> { - Matrix::new(1, self.size(), self.0.clone()) - } -} - -impl<T: Float> VectorMatrix<T> for ColumnVector<T> { - fn to_matrix(&self) -> Matrix<T> { - Matrix::new(self.size(), 1, self.0.clone()) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> Vector<T> for VectorImpl<T, ROWS, COLS> { - fn from_matrix(matrix: &Matrix<T>) -> Self { - Self::new(matrix.data.clone()) - } - - fn new(data: Vec<T>) -> Self { - Self(data) - } - - fn get(&self, index: usize) -> T { - self.0[index] - } - - fn set(&mut self, index: usize, value: T) { - self.0[index] = value; - } - - fn size(&self) -> usize { - self.0.len() - } - - fn dot(&self, other: &Self) -> T { - self.0 - .iter() - .zip(other.0.iter()) - .map(|(a, b)| *a * *b) - .fold(T::zero(), |acc, x| acc + x) - } - - fn norm(&self) -> T { - self.0 - .iter() - .map(|x| *x * *x) - .fold(T::zero(), |acc, x| acc + x) - .sqrt() - } - - fn max(&self) -> T { - *self - .0 - .iter() - .max_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap_or(&T::zero()) - } - - fn min(&self) -> T { - *self - .0 - .iter() - .min_by(|a, b| a.partial_cmp(b).unwrap()) - .unwrap_or(&T::zero()) - } - - fn sum(&self) -> T { - self.0.iter().fold(T::zero(), |acc, x| acc + *x) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> VectorImpl<T, ROWS, COLS> { - pub fn add_to(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> { - if self.size() != other.size() { - return None; - } - - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let sum = self.get(i) + other.get(i); - result.set(i, sum); - } - - Some(result) - } - - pub fn subtract(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> { - if self.size() != other.size() { - return None; - } - - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let sum = self.get(i) - other.get(i); - result.set(i, sum); - } - - Some(result) - } - - pub fn scale(&self, scalar: T) -> VectorImpl<T, ROWS, COLS> { - let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]); - - for i in 0..self.size() { - let product = self.get(i) * scalar; - result.set(i, product); - } - - result - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Index<usize> - for VectorImpl<T, ROWS, COLS> -{ - type Output = T; - - fn index(&self, index: usize) -> &Self::Output { - &self.0[index] - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::IndexMut<usize> - for VectorImpl<T, ROWS, COLS> -{ - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.0[index] - } -} - -impl<T: Float + fmt::Debug> fmt::Debug for RowVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "RowVector({:?})", self.0) - } -} - -impl<T: Float + fmt::Debug> fmt::Debug for ColumnVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "ColumnVector({:?})", self.0) - } -} - -impl<T: Float + fmt::Display> fmt::Display for RowVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "[{}]", - self.0 - .iter() - .map(|x| x.to_string()) - .collect::<Vec<String>>() - .join(", ") - ) - } -} - -impl<T: Float + fmt::Display> fmt::Display for ColumnVector<T> { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "{}", self.to_matrix()) - } -} diff --git a/libpsi-core/src/maths/vector_ops.rs b/libpsi-core/src/maths/vector_ops.rs deleted file mode 100644 index 1e7148c..0000000 --- a/libpsi-core/src/maths/vector_ops.rs +++ /dev/null @@ -1,107 +0,0 @@ -use super::{Float, Matrix}; -use crate::{ColumnVector, RowVector, VectorImpl}; -use core::ops; - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Add<&VectorImpl<T, ROWS, COLS>> - for VectorImpl<T, ROWS, COLS> -{ - type Output = Option<VectorImpl<T, ROWS, COLS>>; - - fn add(self, other: &VectorImpl<T, ROWS, COLS>) -> Self::Output { - self.add_to(other) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Sub<&VectorImpl<T, ROWS, COLS>> - for VectorImpl<T, ROWS, COLS> -{ - type Output = Option<VectorImpl<T, ROWS, COLS>>; - - fn sub(self, other: &VectorImpl<T, ROWS, COLS>) -> Self::Output { - self.subtract(other) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Mul<T> for VectorImpl<T, ROWS, COLS> { - type Output = VectorImpl<T, ROWS, COLS>; - - fn mul(self, scalar: T) -> Self::Output { - self.scale(scalar) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::Div<T> for VectorImpl<T, ROWS, COLS> { - type Output = VectorImpl<T, ROWS, COLS>; - - fn div(self, scalar: T) -> Self::Output { - self.scale(T::one() / scalar) - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::AddAssign<VectorImpl<T, ROWS, COLS>> - for VectorImpl<T, ROWS, COLS> -{ - fn add_assign(&mut self, other: VectorImpl<T, ROWS, COLS>) { - if let Some(result) = self.add_to(&other) { - *self = result; - } - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::SubAssign<VectorImpl<T, ROWS, COLS>> - for VectorImpl<T, ROWS, COLS> -{ - fn sub_assign(&mut self, other: VectorImpl<T, ROWS, COLS>) { - if let Some(result) = self.subtract(&other) { - *self = result; - } - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::MulAssign<T> - for VectorImpl<T, ROWS, COLS> -{ - fn mul_assign(&mut self, scalar: T) { - *self = self.scale(scalar); - } -} - -impl<T: Float, const ROWS: usize, const COLS: usize> ops::DivAssign<T> - for VectorImpl<T, ROWS, COLS> -{ - fn div_assign(&mut self, scalar: T) { - *self = self.scale(T::one() / scalar); - } -} - -impl<T: Float> ops::Mul<&Matrix<T>> for RowVector<T> { - type Output = Option<RowVector<T>>; - - fn mul(self, matrix: &Matrix<T>) -> Self::Output { - self.mul_matrix(matrix) - } -} - -impl<T: Float> ops::Mul<&Matrix<T>> for ColumnVector<T> { - type Output = Option<ColumnVector<T>>; - - fn mul(self, matrix: &Matrix<T>) -> Self::Output { - self.mul_matrix(matrix) - } -} - -impl<T: Float> ops::MulAssign<&Matrix<T>> for RowVector<T> { - fn mul_assign(&mut self, matrix: &Matrix<T>) { - if let Some(result) = self.mul_matrix(matrix) { - *self = result; - } - } -} - -impl<T: Float> ops::MulAssign<&Matrix<T>> for ColumnVector<T> { - fn mul_assign(&mut self, matrix: &Matrix<T>) { - if let Some(result) = self.mul_matrix(matrix) { - *self = result; - } - } -} |