Quantum gates definition and decomposition

#
# Copyright (c) 2020- Beit, Beit.Tech, Beit.Inc
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
import numpy as np
import numpy.linalg as la
import warnings

from numpy import allclose, array, cos, exp, isclose, kron, log2, matmul, pi, set_printoptions, sin, sqrt
from functools import reduce

set_printoptions(
    formatter={
        'int': ('{:+2}').format,
        'float': ('{:+0.2f}').format,
        'complexfloat': ('{:+.2f}').format
    },
    linewidth=120
)
def kron(*xs):
    return reduce(np.kron, xs, array(1))

def kron_str(*xs):
    return '(' + ', '.join(xs) + ')'

def dot(*xs):
    return reduce(np.dot, xs, array(1))

def tr(matrix):
    return np.transpose(matrix)

def cj(matrix):
    return np.conjugate(matrix)

def inv(matrix):
    return la.inv(matrix)

def dagger(matrix):
    return cj(tr(matrix))

def is_real(c, atol=0.001):
    return isclose(complex(c).imag, 0, atol=atol)

def is_integer(f):
    return (isinstance(f, int) or isinstance(f, float)) and isclose(int(f), f)

def is_unitary(matrix, atol=0.001):
    return allclose(matrix @ dagger(matrix), np.eye(matrix.shape[0]), atol=atol)

def is_hermitian(matrix, atol=0.001):
    return allclose(matrix, tr(cj(matrix)), atol=atol)
    
def show(matrix, precision=2):
    set_printoptions(
        formatter={
            'int': ('{: ' + str(precision) + '}').format,
            'float': ('{: 0.' + str(precision) + 'f}').format,
            'complexfloat': ('{: .' + str(precision) + 'f}').format
        },
        linewidth=120
    )

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")

        if np.all(np.vectorize(is_real)(matrix)):
            matrix = matrix.astype(float)

        if np.all(np.vectorize(is_integer)(matrix)):
            matrix = matrix.astype(int)

    print(matrix)
    
    
