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

use std::ops::{Mul, Add};
use nalgebra::{DefaultAllocator, MatrixMN, DimName, U2, U4, Matrix2, Matrix4, Matrix};
use nalgebra::allocator::Allocator;

use {Ket, Complex, SQRT_2_INVERSE};

/// Outer product in bra-ket notation, as well used as a linear operatior.
/// You can multiply it by ket to put it in the observable state defined by the operator.
#[derive(Clone, Debug, PartialEq)]
pub struct Outer<D: DimName>(pub(crate) MatrixMN<Complex, D, D>)
    where DefaultAllocator: Allocator<Complex, D, D>;

impl<D: DimName> Mul<Ket<D>> for Outer<D>
    where DefaultAllocator: Allocator<Complex, D> + Allocator<Complex, D, D>
{
    type Output = Ket<D>;

    fn mul(self, other: Ket<D>) -> Self::Output {
        Ket(self.0 * other.0)
    }
}

impl<D: DimName> ::std::fmt::Display for Outer<D>
    where DefaultAllocator: Allocator<Complex, D, D>,
    DefaultAllocator: Allocator<usize, D, D>
{
    fn fmt(&self, f: &mut ::std::fmt::Formatter) -> ::std::fmt::Result {
        write!(f, "{}", self.0)
    }
}

impl<D: DimName> Add<Outer<D>> for Outer<D>
    where DefaultAllocator: Allocator<Complex, D, D>
{
    type Output = Self;

    fn add(self, other: Outer<D>) -> Self::Output {
        Outer(self.0 + other.0)
    }
}

impl<D: DimName> Outer<D>
    where DefaultAllocator: Allocator<Complex, D, D>
{
    /// Deconstruct the outer product returning the inner matrix.
    pub fn into_matrix(self) -> MatrixMN<Complex, D, D> {
        self.0
    }

    /// H2 (2 dim hadamard) operator
    pub fn h2() -> Outer<U2> {
        Outer::<U2>(Matrix2::<Complex>::new(1.0.into(), 1.0.into(), 1.0.into(), (-1.0).into())) * SQRT_2_INVERSE
    }

    /// Z2 (2 dim Z-gate) operator
    pub fn z2() -> Outer<U2> {
        Outer::<U2>(Matrix2::<Complex>::new(1.0.into(), 0.0.into(), 0.0.into(), (-1.0).into()))
    }

    /// N2 (2 dim N-gate) operator
    pub fn n2() -> Outer<U2> {
        Outer::<U2>(Matrix2::<Complex>::new(0.0.into(), 1.0.into(), 1.0.into(), 0.0.into()))
    }

    /// CNOT (2*2 dim) operator
    pub fn cnot() -> Outer<U4> {
        Outer::<U4>(
            Matrix4::<Complex>::new(
                1.0.into(), 0.0.into(), 0.0.into(), 0.0.into(),
                0.0.into(), 1.0.into(), 0.0.into(), 0.0.into(),
                0.0.into(), 0.0.into(), 0.0.into(), 1.0.into(),
                0.0.into(), 0.0.into(), 1.0.into(), 0.0.into(),
            )
        )
    }

    /// Quantum Fourier Transform (QFT) matrix operator
    pub fn qft() -> Outer<D> {
        let mut matrix: MatrixMN<Complex, D, D> = Matrix::zeros_generic(D::name(), D::name());
        let dim = D::name().value();

        let n = dim as f64;

        let coef = (Complex::from(1.0) / (n as f64)).sqrt();

        for i in 0..dim {
            for j in 0..dim {
                let power = (Complex::from(2.0) * ::std::f64::consts::PI * Complex::i() / n)
                    * (i as f64) * (j as f64);

                *matrix.get_mut((i, j)).expect("(i, j) in (dim, dim) range") = power.exp() * coef;
            }
        }

        Outer(matrix)
    }
}

impl<D: DimName> From<MatrixMN<Complex, D, D>> for Outer<D>
    where DefaultAllocator: Allocator<Complex, D, D>
{
    fn from(v: MatrixMN<Complex, D, D>) -> Self {
        Outer(v)
    }
}

impl<D: DimName> Mul<f64> for Outer<D>
    where DefaultAllocator: Allocator<Complex, D, D>
{
    type Output = Self;

    fn mul(self, other: f64) -> Self::Output {
        Outer(self.0 * Complex::new(other, 0.0))
    }
}