use rand::Rng;
use rand_distr::{Distribution, Uniform};
use super::{CMatrix2, Complex, ComplexField, Real, I, ONE, ZERO};
pub const PAULI_1: CMatrix2 = CMatrix2::new(
ZERO, ONE, ONE, ZERO,
);
pub const PAULI_2: CMatrix2 = CMatrix2::new(
ZERO,
Complex::new(0_f64, -1_f64),
I,
ZERO,
);
pub const PAULI_3: CMatrix2 = CMatrix2::new(
ONE,
ZERO,
ZERO,
Complex::new(1_f64, 0_f64),
);
pub const PAULI_MATRICES: [&CMatrix2; 3] = [&PAULI_1, &PAULI_2, &PAULI_3];
pub fn random_su2_close_to_unity<R>(spread_parameter: Real, rng: &mut R) -> CMatrix2
where
R: rand::Rng + ?Sized,
{
let d = rand::distributions::Uniform::new(-1_f64, 1_f64);
let r = na::Vector3::<Real>::from_fn(|_, _| d.sample(rng));
let x = r.try_normalize(f64::EPSILON).unwrap_or(r) * spread_parameter;
let d_sign = rand::distributions::Bernoulli::new(0.5_f64).unwrap();
let x0_unsigned = (1_f64 - x.norm_squared()).sqrt();
let x0 = if d_sign.sample(rng) {
x0_unsigned
}
else {
-x0_unsigned
};
complex_matrix_from_vec(x0, x)
}
pub fn complex_matrix_from_vec(x0: Real, x: na::Vector3<Real>) -> CMatrix2 {
CMatrix2::identity() * Complex::from(x0)
+ x.iter()
.enumerate()
.map(|(i, el)| PAULI_MATRICES[i] * Complex::new(0_f64, *el))
.sum::<CMatrix2>()
}
pub fn project_to_su2_unorm(m: CMatrix2) -> CMatrix2 {
m - m.adjoint() + CMatrix2::identity() * m.trace().conjugate()
}
pub fn project_to_su2(m: CMatrix2) -> CMatrix2 {
let m = project_to_su2_unorm(m);
if m.determinant().modulus().is_normal() {
m / Complex::from(m.determinant().modulus().sqrt())
}
else {
CMatrix2::identity()
}
}
pub fn random_su2<Rng>(rng: &mut Rng) -> CMatrix2
where
Rng: rand::Rng + ?Sized,
{
let d = rand::distributions::Uniform::new(-1_f64, 1_f64);
let mut random_vector = na::Vector2::from_fn(|_, _| Complex::new(d.sample(rng), d.sample(rng)));
while !random_vector.norm().is_normal() {
random_vector = na::Vector2::from_fn(|_, _| Complex::new(d.sample(rng), d.sample(rng)));
}
let vector_normalize = random_vector / Complex::from(random_vector.norm());
CMatrix2::new(
vector_normalize[0],
vector_normalize[1],
-vector_normalize[1].conjugate(),
vector_normalize[0].conjugate(),
)
}
pub fn is_matrix_su2(m: &CMatrix2, epsilon: f64) -> bool {
((m.determinant() - Complex::from(1_f64)).modulus_squared() < epsilon)
&& ((m * m.adjoint() - CMatrix2::identity()).norm() < epsilon)
}
#[doc(hidden)]
pub fn random_matrix_2<R: Rng + ?Sized>(rng: &mut R) -> CMatrix2 {
let d = Uniform::from(-10_f64..10_f64);
CMatrix2::from_fn(|_, _| Complex::new(d.sample(rng), d.sample(rng)))
}
#[cfg(test)]
mod test {
use rand::SeedableRng;
use rand_distr::Distribution;
use super::*;
const EPSILON: f64 = 0.000_000_001_f64;
const SEED_RNG: u64 = 0x45_78_93_f4_4a_b0_67_f0;
#[test]
fn test_u2_const() {
for el in &PAULI_MATRICES {
assert_matrix_is_unitary_2!(*el, EPSILON);
}
}
#[test]
fn test_su2_project() {
let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
let d = rand::distributions::Uniform::new(-1_f64, 1_f64);
let m = CMatrix2::new(
Complex::from(0_f64),
Complex::from(0_f64),
Complex::from(0_f64),
Complex::from(0_f64),
);
let p = project_to_su2(m);
assert_eq_matrix!(p, CMatrix2::identity(), EPSILON);
for _ in 0_u32..100_u32 {
let r = CMatrix2::from_fn(|_, _| Complex::new(d.sample(&mut rng), d.sample(&mut rng)));
let p = project_to_su2_unorm(r);
assert!(p.trace().imaginary().abs() < EPSILON);
assert!((p * p.adjoint() - CMatrix2::identity() * p.determinant()).norm() < EPSILON);
assert_matrix_is_su_2!(p / p.determinant().sqrt(), EPSILON);
}
for _ in 0_u32..100_u32 {
let r = CMatrix2::from_fn(|_, _| Complex::new(d.sample(&mut rng), d.sample(&mut rng)));
let p = project_to_su2(r);
assert_matrix_is_su_2!(p, EPSILON);
}
}
#[test]
fn random_su2_t() {
let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
for _ in 0_u32..100_u32 {
let m = random_su2(&mut rng);
assert_matrix_is_su_2!(m, EPSILON);
}
for _ in 0_u32..100_u32 {
let m = random_su2(&mut rng);
assert!(is_matrix_su2(&m, EPSILON));
}
for _ in 0_u32..100_u32 {
let m = random_su2(&mut rng) * Complex::new(1.5_f64, 0.7_f64);
assert!(!is_matrix_su2(&m, EPSILON));
}
}
}