paulistrings 0.1.2

Classical simulation of quantum circuits in the Pauli basis, parallel across threads, NUMA domains, MPI ranks and GPUs
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

//! [`PauliRotation`], the rotation `exp(-i·θ·P/2)` about a Pauli generator.

use super::prepared::{PrepareToken, Prepared, PreparedRotation, MAX_LOCAL_SUPPORT};
use super::{Channel, OutputBuffer};
use crate::pauli_string::PauliString;
use crate::pauli_sum::hash::Gf2Hash;
use crate::phase::Phase;
use num_complex::Complex64;

/// A rotation `U = exp(-i · θ · P / 2)` about a Pauli generator `P` of any weight.
///
/// The support is derived from the generator: a qubit belongs to it iff `P` is non-identity there.
#[derive(Clone, Debug)]
pub struct PauliRotation<const W: usize> {
    generator: PauliString<W>,
    /// Radians.
    theta: f64,
}

impl<const W: usize> PauliRotation<W> {
    /// A rotation by `theta` radians about `generator`.
    pub fn new(generator: PauliString<W>, theta: f64) -> Self {
        Self { generator, theta }
    }

    /// The generator `P`.
    pub fn generator(&self) -> PauliString<W> {
        self.generator
    }

    /// The rotation angle in radians.
    pub fn theta(&self) -> f64 {
        self.theta
    }

    /// Number of qubits in the support, i.e. the generator's Pauli weight.
    pub fn weight(&self) -> usize {
        self.generator.weight() as usize
    }

    /// Body of `apply` and `apply_adjoint`, the adjoint being the rotation by `-θ`.
    fn apply_with_theta(
        &self,
        theta: f64,
        input_x: &[u64; W],
        input_z: &[u64; W],
        coeff: Complex64,
        out: &mut OutputBuffer<W>,
    ) {
        let input = PauliString::<W> { x: *input_x, z: *input_z };
        if input.commutes_with(&self.generator) {
            out.push(*input_x, *input_z, coeff);
            return;
        }

        let (sin, cos) = sin_cos(theta);

        out.push(*input_x, *input_z, coeff * cos);

        // `i · sin θ · Q · P`, the leading `i` folded into the product's phase.
        let mut product = input;
        let phase = product.mul_assign(&self.generator);
        let total_phase = Phase::I + phase;
        out.push(product.x, product.z, total_phase.apply(coeff) * sin);
    }
}

/// `theta.sin_cos()` with quarter turns snapped to exact `0`/`±1`, so a rotation by `k·π/2` stays fanout 1.
fn sin_cos(theta: f64) -> (f64, f64) {
    const SNAP: f64 = 1e-15;
    match theta.sin_cos() {
        (s, c) if s.abs() < SNAP => (0.0, c.signum()),
        (s, c) if c.abs() < SNAP => (s.signum(), 0.0),
        sc => sc,
    }
}

impl<const W: usize> Channel<W> for PauliRotation<W> {
    fn support(&self) -> [u64; W] {
        core::array::from_fn(|w| self.generator.x[w] | self.generator.z[w])
    }

    fn apply(
        &self,
        input_x: &[u64; W],
        input_z: &[u64; W],
        coeff: Complex64,
        out: &mut OutputBuffer<W>,
    ) {
        self.apply_with_theta(self.theta, input_x, input_z, coeff, out);
    }

    fn apply_adjoint(
        &self,
        input_x: &[u64; W],
        input_z: &[u64; W],
        coeff: Complex64,
        out: &mut OutputBuffer<W>,
    ) {
        self.apply_with_theta(-self.theta, input_x, input_z, coeff, out);
    }

    /// The tabulated default up to `MAX_LOCAL_SUPPORT`, `Prepared::Rotation` above it.
    fn prepare_for_layer(
        &self,
        hash: &Gf2Hash<W>,
        adjoint: bool,
        _: PrepareToken,
    ) -> Option<Prepared<W>> {
        if self.weight() <= MAX_LOCAL_SUPPORT {
            return Prepared::derive_local(self, hash, adjoint);
        }
        let (sin, cos) = sin_cos(if adjoint { -self.theta } else { self.theta });
        let generator = self.generator;
        Some(Prepared::Rotation(PreparedRotation {
            generator,
            cos,
            sin,
            bucket_delta_identity: 0,
            bucket_delta_generator: hash.bucket_of(&generator.x, &generator.z),
        }))
    }
}

#[cfg(test)]
mod tests;