1use crate::prelude::SimulatorError;
10use scirs2_core::ndarray::Array2;
11use scirs2_core::parallel_ops::*;
12use scirs2_core::Complex64;
13use std::collections::HashMap;
14use std::fmt;
15
16use crate::error::Result;
17use crate::trotter::{DynamicCircuit, Hamiltonian, TrotterDecomposer, TrotterMethod};
18use quantrs2_core::gate::{multi::*, single::*, GateOp};
19use quantrs2_core::qubit::QubitId;
20
21#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
23pub enum PauliOperator {
24 I,
26 X,
28 Y,
30 Z,
32}
33
34impl fmt::Display for PauliOperator {
35 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
36 match self {
37 Self::I => write!(f, "I"),
38 Self::X => write!(f, "X"),
39 Self::Y => write!(f, "Y"),
40 Self::Z => write!(f, "Z"),
41 }
42 }
43}
44
45impl PauliOperator {
46 pub fn from_str(s: &str) -> Result<Self> {
48 match s.to_uppercase().as_str() {
49 "I" => Ok(Self::I),
50 "X" => Ok(Self::X),
51 "Y" => Ok(Self::Y),
52 "Z" => Ok(Self::Z),
53 _ => Err(SimulatorError::InvalidInput(format!(
54 "Invalid Pauli operator: {s}"
55 ))),
56 }
57 }
58
59 pub fn matrix(&self) -> Array2<Complex64> {
61 match self {
62 Self::I => Array2::from_shape_vec(
63 (2, 2),
64 vec![
65 Complex64::new(1.0, 0.0),
66 Complex64::new(0.0, 0.0),
67 Complex64::new(0.0, 0.0),
68 Complex64::new(1.0, 0.0),
69 ],
70 )
71 .unwrap(),
72 Self::X => Array2::from_shape_vec(
73 (2, 2),
74 vec![
75 Complex64::new(0.0, 0.0),
76 Complex64::new(1.0, 0.0),
77 Complex64::new(1.0, 0.0),
78 Complex64::new(0.0, 0.0),
79 ],
80 )
81 .unwrap(),
82 Self::Y => Array2::from_shape_vec(
83 (2, 2),
84 vec![
85 Complex64::new(0.0, 0.0),
86 Complex64::new(0.0, -1.0),
87 Complex64::new(0.0, 1.0),
88 Complex64::new(0.0, 0.0),
89 ],
90 )
91 .unwrap(),
92 Self::Z => Array2::from_shape_vec(
93 (2, 2),
94 vec![
95 Complex64::new(1.0, 0.0),
96 Complex64::new(0.0, 0.0),
97 Complex64::new(0.0, 0.0),
98 Complex64::new(-1.0, 0.0),
99 ],
100 )
101 .unwrap(),
102 }
103 }
104
105 pub fn commutes_with(&self, other: &Self) -> bool {
107 match (self, other) {
108 (Self::I, _) | (_, Self::I) => true,
109 (a, b) if a == b => true,
110 _ => false,
111 }
112 }
113
114 pub const fn multiply(&self, other: &Self) -> (Self, Complex64) {
116 match (self, other) {
117 (Self::I, p) | (p, Self::I) => (*p, Complex64::new(1.0, 0.0)),
118 (Self::X, Self::X) | (Self::Y, Self::Y) | (Self::Z, Self::Z) => {
119 (Self::I, Complex64::new(1.0, 0.0))
120 }
121 (Self::X, Self::Y) => (Self::Z, Complex64::new(0.0, 1.0)),
122 (Self::Y, Self::X) => (Self::Z, Complex64::new(0.0, -1.0)),
123 (Self::Y, Self::Z) => (Self::X, Complex64::new(0.0, 1.0)),
124 (Self::Z, Self::Y) => (Self::X, Complex64::new(0.0, -1.0)),
125 (Self::Z, Self::X) => (Self::Y, Complex64::new(0.0, 1.0)),
126 (Self::X, Self::Z) => (Self::Y, Complex64::new(0.0, -1.0)),
127 }
128 }
129}
130
131#[derive(Debug, Clone)]
133pub struct PauliString {
134 pub operators: Vec<PauliOperator>,
136 pub coefficient: Complex64,
138 pub num_qubits: usize,
140}
141
142impl PauliString {
143 pub fn new(num_qubits: usize) -> Self {
145 Self {
146 operators: vec![PauliOperator::I; num_qubits],
147 coefficient: Complex64::new(1.0, 0.0),
148 num_qubits,
149 }
150 }
151
152 pub fn from_string(pauli_str: &str, coefficient: Complex64) -> Result<Self> {
154 let operators: Result<Vec<_>> = pauli_str
155 .chars()
156 .map(|c| PauliOperator::from_str(&c.to_string()))
157 .collect();
158
159 Ok(Self {
160 operators: operators?,
161 coefficient,
162 num_qubits: pauli_str.len(),
163 })
164 }
165
166 pub fn from_ops(
168 num_qubits: usize,
169 ops: &[(usize, PauliOperator)],
170 coefficient: Complex64,
171 ) -> Result<Self> {
172 let mut pauli_string = Self::new(num_qubits);
173 pauli_string.coefficient = coefficient;
174
175 for &(qubit, op) in ops {
176 if qubit >= num_qubits {
177 return Err(SimulatorError::IndexOutOfBounds(qubit));
178 }
179 pauli_string.operators[qubit] = op;
180 }
181
182 Ok(pauli_string)
183 }
184
185 pub fn set_operator(&mut self, qubit: usize, op: PauliOperator) -> Result<()> {
187 if qubit >= self.num_qubits {
188 return Err(SimulatorError::IndexOutOfBounds(qubit));
189 }
190 self.operators[qubit] = op;
191 Ok(())
192 }
193
194 pub fn get_operator(&self, qubit: usize) -> Result<PauliOperator> {
196 if qubit >= self.num_qubits {
197 return Err(SimulatorError::IndexOutOfBounds(qubit));
198 }
199 Ok(self.operators[qubit])
200 }
201
202 pub fn non_identity_ops(&self) -> Vec<(usize, PauliOperator)> {
204 self.operators
205 .iter()
206 .enumerate()
207 .filter(|(_, &op)| op != PauliOperator::I)
208 .map(|(i, &op)| (i, op))
209 .collect()
210 }
211
212 pub fn commutes_with(&self, other: &Self) -> bool {
214 if self.num_qubits != other.num_qubits {
215 return false;
216 }
217
218 let mut anti_commute_count = 0;
219 for i in 0..self.num_qubits {
220 if !self.operators[i].commutes_with(&other.operators[i]) {
221 anti_commute_count += 1;
222 }
223 }
224
225 anti_commute_count % 2 == 0
227 }
228
229 pub fn multiply(&self, other: &Self) -> Result<Self> {
231 if self.num_qubits != other.num_qubits {
232 return Err(SimulatorError::DimensionMismatch(format!(
233 "Pauli strings have different lengths: {} vs {}",
234 self.num_qubits, other.num_qubits
235 )));
236 }
237
238 let mut result = Self::new(self.num_qubits);
239 let mut total_phase = self.coefficient * other.coefficient;
240
241 for i in 0..self.num_qubits {
242 let (op, phase) = self.operators[i].multiply(&other.operators[i]);
243 result.operators[i] = op;
244 total_phase *= phase;
245 }
246
247 result.coefficient = total_phase;
248 Ok(result)
249 }
250
251 pub fn weight(&self) -> usize {
253 self.operators
254 .iter()
255 .filter(|&&op| op != PauliOperator::I)
256 .count()
257 }
258
259 pub fn pauli_string(&self) -> String {
261 self.operators.iter().map(|op| op.to_string()).collect()
262 }
263
264 pub fn evolution_circuit(&self, time: f64) -> Result<DynamicCircuit> {
266 let mut circuit = DynamicCircuit::new(self.num_qubits);
267
268 if self.weight() == 0 {
269 return Ok(circuit);
271 }
272
273 let non_identity = self.non_identity_ops();
274 let angle = -2.0 * self.coefficient.re * time;
275
276 if non_identity.len() == 1 {
277 let (qubit, op) = non_identity[0];
279 match op {
280 PauliOperator::X => circuit.add_gate(Box::new(RotationX {
281 target: QubitId::new(qubit as u32),
282 theta: angle,
283 }))?,
284 PauliOperator::Y => circuit.add_gate(Box::new(RotationY {
285 target: QubitId::new(qubit as u32),
286 theta: angle,
287 }))?,
288 PauliOperator::Z => circuit.add_gate(Box::new(RotationZ {
289 target: QubitId::new(qubit as u32),
290 theta: angle,
291 }))?,
292 PauliOperator::I => {} }
294 } else {
295 for &(qubit, op) in &non_identity {
299 match op {
300 PauliOperator::X => circuit.add_gate(Box::new(Hadamard {
301 target: QubitId::new(qubit as u32),
302 }))?,
303 PauliOperator::Y => {
304 circuit.add_gate(Box::new(Hadamard {
305 target: QubitId::new(qubit as u32),
306 }))?;
307 circuit.add_gate(Box::new(Phase {
308 target: QubitId::new(qubit as u32),
309 }))?;
310 }
311 PauliOperator::Z => {} PauliOperator::I => {} }
314 }
315
316 for i in 0..non_identity.len() - 1 {
318 circuit.add_gate(Box::new(CNOT {
319 control: QubitId::new(non_identity[i].0 as u32),
320 target: QubitId::new(non_identity[i + 1].0 as u32),
321 }))?;
322 }
323
324 circuit.add_gate(Box::new(RotationZ {
326 target: QubitId::new(non_identity[non_identity.len() - 1].0 as u32),
327 theta: angle,
328 }))?;
329
330 for i in (0..non_identity.len() - 1).rev() {
332 circuit.add_gate(Box::new(CNOT {
333 control: QubitId::new(non_identity[i].0 as u32),
334 target: QubitId::new(non_identity[i + 1].0 as u32),
335 }))?;
336 }
337
338 for &(qubit, op) in non_identity.iter().rev() {
340 match op {
341 PauliOperator::X => circuit.add_gate(Box::new(Hadamard {
342 target: QubitId::new(qubit as u32),
343 }))?,
344 PauliOperator::Y => {
345 circuit.add_gate(Box::new(PhaseDagger {
346 target: QubitId::new(qubit as u32),
347 }))?;
348 circuit.add_gate(Box::new(Hadamard {
349 target: QubitId::new(qubit as u32),
350 }))?;
351 }
352 PauliOperator::Z => {} PauliOperator::I => {} }
355 }
356 }
357
358 Ok(circuit)
359 }
360}
361
362impl fmt::Display for PauliString {
363 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
364 write!(f, "({}) {}", self.coefficient, self.pauli_string())
365 }
366}
367
368#[derive(Debug, Clone)]
370pub struct PauliOperatorSum {
371 pub terms: Vec<PauliString>,
373 pub num_qubits: usize,
375}
376
377impl PauliOperatorSum {
378 pub const fn new(num_qubits: usize) -> Self {
380 Self {
381 terms: Vec::new(),
382 num_qubits,
383 }
384 }
385
386 pub fn add_term(&mut self, pauli_string: PauliString) -> Result<()> {
388 if pauli_string.num_qubits != self.num_qubits {
389 return Err(SimulatorError::DimensionMismatch(format!(
390 "Pauli string has {} qubits, expected {}",
391 pauli_string.num_qubits, self.num_qubits
392 )));
393 }
394 self.terms.push(pauli_string);
395 Ok(())
396 }
397
398 pub fn simplify(&mut self) {
400 let mut simplified_terms: HashMap<String, Complex64> = HashMap::new();
401
402 for term in &self.terms {
403 let key = format!("{term}");
404 *simplified_terms
405 .entry(key)
406 .or_insert(Complex64::new(0.0, 0.0)) += term.coefficient;
407 }
408
409 self.terms.clear();
410 for (pauli_str, coeff) in simplified_terms {
411 if coeff.norm() > 1e-15 {
412 if let Ok(term) = PauliString::from_string(&pauli_str, coeff) {
413 self.terms.push(term);
414 }
415 }
416 }
417 }
418
419 pub fn to_hamiltonian(&self) -> Hamiltonian {
421 let mut ham = Hamiltonian::new(self.num_qubits);
422
423 for term in &self.terms {
424 let non_identity = term.non_identity_ops();
425
426 match non_identity.len() {
427 0 => {} 1 => {
429 let (qubit, op) = non_identity[0];
430 let pauli_str = match op {
431 PauliOperator::X => "X",
432 PauliOperator::Y => "Y",
433 PauliOperator::Z => "Z",
434 PauliOperator::I => continue,
435 };
436 let _ = ham.add_single_pauli(qubit, pauli_str, term.coefficient.re);
437 }
438 2 => {
439 let (q1, op1) = non_identity[0];
440 let (q2, op2) = non_identity[1];
441 let pauli1 = match op1 {
442 PauliOperator::X => "X",
443 PauliOperator::Y => "Y",
444 PauliOperator::Z => "Z",
445 PauliOperator::I => continue,
446 };
447 let pauli2 = match op2 {
448 PauliOperator::X => "X",
449 PauliOperator::Y => "Y",
450 PauliOperator::Z => "Z",
451 PauliOperator::I => continue,
452 };
453 let _ = ham.add_two_pauli(q1, q2, pauli1, pauli2, term.coefficient.re);
454 }
455 _ => {
456 let qubits: Vec<usize> = non_identity.iter().map(|&(q, _)| q).collect();
458 let paulis: Vec<String> = non_identity
459 .iter()
460 .map(|&(_, op)| match op {
461 PauliOperator::X => "X".to_string(),
462 PauliOperator::Y => "Y".to_string(),
463 PauliOperator::Z => "Z".to_string(),
464 PauliOperator::I => "I".to_string(),
465 })
466 .collect();
467 let _ = ham.add_pauli_string(qubits, paulis, term.coefficient.re);
468 }
469 }
470 }
471
472 ham
473 }
474
475 pub fn time_evolution_circuit(
477 &self,
478 time: f64,
479 trotter_steps: usize,
480 method: TrotterMethod,
481 ) -> Result<DynamicCircuit> {
482 let hamiltonian = self.to_hamiltonian();
483 let decomposer = TrotterDecomposer::new(method, trotter_steps);
484 decomposer.decompose(&hamiltonian, time)
485 }
486
487 pub fn exact_evolution_circuit(&self, time: f64) -> Result<DynamicCircuit> {
489 let mut circuit = DynamicCircuit::new(self.num_qubits);
490
491 for term in &self.terms {
492 let term_circuit = term.evolution_circuit(time)?;
493 for gate in term_circuit.gates() {
494 circuit.add_gate(gate.clone())?;
495 }
496 }
497
498 Ok(circuit)
499 }
500}
501
502impl fmt::Display for PauliOperatorSum {
503 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
504 if self.terms.is_empty() {
505 write!(f, "0")
506 } else {
507 let term_strs: Vec<String> = self.terms.iter().map(|t| format!("{t}")).collect();
508 write!(f, "{}", term_strs.join(" + "))
509 }
510 }
511}
512
513pub struct PauliUtils;
515
516impl PauliUtils {
517 pub fn single_qubit(
519 num_qubits: usize,
520 qubit: usize,
521 op: PauliOperator,
522 coeff: Complex64,
523 ) -> Result<PauliString> {
524 PauliString::from_ops(num_qubits, &[(qubit, op)], coeff)
525 }
526
527 pub fn all_x(num_qubits: usize, coeff: Complex64) -> PauliString {
529 let mut pauli = PauliString::new(num_qubits);
530 pauli.coefficient = coeff;
531 for i in 0..num_qubits {
532 pauli.operators[i] = PauliOperator::X;
533 }
534 pauli
535 }
536
537 pub fn all_z(num_qubits: usize, coeff: Complex64) -> PauliString {
539 let mut pauli = PauliString::new(num_qubits);
540 pauli.coefficient = coeff;
541 for i in 0..num_qubits {
542 pauli.operators[i] = PauliOperator::Z;
543 }
544 pauli
545 }
546
547 pub fn random(num_qubits: usize, weight: usize, coeff: Complex64) -> Result<PauliString> {
549 if weight > num_qubits {
550 return Err(SimulatorError::InvalidInput(
551 "Weight cannot exceed number of qubits".to_string(),
552 ));
553 }
554
555 let mut pauli = PauliString::new(num_qubits);
556 pauli.coefficient = coeff;
557
558 let mut positions: Vec<usize> = (0..num_qubits).collect();
560 fastrand::shuffle(&mut positions);
561
562 let ops = [PauliOperator::X, PauliOperator::Y, PauliOperator::Z];
563
564 for i in 0..weight {
565 pauli.operators[positions[i]] = ops[fastrand::usize(0..3)];
566 }
567
568 Ok(pauli)
569 }
570
571 pub fn are_mutually_commuting(pauli_strings: &[PauliString]) -> bool {
573 for i in 0..pauli_strings.len() {
574 for j in i + 1..pauli_strings.len() {
575 if !pauli_strings[i].commutes_with(&pauli_strings[j]) {
576 return false;
577 }
578 }
579 }
580 true
581 }
582
583 pub fn maximal_commuting_set(pauli_strings: &[PauliString]) -> Vec<usize> {
585 let mut commuting_set = Vec::new();
586
587 for (i, pauli) in pauli_strings.iter().enumerate() {
588 let mut can_add = true;
589 for &j in &commuting_set {
590 if !pauli.commutes_with(&pauli_strings[j]) {
591 can_add = false;
592 break;
593 }
594 }
595
596 if can_add {
597 commuting_set.push(i);
598 }
599 }
600
601 commuting_set
602 }
603}
604
605#[cfg(test)]
606mod tests {
607 use super::*;
608
609 #[test]
610 fn test_pauli_operator_multiply() {
611 let (result, phase) = PauliOperator::X.multiply(&PauliOperator::Y);
612 assert_eq!(result, PauliOperator::Z);
613 assert_eq!(phase, Complex64::new(0.0, 1.0));
614
615 let (result, phase) = PauliOperator::Y.multiply(&PauliOperator::X);
616 assert_eq!(result, PauliOperator::Z);
617 assert_eq!(phase, Complex64::new(0.0, -1.0));
618 }
619
620 #[test]
621 fn test_pauli_string_creation() {
622 let pauli = PauliString::from_string("XYZ", Complex64::new(1.0, 0.0)).unwrap();
623 assert_eq!(pauli.num_qubits, 3);
624 assert_eq!(pauli.operators[0], PauliOperator::X);
625 assert_eq!(pauli.operators[1], PauliOperator::Y);
626 assert_eq!(pauli.operators[2], PauliOperator::Z);
627 }
628
629 #[test]
630 fn test_pauli_string_multiply() {
631 let p1 = PauliString::from_string("XY", Complex64::new(1.0, 0.0)).unwrap();
632 let p2 = PauliString::from_string("YZ", Complex64::new(1.0, 0.0)).unwrap();
633
634 let result = p1.multiply(&p2).unwrap();
635 assert_eq!(result.operators[0], PauliOperator::Z);
636 assert_eq!(result.operators[1], PauliOperator::X);
637 assert_eq!(result.coefficient, Complex64::new(-1.0, 0.0));
638 }
639
640 #[test]
641 fn test_pauli_string_commutation() {
642 let p1 = PauliString::from_string("XY", Complex64::new(1.0, 0.0)).unwrap();
643 let p2 = PauliString::from_string("ZI", Complex64::new(1.0, 0.0)).unwrap();
644 let p3 = PauliString::from_string("XI", Complex64::new(1.0, 0.0)).unwrap();
645
646 assert!(!p1.commutes_with(&p2)); assert!(p1.commutes_with(&p3)); }
649
650 #[test]
651 fn test_pauli_string_weight() {
652 let p1 = PauliString::from_string("XIYZ", Complex64::new(1.0, 0.0)).unwrap();
653 assert_eq!(p1.weight(), 3);
654
655 let p2 = PauliString::from_string("IIII", Complex64::new(1.0, 0.0)).unwrap();
656 assert_eq!(p2.weight(), 0);
657 }
658
659 #[test]
660 fn test_pauli_operator_sum() {
661 let mut sum = PauliOperatorSum::new(2);
662
663 let p1 = PauliString::from_string("XX", Complex64::new(1.0, 0.0)).unwrap();
664 let p2 = PauliString::from_string("YY", Complex64::new(0.5, 0.0)).unwrap();
665
666 sum.add_term(p1).unwrap();
667 sum.add_term(p2).unwrap();
668
669 assert_eq!(sum.terms.len(), 2);
670 }
671
672 #[test]
673 fn test_evolution_circuit_single_qubit() {
674 let pauli = PauliString::from_string("X", Complex64::new(1.0, 0.0)).unwrap();
675 let circuit = pauli.evolution_circuit(1.0).unwrap();
676 assert!(circuit.gate_count() > 0);
677 }
678
679 #[test]
680 fn test_evolution_circuit_multi_qubit() {
681 let pauli = PauliString::from_string("XYZ", Complex64::new(1.0, 0.0)).unwrap();
682 let circuit = pauli.evolution_circuit(1.0).unwrap();
683 assert!(circuit.gate_count() > 0);
684 }
685
686 #[test]
687 fn test_utils_single_qubit() {
688 let pauli =
689 PauliUtils::single_qubit(3, 1, PauliOperator::X, Complex64::new(1.0, 0.0)).unwrap();
690 assert_eq!(pauli.operators[0], PauliOperator::I);
691 assert_eq!(pauli.operators[1], PauliOperator::X);
692 assert_eq!(pauli.operators[2], PauliOperator::I);
693 }
694
695 #[test]
696 fn test_mutually_commuting() {
697 let p1 = PauliString::from_string("XX", Complex64::new(1.0, 0.0)).unwrap();
698 let p2 = PauliString::from_string("ZZ", Complex64::new(1.0, 0.0)).unwrap();
699 let p3 = PauliString::from_string("YY", Complex64::new(1.0, 0.0)).unwrap();
700
701 let paulis = vec![p1, p2, p3];
702 assert!(PauliUtils::are_mutually_commuting(&paulis));
703 }
704}