class program:
    def __init__(self, *moments):
        moments = [kron(*moment) if isinstance(moment, tuple) else moment for moment in moments]
        self._matrix = reduce(matmul, moments[::-1])
        
    def _get_psi(self, psi):
        if psi is None:
            size = self._matrix.shape[0]
            return [[1]] + [[0] for _ in range(size-1)]
        elif isinstance(psi, tuple):
            return kron(*psi)
        else:
            return psi

    def unitary(self):
        return self._matrix
        
    def run(self, psi=None) -> array:
        return self.unitary() @ self._get_psi(psi)
    
    def show(self, psi=None, precision=2):
        set_printoptions(
            formatter={
                'int': ('{: ' + str(precision) + '}').format,
                'float': ('{: 0.' + str(precision) + 'f}').format,
                'complexfloat': ('{: .' + str(precision) + 'f}').format
            },
            linewidth=120
        )

        matrix = self.unitary()
        psi0 = self._get_psi(psi)
        psiN = matrix @ psi0
        size = self._matrix.shape[0]

        with warnings.catch_warnings():
            warnings.simplefilter("ignore")

            if np.all(np.vectorize(is_real)(matrix)):
                matrix = matrix.astype(float)

            if np.all(np.vectorize(is_real)(psi0)):
                psi0 = psi0.astype(float)

            if np.all(np.vectorize(is_real)(psiN)):
                psiN = psiN.astype(float)            
        
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")

            if np.all(np.vectorize(is_integer)(matrix)):
                matrix = matrix.astype(int)

            if np.all(np.vectorize(is_integer)(psi0)):
                psi0 = psi0.astype(int)

            if np.all(np.vectorize(is_integer)(psiN)):
                psiN = psiN.astype(int) 
        
        for (id_, (row, psi0_i, psiN_i)) in enumerate(zip(matrix, psi0, psiN)):
            print((
                '{}' + (' @ ' if id_ == size//2 else '   ') +
                '[{:.' + str(precision) + 'f}]' + (' = ' if id_ == size//2 else '   ') +
                '[{: .' + str(precision) + 'f}]' + 
                '  [{:3.0f}%]' + 
                '  |{:0' + str(int(log2(matrix.shape[0]))) + 'b}>'
            ).format(
                row, 
                psi0_i[0], 
                psiN_i[0], 
                100*abs(psiN_i[0])**2,
                id_
            )
        )

Quantum Gates Definition

#
## States
#

# ket zero state
q0 = array([ # |0>
    [1], 
    [0]
])

# ket one state
q1 = array([  # X |0>
    [0], 
    [1]
])

# ket plus state
qPlus = 0.5**0.5 * array([ # H |0>
    [1 ],
    [1j]
])

# ket minus state
qMinus = 0.5**0.5 * array([  # H X |0>
    [ 1 ],
    [-1j]
])

# ket plus Y state
qPlusY = 0.5**0.5 * array([  # S H |0>
    [1],
    [1]
])

# ket minus Y state
qMinusY = 0.5**0.5 * array([  # Sdg H |0>
    [ 1],
    [-1]
])
# I - Id - Iden - Identity gate
I = array([
    [1, 0], 
    [0, 1]
])

# Two qubit identity gate
I4 = array([
    [1, 0, 0, 0], 
    [0, 1, 0, 0], 
    [0, 0, 1, 0], 
    [0, 0, 0, 1]
])
# X - Pauli-X - Not - gate matrix
X = array([
    [0, 1], 
    [1, 0]
])

Not = X

# CX - Controlled X - CNot - Controlled Not - gate matrix
CX = array([
    [1, 0, 0, 0], 
    [0, 1, 0, 0], 
    [0, 0, 0, 1], 
    [0, 0, 1, 0]
])

CNot = CX

# XC - X Controlled - Not C - Not Controlled - gate matrix
XC = array([
    [1, 0, 0, 0], 
    [0, 0, 0, 1], 
    [0, 0, 1, 0], 
    [0, 1, 0, 0]
])

NotC = XC

# CCX - Controlled CX - Controlled Controlled X - CCNot - Controlled CNot - Controlled Controlled Not - gate matrix
CCX = array([
    [1, 0, 0, 0, 0, 0, 0, 0],
    [0, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 1, 0]
])

CCNot = CCX

# Square root of Pauli-X - square root of X - square root of Not - gate matrix
XSR = 0.5 * array([
    [1+1j, 1-1j],
    [1-1j, 1+1j]
])

# Controlled square root of Pauli-X - Controlled square root of X - Controlled square root of Not - gate matrix
CXsr = 0.5 * array([
    [2, 0,    0,    0],
    [0, 2,    0,    0],
    [0, 0, 1+1j, 1-1j],
    [0, 0, 1-1j, 1+1j]
])
# Y - Pauli-Y - gate matrix
Y = array([
    [0, -1j], 
    [1j,  0]
])

# CY - Controlled Y - gate matrix
CY = array([
    [1, 0,  0,   0], 
    [0, 1,  0,   0], 
    [0, 0,  0, -1j], 
    [0, 0, 1j,   0]
])

# YC - Y Controlled - gate matrix
YC = array([
    [1,  0, 0,   0], 
    [0,  0, 0, -1j], 
    [0,  0, 1,   0], 
    [0, 1j, 0,   0]
])

# CCY - Controlled Controlled Y - gate matrix
CCY = array([
    [1, 0, 0, 0, 0, 0,  0,   0],
    [0, 1, 0, 0, 0, 0,  0,   0],
    [0, 0, 1, 0, 0, 0,  0,   0],
    [0, 0, 0, 1, 0, 0,  0,   0],
    [0, 0, 0, 0, 1, 0,  0,   0],
    [0, 0, 0, 0, 0, 1,  0,   0],
    [0, 0, 0, 0, 0, 0,  0, -1j],
    [0, 0, 0, 0, 0, 0, 1j,   0]
])
# Z - Pauli-Z - gate matrix
Z = array([
    [1,  0], 
    [0, -1]
])

# CZ - Controlled Z - gate matrix
CZ = array([
    [1, 0, 0,  0], 
    [0, 1, 0,  0], 
    [0, 0, 1,  0], 
    [0, 0, 0, -1]
])

# ZC - Z Controlled - gate matrix
ZC = CZ

# CCZ - Controlled Controlled Z - gate matrix
CCZ = array([
    [1, 0, 0, 0, 0, 0, 0,  0],
    [0, 1, 0, 0, 0, 0, 0,  0],
    [0, 0, 1, 0, 0, 0, 0,  0],
    [0, 0, 0, 1, 0, 0, 0,  0],
    [0, 0, 0, 0, 1, 0, 0,  0],
    [0, 0, 0, 0, 0, 1, 0,  0],
    [0, 0, 0, 0, 0, 0, 1,  0],
    [0, 0, 0, 0, 0, 0, 0, -1]
])
# H - Hadamard - gate matrix
H = 0.5**0.5 * \
    array([
        [1,  1], 
        [1, -1]
])

# CH - Controlled Hadamard - gate matrix
CH = array([
    [1, 0,        0,         0], 
    [0, 1,        0,         0], 
    [0, 0, 0.5**0.5,  0.5**0.5], 
    [0, 0, 0.5**0.5, -0.5**0.5]
])

# HC - Hadamard Controlled - gate matrix
HC = array([
    [1,        0, 0,         0], 
    [0, 0.5**0.5, 0,  0.5**0.5], 
    [0,        0, 1,         0], 
    [0, 0.5**0.5, 0, -0.5**0.5]
])

# CCH - Controlled Controlled Hadamard - gate matrix
CCH = array([
    [1, 0, 0, 0, 0, 0,        0,         0],
    [0, 1, 0, 0, 0, 0,        0,         0],
    [0, 0, 1, 0, 0, 0,        0,         0],
    [0, 0, 0, 1, 0, 0,        0,         0],
    [0, 0, 0, 0, 1, 0,        0,         0],
    [0, 0, 0, 0, 0, 1,        0,         0],
    [0, 0, 0, 0, 0, 0, 0.5**0.5,  0.5**0.5],
    [0, 0, 0, 0, 0, 0, 0.5**0.5, -0.5**0.5]
])
# Swap gate matrix
Swap = array([
    [1, 0, 0, 0], 
    [0, 0, 1, 0], 
    [0, 1, 0, 0], 
    [0, 0, 0, 1]
])

# Controlled Swap - Fredkin - gate matrix
CSwap = array([
    [1, 0, 0, 0, 0, 0, 0, 0],
    [0, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 1]
])

Fredkin = CSwap
# U3 - Universal gate - matrix
def U3(theta, phi, lambda_):
    return array([ # U3(theta, phi, lambda) = RZ(phi+3pi) RX(pi/2) RZ(theta+pi) RX(pi/2) RZ(lambda)
        [              cos(theta/2), -exp(1j *      lambda_ ) * sin(theta/2)],
        [exp(1j*phi) * sin(theta/2),  exp(1j * (phi+lambda_)) * cos(theta/2)]
    ])

# CU3 - Controlled Universal gate - matrix
def CU3(theta, phi, lambda_):
    return array([
        [1, 0,                          0,                                       0],
        [0, 1,                          0,                                       0],
        [0, 0,               cos(theta/2), -exp(1j *      lambda_ ) * sin(theta/2)],
        [0, 0, exp(1j*phi) * sin(theta/2),  exp(1j * (phi+lambda_)) * cos(theta/2)]
    ])

# U2 - Universal gate - matrix
def U2(phi, lambda_):
    return array([
        [              cos(pi/2/2), -exp(1j *      lambda_ ) * sin(pi/2/2)],
        [exp(1j*phi) * sin(pi/2/2),  exp(1j * (phi+lambda_)) * cos(pi/2/2)]
    ])

# U1 - Universal gate - matrix
def U1(lambda_):
    return array([
        [1,                  0],
        [0,  exp(1j * lambda_)]
    ])
# S gate matrix
S = array([
    [1,  0], 
    [0, 1j]
])

# CS - Controlled S - gate matrix
CS = array([
    [1, 0, 0,  0],
    [0, 1, 0,  0],
    [0, 0, 1,  0], 
    [0, 0, 0, 1j]
])

# Sdg - S dagger - gate matrix
Sdg = array([
    [1,   0], 
    [0, -1j]
])

# CSdg - Controlled Sdg - gate matrix
CSdg = array([
    [1, 0, 0,   0],
    [0, 1, 0,   0],
    [0, 0, 1,   0], 
    [0, 0, 0, -1j]
])

# T gate matrix
T = array([
    [1,              0], 
    [0, exp(2*pi*1j/8)]
])

# CT - Controlled T - gate matrix
CT = array([
    [1, 0, 0,              0],
    [0, 1, 0,              0],
    [0, 0, 1,              0], 
    [0, 0, 0, exp(2*pi*1j/8)]
])

# Tdg - T dagger - gate matrix
Tdg = array([
    [1,               0], 
    [0, exp(-2*pi*1j/8)]
])

# CTdg - Controlled Tdg - gate matrix
CTdg = array([
    [1, 0, 0,               0],
    [0, 1, 0,               0],
    [0, 0, 1,               0], 
    [0, 0, 0, exp(-2*pi*1j/8)]
])
# RX - Rotation around X - gate matrix
def RX(theta):
    return array([
        [      cos(theta/2), -1j * sin(theta/2)],
        [-1j * sin(theta/2),       cos(theta/2)]
    ])

# CRX - Controlled RX - Controlled Rotation around X - gate matrix
def CRX(theta):
    return array([
        [1, 0,                  0,                  0],
        [0, 1,                  0,                  0],
        [0, 0,       cos(theta/2), -1j * sin(theta/2)],
        [0, 0, -1j * sin(theta/2),       cos(theta/2)]
    ])

# RY gate matrix
def RY(theta):
    return array([
        [cos(theta/2), -sin(theta/2)],
        [sin(theta/2),  cos(theta/2)]
    ])

# CRY - Controlled RY - Controlled Rotation around Y - gate matrix
def CRY(theta):
    return array([
        [1, 0,            0,             0],
        [0, 1,            0,             0],
        [0, 0, cos(theta/2), -sin(theta/2)],
        [0, 0, sin(theta/2),  cos(theta/2)]
    ])

# RZ gate matrix
def RZ(theta):
    return array([
        [exp(-1j * theta/2),                 0],
        [                 0, exp(1j * theta/2)]
    ])

# CRZ - Controll RZ - Controlled Rotation Z - gate matrix
def RZ(theta):
    return array([
        [1, 0,                  0,                0],
        [0, 1,                  0,                0],
        [0, 0, exp(-1j * theta/2),                0],
        [0, 0,                  0, exp(1j * theta/2)]
    ])

# R gate matrix
def R(theta, phi): # Rotation theta around the cos(psi)x + sin(psi)y axis
    return array([
        [                      cos(theta/2), -1j * exp(-1j * phi) * sin(theta/2)],
        [-1j * exp(1j * phi) * sin(theta/2),                        cos(theta/2)]
    ])

# CR - Controlled R - gate matrix
def R(theta, phi): # Rotation theta around the cos(psi)x + sin(psi)y axis
    return array([ 
        [1, 0,                                  0,                                   0],
        [0, 1,                                  0,                                   0],
        [0, 0,                       cos(theta/2), -1j * exp(-1j * phi) * sin(theta/2)],
        [0, 0, -1j * exp(1j * phi) * sin(theta/2),                        cos(theta/2)]
    ])

# Phase gate matrix
def Phase(theta):
    return array([
        [1,               0],
        [0, exp(1j * theta)]
    ])

# CPhase - Controlled Phase - gate matrix
def Phase(theta):
    return array([
        [1, 0, 0,               0],
        [0, 1, 0,               0],
        [0, 0, 1,               0],
        [0, 0, 0, exp(1j * theta)]
    ])

Quantum Gates Decomposition

# X - Pauli-X - gate decomposition
X = U3(pi, 0, pi)

# CX - Controlled X - CNot - Controlled Not - gate decomposition
CX = program(
    (I, H),
    CZ,
    (I, H)
).unitary()

# X Controlled gate decomposition
XC = program(
    (H, H),
    CX,
    (H, H)
).unitary()

# CCX - Controlled CX - Controlled Controlled X - CCNot - Controlled CNot - Controlled Controlled Not - gate decomposition
CCX = program(
    (I, CX),
    (Swap, I), CCX, (Swap, I),
    (I, CX)
).unitary()

CCX = program(
    (T, I, I),
    (I, I, H),
    (CX, I),
    (I, Tdg, I),
    (CX, I),
    (I, T, I),
    (I, CX),
    (I, I, Tdg),
    (Swap, I), (I, CX), (Swap, I),  # (C, I, X))
    (I, I, T),
    (I, CX),
    (I, I, Tdg),
    (Swap, I), (I, CX), (Swap, I),  # (C, I, X)
    (I, I, T),
    (I, I, H)
).unitary()

CCX = program(
  (I, I, H),
  (I, CX),
  (I, I, Tdg),
  (Swap, I), (I, CX), (Swap, I),  # (C, I, X)
  (I, I, T),
  (I, CX),
  (I, I, Tdg),
  (Swap, I), (I, CX), (Swap, I),  # (C, I, X)
  (I, T, I),
  (I, I, T),
  (I, I, H),
  (CX, I),
  (T, I, I),
  (I, Tdg, I),
  (CX, I)
).unitary()

# Square root of Pauli-X - Square root of X - Square root of Not - gate decomposition
XSR = U3(0.5*pi, -0.5*pi, 0.5*pi)

# Controled square root of Pauli-X - controlled square root of X - controlled square root of Not
CXSR = program(
    (I, H),
    (T, T),
    CX,
    (I, Tdg),
    CX,
    (I, H)
).unitary()

# Controled square root of Pauli-X - controlled square root of X - controlled square root of Not
CXSRdg = program(
    (I, H),
    CX,
    (I, T),
    CX,
    (Tdg, Tdg),
    (I, H)
).unitary()

# Margolus gate decomposition
Margolus = program(
    (I, I, RY(pi/4)),
    (I, CX),
    (I, I, RY(pi/4)),
    (Swap, I), (I, CX), (Swap, I),  # (C, I, X)
    (I, I, RY(-pi/4)),
    (I, CX),
    (I, I, RY(-pi/4))
).unitary()
# Y - Pauli-Y - gate decomposition
Y = U3(pi, pi/2, pi/2)

# CY - Controlled Y - gate decomposition
CY = program(
    (I, Sdg),
    CX,
    (I, S)
).unitary()
# Z - Pauli Z - gate decomposition
Z = U3(0, 0, pi)

# CZ - Controlled Z - gate decomposition
CZ = CU3(0, 0, pi)

CZ = program(
    (I, H),
    CX,
    (I, H)
).unitary()
# H - Hadamard - gate decomposition
H = U3(pi/2, 0, pi)

H = 1j * program(
    RY(pi/2),
    RX(pi)
).unitary()

# CH - Controlled H - Controlled Hadamard - gate decomposition
CH = exp(-1j*pi/4) * program(  # #11
    (I, H),
    (I, Sdg),
    CX,
    (I, H),
    (I, T),
    CX,
    (I, T),
    (I, H),
    (I, S),
    (I, X),
    (S, I)
).unitary()

# CH - Controlled H - Controlled Hadamard - gate decomposition
CH = CU3(pi/2, 0, pi)
# S gate decomposition
S = U3(0, 0, pi/2)

# Sdg - S dagger - gate decomposition
Sdg = U3(0, 0, -pi/2)

# T gate decomposition
T = U3(0, 0, pi/4)

# Tdg - T dagger - gate decomposition
Tdg = U3(0, 0, -pi/4)
# U2 - Universal - gate decomposition
def U2(phi, lambda_): 
    return U3(pi/2, phi, lambda_)

# U1 - Universal - gate decomposition
def U1(lambda_):
    return U3(0, 0, lambda_)
# Phase gate decomposition
Phase = U1
    
# CPhase - Controlled Phase - gate decomposition
def CPhase(theta):
    return program(  # #5
        (I, U3(0, 0, theta/2)),
        CX,
        (I, U3(0, 0, -theta/2)),
        CX,
        (U3(0, 0, theta/2), I)
    ).unitary()
# CS - Controlled S - gate decomposition
CS = CU3(0, 0, pi/2)

# CSdg - Controlled Sdg - Controlled S dagger - gate decomposition
CSdg = CU3(0, 0, -pi/2)

# CT - Controlled T - gate decomposition
CT = CU3(0, 0, pi/4)

# CTdg - Controlled Tdg - Controlled T dagger - gate decomposition
CTdg = CU3(0, 0, -pi/4)
# Swap gate decomposition
Swap = program(
    CX,
    XC,
    CX
).unitary()

# Controlled Swap - Fredkin - gate decomposition
CSwap = program(  # #7
    (I, Swap), (CCX), (I, Swap),  # (C, X, C)
    CCX,
    (I, Swap), (CCX), (I, Swap)  # (C, X, C)
).unitary()

CSwap = program(  # #9
    (I, XC),
    (I, CXSR),
    (Swap, I), (I, CXSR), (Swap, I),  # (C, I, XSR)
    (CX, I),
    (I, CXSRdg),
    (I, XC),
    (CX, I)
).unitary()
# RX gate decomposition
def RX(theta):
    return U3(theta, -pi/2, pi/2)

# RY gate decomposition
def RY(theta):
    return U3(theta, 0, 0)

# RZ gate decomposition
def RZ(theta): # Quasar Implementation
    return program(  # #4
        X,
        U3(0, 0, -theta/2),
        X,
        U3(0, 0, theta/2)
    ).unitary()

Two qubits oracles

# 00 ->  00
# 01 -> -01
# 10 ->  10
# 11 -> -11
program(
    (I, Z)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0 -1  0  0]   [0.50]   [-0.50]  [ 25%]  |01>
[ 0  0  1  0] @ [0.50] = [ 0.50]  [ 25%]  |10>
[ 0  0  0 -1]   [0.50]   [-0.50]  [ 25%]  |11>
# 00 ->  00
# 01 ->  01
# 10 -> -10
# 11 -> -11
program(
    (Z, I)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0  1  0  0]   [0.50]   [ 0.50]  [ 25%]  |01>
[ 0  0 -1  0] @ [0.50] = [-0.50]  [ 25%]  |10>
[ 0  0  0 -1]   [0.50]   [-0.50]  [ 25%]  |11>
# 00 ->  00
# 01 -> -01
# 10 -> -10
# 11 ->  11
program(
    (Z, Z)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0 -1  0  0]   [0.50]   [-0.50]  [ 25%]  |01>
[ 0  0 -1  0] @ [0.50] = [-0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# 00 -> -00
# 01 ->  01
# 10 ->  10
# 11 ->  11
program(
    (X, X),
    (CZ),
    (X, X)
).show(kron(H, H))
[-1  0  0  0]   [0.50]   [-0.50]  [ 25%]  |00>
[ 0  1  0  0]   [0.50]   [ 0.50]  [ 25%]  |01>
[ 0  0  1  0] @ [0.50] = [ 0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# 00 ->  00
# 01 -> -01
# 10 ->  10
# 11 ->  11
program(
    (X, I),
    (CZ),
    (X, I)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0 -1  0  0]   [0.50]   [-0.50]  [ 25%]  |01>
[ 0  0  1  0] @ [0.50] = [ 0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# 00 ->  00
# 01 -> -01
# 10 ->  10
# 11 ->  11
program(  # #5
    (I, S),
    (I, H),
    CX,
    (I, H),
    (I, S)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0 -1  0  0]   [0.50]   [-0.50]  [ 25%]  |01>
[ 0  0  1  0] @ [0.50] = [ 0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# 00 ->  00
# 01 ->  01
# 10 -> -10
# 11 ->  11
program(
    (I, X),
    (CZ),
    (I, X)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0  1  0  0]   [0.50]   [ 0.50]  [ 25%]  |01>
[ 0  0 -1  0] @ [0.50] = [-0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# 00 ->  00
# 01 ->  01
# 10 -> -10
# 11 ->  11
program(  # #5
    (S, I),
    (I, H),
    CX,
    (I, H),
    (S, I)
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0  1  0  0]   [0.50]   [ 0.50]  [ 25%]  |01>
[ 0  0 -1  0] @ [0.50] = [-0.50]  [ 25%]  |10>
[ 0  0  0  1]   [0.50]   [ 0.50]  [ 25%]  |11>
# Two qubits AND oracle
# 00 ->  00
# 01 ->  01
# 10 ->  10
# 11 -> -11
program(
    CZ
).show(kron(H, H))
[ 1  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |00>
[ 0  1  0  0]   [0.50]   [ 0.50]  [ 25%]  |01>
[ 0  0  1  0] @ [0.50] = [ 0.50]  [ 25%]  |10>
[ 0  0  0 -1]   [0.50]   [-0.50]  [ 25%]  |11>

Two qubits with ancilla oracles

# Two qubits AND oracle with ancilla
# 00 ->  00
# 01 ->  01
# 10 ->  10
# 11 -> -11
program(
    CCX,
    (I, I, Z),
    CCX  # uncompute
).show(kron(H, H, I))
[ 1  0  0  0  0  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |000>
[ 0 -1  0  0  0  0  0  0]   [0.00]   [-0.00]  [  0%]  |001>
[ 0  0  1  0  0  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |010>
[ 0  0  0 -1  0  0  0  0]   [0.00]   [ 0.00]  [  0%]  |011>
[ 0  0  0  0  1  0  0  0] @ [0.50] = [ 0.50]  [ 25%]  |100>
[ 0  0  0  0  0 -1  0  0]   [0.00]   [-0.00]  [  0%]  |101>
[ 0  0  0  0  0  0 -1  0]   [0.50]   [-0.50]  [ 25%]  |110>
[ 0  0  0  0  0  0  0  1]   [0.00]   [ 0.00]  [  0%]  |111>
# Two qubits NAND oracle with ancilla - NAND == NOT(AND)
# 00 -> -00
# 01 -> -01
# 10 -> -10
# 11 ->  11
program(  # #5
    CCX, 
    (I, I, X),
    (I, I, Z),
    (I, I, X), 
    CCX  # uncompute
).show(kron(H, H, I), precision=0)
[-1  0  0  0  0  0  0  0]   [0]   [-1]  [ 25%]  |000>
[ 0  1  0  0  0  0  0  0]   [0]   [ 0]  [  0%]  |001>
[ 0  0 -1  0  0  0  0  0]   [0]   [-1]  [ 25%]  |010>
[ 0  0  0  1  0  0  0  0]   [0]   [ 0]  [  0%]  |011>
[ 0  0  0  0 -1  0  0  0] @ [0] = [-1]  [ 25%]  |100>
[ 0  0  0  0  0  1  0  0]   [0]   [ 0]  [  0%]  |101>
[ 0  0  0  0  0  0  1  0]   [1]   [ 1]  [ 25%]  |110>
[ 0  0  0  0  0  0  0 -1]   [0]   [ 0]  [  0%]  |111>
# Two qubits OR oracle with ancilla - OR == NOT(ZERO)
# 00 ->  00
# 01 -> -01
# 10 -> -10
# 11 -> -11
program(  # #9
    (X, X, I),
    CCX,
    (X, X, I),
    (I, I, X),
    (I, I, Z),
    (I, I, X),
    (X, X, I),
    CCX,
    (X, X, I)
).show(kron(H, H, I))
[ 1  0  0  0  0  0  0  0]   [0.50]   [ 0.50]  [ 25%]  |000>
[ 0 -1  0  0  0  0  0  0]   [0.00]   [ 0.00]  [  0%]  |001>
[ 0  0 -1  0  0  0  0  0]   [0.50]   [-0.50]  [ 25%]  |010>
[ 0  0  0  1  0  0  0  0]   [0.00]   [ 0.00]  [  0%]  |011>
[ 0  0  0  0 -1  0  0  0] @ [0.50] = [-0.50]  [ 25%]  |100>
[ 0  0  0  0  0  1  0  0]   [0.00]   [ 0.00]  [  0%]  |101>
[ 0  0  0  0  0  0 -1  0]   [0.50]   [-0.50]  [ 25%]  |110>
[ 0  0  0  0  0  0  0  1]   [0.00]   [ 0.00]  [  0%]  |111>

Circuit Generators

#
## Swap any two qubits on linear topology
#

def _gen_swap(qubits, i, j, iden, swap, kron):
    assert(i < qubits)
    assert(j < qubits)

    result = []

    if (i == j):
        return result

    if (i > j):
        return _gen_swap(qubits, j, i, iden, swap, kron)

    for k in range(i, j):
        cmd = (k * [iden]) + [swap] + ((qubits - k - 2) * [iden])
        result.append(kron(*cmd))

    for k in range(j - 2, i - 1, -1):
        cmd = (k * [iden]) + [swap] + ((qubits - k - 2) * [iden])
        result.append(kron(*cmd))

    return result

def get_swap(qubits, i, j):
    moments = _gen_swap(qubits, i, j, I, Swap, kron)
    
    if len(moments) == 0:
        return np.eye(2**qubits)
    
    return dot(*moments)

def gen_swap(qubits, i, j):
    return _gen_swap(qubits, i, j, 'I', 'Swap', kron_str)
    
assert(gen_swap(4, 0, 3) == ['(Swap, I, I)', '(I, Swap, I)', '(I, I, Swap)', '(I, Swap, I)', '(Swap, I, I)'])
#
## Controlled quantum gate on any qubits on linear topology
#

def _gen_controlled_1(qubits, control, target, ccmd, iden, swap, kron):
    result = []
    
    (left, pivot) = (min(control, target), max(control, target))
    
    if left != pivot - 1:
        result.extend(_gen_swap(qubits, left, pivot - 1, iden, swap, kron))
            
    if left != control:
        result.extend(_gen_swap(qubits, pivot - 1, pivot, iden, swap, kron))
                      
    cmd = (pivot - 1) * [iden] + [ccmd] + (qubits - pivot - 1) * [iden]
    cmd = kron(*cmd)

    return result + [cmd] + result[::-1]
    
def get_controlled_1(qubits, control, target, ccmd):
    moments = _gen_controlled_1(qubits, control, target, ccmd, I, Swap, kron)
    
    if len(moments) == 0:
        return np.eye(2**qubits)
    
    return dot(*moments)

def gen_controlled_1(qubits, control, target, ccmd):
    return _gen_controlled_1(qubits, control, target, ccmd, 'I', 'Swap', kron_str)
    
assert(gen_controlled_1(2, 1, 0, 'CX') == ['(Swap)', '(CX)', '(Swap)'])
#
## Controlled controlled quantum gate on any qubits on linear topology
#

def _gen_controlled_2(qubits, control1, control2, target, ccmd, iden, swap, kron):
    result = []
    
    (left, pivot, right) = (sorted([control1, control2, target]))
    
    if left != pivot - 1:
        result.extend(_gen_swap(qubits, left, pivot - 1, iden, swap, kron))

    if right != pivot + 1:
        result.extend(_gen_swap(qubits, pivot + 1, right, iden, swap, kron))
        
    if left != control1:
        if pivot == control1:
            result.extend(_gen_swap(qubits, pivot - 1, pivot, iden, swap, kron))
        else:
            result.extend(_gen_swap(qubits, pivot - 1, pivot + 1, iden, swap, kron))
            
    if pivot != control2:
        result.extend(_gen_swap(qubits, pivot, pivot + 1, iden, swap, kron))
                      
    cmd = (pivot - 1) * [iden] + [ccmd] + (qubits - pivot - 2) * [iden]
    cmd = kron(*(cmd))
    
    return result + [cmd] + result[:: -1]
    
def get_controlled_2(qubits, control1, control2, target, ccmd):
    moments = _gen_controlled_2(qubits, control1, control2, target, ccmd, I, Swap, kron)
    
    if len(moments) == 0:
        return np.eye(2**qubits)
    
    return dot(*moments)

def gen_controlled_2(qubits, control1, control2, target, ccmd):
    return _gen_controlled_2(qubits, control1, control2, target, ccmd, 'I', 'Swap', kron_str)

assert(gen_controlled_2(3, 0, 2, 1, 'CCX') == ['(I, Swap)', '(CCX)', '(I, Swap)'])
#
## Controlled Controlled Quantum Gate on Linear Topology
#

class _Controlled:
    pass

C = _Controlled()

def is_controlled(c):
    return isinstance(c, _Controlled)

def is_identity(i):
    return (i.shape == (2, 2)) and (allclose(i, I))

def is_gate(g):
    return not(is_controlled(g) or is_identity(g))

def get_index(xs, pred):
    if len(xs) == 0:
        return -1
    
    if pred(xs[0]):
        return 0
    
    return 1 + get_index(xs[1:], pred)

def get_hash(obj):
    return hash(str(obj))

def _decompose(ms, cu_mapping, ccu_mapping, gen_controlled_1_cmd, gen_controlled_2_cmd, iden, swap, kron): 
    controlled_counter = sum(map(is_controlled, ms))
    assert(controlled_counter < 3)
    
    gate_counter = sum(map(is_gate, ms))
    if gate_counter != 1:
        return [kron(*ms)]
    
    if controlled_counter == 0:
        return [kron(*ms)]
    elif controlled_counter == 1:
        control = get_index(ms, is_controlled)
        target = get_index(ms, is_gate)
        ccmd = cu_mapping[get_hash(ms[target])]
        return gen_controlled_1_cmd(len(ms), control, target, ccmd, iden, swap, kron)
    else:
        control_1 = get_index(ms, is_controlled)
        control_2 = (control_1 + 1) + get_index(ms[control_1 + 1: ], is_controlled)
        target = get_index(ms, is_gate)
        ccmd = ccu_mapping[get_hash(ms[target])]        
        return gen_controlled_2_cmd(len(ms), control_1, control_2, target, ccmd, iden, swap, kron)

def gen_decomposition(ms):
    cu_mapping = {
        get_hash(X)  : 'CX',
        get_hash(Y)  : 'CY',
        get_hash(Z)  : 'CZ',
        get_hash(H)  : 'CH',
        get_hash(S)  : 'CS',
        get_hash(Sdg): 'CSdg',
        get_hash(T)  : 'CT',
        get_hash(Tdg): 'CSdg',
    }
    
    ccu_mapping = {
        get_hash(X)  : 'CCX',
    }
    
    return _decompose(
        ms, 
        cu_mapping, 
        ccu_mapping, 
        _gen_controlled_1, 
        _gen_controlled_2, 
        'I', 
        'Swap', 
        kron_str
    )
    
def get_decomposition(ms):
    cu_mapping = {
        get_hash(X)  : CX,
        get_hash(Y)  : CY,
        get_hash(Z)  : CZ,
        get_hash(H)  : CH,
        get_hash(S)  : CS,
        get_hash(Sdg): CSdg,
        get_hash(T)  : CT,
        get_hash(Tdg): CSdg,
    }
    
    ccu_mapping = {
        get_hash(X)  : CCX,
    }
    
    return dot(*_decompose(
        ms, 
        cu_mapping, 
        ccu_mapping, 
        _gen_controlled_1, 
        _gen_controlled_2, 
        I, 
        Swap, 
        kron
    ))

assert(gen_decomposition((I, C, I, C, X)) == ['(I, Swap, I, I)', '(I, I, CCX)', '(I, Swap, I, I)'])