Код:
#!/usr/bin/env python
# Author: corbett@caltech.edu

import numpy as np
import unittest
import re
import random
import exceptions
import itertools
from math import sqrt,pi,e,log
import time
####
## Gates
####
class Gate(object):
	i_=np.complex(0,1)
	## One qubit gates
	# Hadamard gate
	H=1./sqrt(2)*np.matrix('1 1; 1 -1') 
	# Pauli gates
	X=np.matrix('0 1; 1 0')
	Y=np.matrix([[0, -i_],[i_, 0]])
	Z=np.matrix([[1,0],[0,-1]])

	# Defined as part of the Bell state experiment
	W=1/sqrt(2)*(X+Z)
	V=1/sqrt(2)*(-X+Z)
	
	# Other useful gates
	eye=np.eye(2,2)

	S=np.matrix([[1,0],[0,i_]])
	Sdagger=np.matrix([[1,0],[0,-i_]]) # convenience Sdagger = S.conjugate().transpose()

	T=np.matrix([[1,0],[0, e**(i_*pi/4.)]])
	Tdagger=np.matrix([[1,0],[0, e**(-i_*pi/4.)]]) # convenience Tdagger= T.conjugate().transpose()


	# TODO: for CNOT gates define programatically instead of the more manual way below
	## Two qubit gates
	# CNOT Gate (control is qubit 0, target qubit 1), this is the default CNOT gate
	CNOT2_01=np.matrix('1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0')
	# control is qubit 1 target is qubit 0 
	CNOT2_10=np.kron(H,H)*CNOT2_01*np.kron(H,H) #=np.matrix('1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0') 

	# operates on 2 out of 3 entangled qubits, control is first subscript, target second
	CNOT3_01=np.kron(CNOT2_01,eye)
	CNOT3_10=np.kron(CNOT2_10,eye)
	CNOT3_12=np.kron(eye,CNOT2_01)
	CNOT3_21=np.kron(eye,CNOT2_10)
	CNOT3_02=np.matrix('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 0 1; 0 0 0 0 0 0 1 0')
	CNOT3_20=np.matrix('1 0 0 0 0 0 0 0; 0 0 0 0 0 1 0 0; 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 1; 0 0 0 0 1 0 0 0; 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 1 0; 0 0 0 1 0 0 0 0')

	# operates on 2 out of 4 entangled qubits, control is first subscript, target second
	CNOT4_01=np.kron(CNOT3_01,eye)
	CNOT4_10=np.kron(CNOT3_10,eye)
	CNOT4_12=np.kron(CNOT3_12,eye)
	CNOT4_21=np.kron(CNOT3_21,eye)
	CNOT4_13=np.kron(eye,CNOT3_02)
	CNOT4_31=np.kron(eye,CNOT3_20)
	CNOT4_02=np.kron(CNOT3_02,eye)
	CNOT4_20=np.kron(CNOT3_20,eye)
	CNOT4_23=np.kron(eye,CNOT3_12)
	CNOT4_32=np.kron(eye,CNOT3_21)
	CNOT4_03=np.eye(16,16)
	CNOT4_03[np.array([8,9])]=CNOT4_03[np.array([9,8])]
	CNOT4_03[np.array([10,11])]=CNOT4_03[np.array([11,10])]
	CNOT4_03[np.array([12,13])]=CNOT4_03[np.array([13,12])]
	CNOT4_03[np.array([14,15])]=CNOT4_03[np.array([15,14])]
	CNOT4_30=np.eye(16,16)
	CNOT4_30[np.array([1,9])]=CNOT4_30[np.array([9,1])]
	CNOT4_30[np.array([3,11])]=CNOT4_30[np.array([11,3])]
	CNOT4_30[np.array([5,13])]=CNOT4_30[np.array([13,5])]
	CNOT4_30[np.array([7,15])]=CNOT4_30[np.array([15,7])]

	# operates on 2 out of 5 entangled qubits, control is first subscript, target second
	CNOT5_01=np.kron(CNOT4_01,eye)
	CNOT5_10=np.kron(CNOT4_10,eye)
	CNOT5_02=np.kron(CNOT4_02,eye)
	CNOT5_20=np.kron(CNOT4_20,eye)
	CNOT5_03=np.kron(CNOT4_03,eye)
	CNOT5_30=np.kron(CNOT4_30,eye)
	CNOT5_12=np.kron(CNOT4_12,eye)
	CNOT5_21=np.kron(CNOT4_21,eye)
	CNOT5_13=np.kron(CNOT4_13,eye)
	CNOT5_31=np.kron(CNOT4_31,eye)
	CNOT5_14=np.kron(eye,CNOT4_03)
	CNOT5_41=np.kron(eye,CNOT4_30)
	CNOT5_23=np.kron(CNOT4_23,eye)
	CNOT5_32=np.kron(CNOT4_32,eye)
	CNOT5_24=np.kron(eye,CNOT4_13)
	CNOT5_42=np.kron(eye,CNOT4_31)
	CNOT5_34=np.kron(eye,CNOT4_23)
	CNOT5_43=np.kron(eye,CNOT4_32)
	CNOT5_04=np.eye(32,32)
	CNOT5_04[np.array([16,17])]=CNOT5_04[np.array([17,16])]
	CNOT5_04[np.array([18,19])]=CNOT5_04[np.array([19,18])]
	CNOT5_04[np.array([20,21])]=CNOT5_04[np.array([21,20])]
	CNOT5_04[np.array([22,23])]=CNOT5_04[np.array([23,22])]
	CNOT5_04[np.array([24,25])]=CNOT5_04[np.array([25,24])]
	CNOT5_04[np.array([26,27])]=CNOT5_04[np.array([27,26])]
	CNOT5_04[np.array([28,29])]=CNOT5_04[np.array([29,28])]
	CNOT5_04[np.array([30,31])]=CNOT5_04[np.array([31,30])]
	CNOT5_40=np.eye(32,32)
	CNOT5_40[np.array([1,17])]=CNOT5_40[np.array([17,1])]
	CNOT5_40[np.array([3,19])]=CNOT5_40[np.array([19,3])]
	CNOT5_40[np.array([5,21])]=CNOT5_40[np.array([21,5])]
	CNOT5_40[np.array([7,23])]=CNOT5_40[np.array([23,7])]
	CNOT5_40[np.array([9,25])]=CNOT5_40[np.array([25,9])]
	CNOT5_40[np.array([11,27])]=CNOT5_40[np.array([27,11])]
	CNOT5_40[np.array([13,29])]=CNOT5_40[np.array([29,13])]
	CNOT5_40[np.array([15,31])]=CNOT5_40[np.array([31,15])]

####
## States
####
class State(object):
	i_=np.complex(0,1)
	## One qubit states (basis)
	# standard basis (z)
	zero_state=np.matrix('1; 0')
	one_state=np.matrix('0; 1')
	# diagonal basis (x)
	plus_state=1/sqrt(2)*np.matrix('1; 1')
	minus_state=1/sqrt(2)*np.matrix('1; -1')
	# circular basis (y)
	plusi_state=1/sqrt(2)*np.matrix([[1],[i_]])    # also known as clockwise arrow state
	minusi_state=1/sqrt(2)*np.matrix([[1],[-i_]])  # also known as counterclockwise arrow state

	# 2-qubit states
	bell_state=1/sqrt(2)*np.matrix('1; 0; 0; 1')
	@staticmethod
	def change_to_x_basis(state):
		return Gate.H*state

	@staticmethod
	def change_to_y_basis(state):
		return Gate.H*Gate.Sdagger*state

	@staticmethod
	def change_to_w_basis(state):
		# W=1/sqrt(2)*(X+Z)
		return Gate.H*Gate.T*Gate.H*Gate.S*state

	@staticmethod
	def change_to_v_basis(state):
		# V=1/sqrt(2)*(-X+Z)
		return Gate.H*Gate.Tdagger*Gate.H*Gate.S*state

	@staticmethod 
	def is_fully_separable(qubit_state):
		try:
			separated_state=State.separate_state(qubit_state)
			for state in separated_state:
				State.string_from_state(state)
			return True
		except StateNotSeparableException, e:
			return False

	@staticmethod
	def get_first_qubit(qubit_state):
		return State.separate_state(qubit_state)[0]

	@staticmethod
	def get_second_qubit(qubit_state):
		return State.separate_state(qubit_state)[1]

	@staticmethod
	def get_third_qubit(qubit_state):
		return State.separate_state(qubit_state)[2]

	@staticmethod
	def get_fourth_qubit(qubit_state):
		return State.separate_state(qubit_state)[3]

	@staticmethod
	def get_fifth_qubit(qubit_state):
		return State.separate_state(qubit_state)[4]

	@staticmethod 
	def all_state_strings(n_qubits):
		return [''.join(map(str,state_desc)) for state_desc in itertools.product([0, 1], repeat=n_qubits)]

	@staticmethod
	def state_from_string(qubit_state_string):
		if not all(x in '01' for x in qubit_state_string):
			raise Exception("Description must be a string in binary")
		state=None
		for qubit in qubit_state_string:
			if qubit=='0':
				new_contrib=State.zero_state
			elif qubit=='1':
				new_contrib=State.one_state
			if state==None:
				state=new_contrib
			else:
				state=np.kron(state,new_contrib)
		return state

	@staticmethod
	def string_from_state(qubit_state):
		separated=State.separate_state(qubit_state)
		desc=''
		for state in separated:
			if np.allclose(state,State.zero_state):
				desc+='0'
			elif np.allclose(state,State.one_state):
				desc+='1'
			else: 
				raise StateNotSeparableException("State is not separable")
		return desc

	@staticmethod
	def separate_state(qubit_state):
		"""This only works if the state is fully separable at present

		Throws exception if not a separable state"""
		n_entangled=QuantumRegister.num_qubits(qubit_state)
		if list(qubit_state.flat).count(1)==1:
			separated_state=[]
			idx_state=list(qubit_state.flat).index(1)
			add_factor=0
			size=qubit_state.shape[0]
			while(len(separated_state)<n_entangled):
				size=size/2
				if idx_state<(add_factor+size):
					separated_state+=[State.zero_state]
					add_factor+=0
				else:
					separated_state+=[State.one_state]
					add_factor+=size
			return separated_state
		else:
			# Try a few naive separations before giving up
			cardinal_states=[State.zero_state,State.one_state,State.plus_state,State.minus_state,State.plusi_state,State.minusi_state]
			for separated_state in itertools.product(cardinal_states, repeat=n_entangled):
				candidate_state=reduce(lambda x,y:np.kron(x,y),separated_state)
				if np.allclose(candidate_state,qubit_state):
					return separated_state
			# TODO: more general separation methods
			raise StateNotSeparableException("TODO: Entangled qubits not represented yet in quantum computer implementation. Can currently do manual calculations; see TestBellState for Examples")

	@staticmethod
	def measure(state):
		"""finally some probabilities, whee. To properly use, set the qubit you measure to the result of this function
			to collapse it. state=measure(state). Currently supports only up to three entangled qubits """
		state_z=state
		n_qubits=QuantumRegister.num_qubits(state)
		probs=Probability.get_probabilities(state_z)
		rand=random.random()
		for idx,state_desc in enumerate(State.all_state_strings(n_qubits)):
			if rand < sum(probs[0:(idx+1)]):
				return State.state_from_string(state_desc)

	@staticmethod
	def get_bloch(state):
		return np.array((Probability.expectation_x(state),Probability.expectation_y(state),Probability.expectation_z(state)))

	@staticmethod
	def pretty_print_gate_action(gate,n_qubits):
		for s in list(itertools.product([0,1], repeat=n_qubits)):
			sname=('%d'*n_qubits)%s
			state=State.state_from_string(sname)
			print sname,'->',State.string_from_state(gate*state)

class StateNotSeparableException(exceptions.Exception):
	def __init__(self,args=None):
		self.args=args

class Probability(object):
	@staticmethod
	def get_probability(coeff):
		return (coeff*coeff.conjugate()).real

	@staticmethod
	def get_probabilities(state):
		return [Probability.get_probability(x) for x in state.flat]

	@staticmethod
	def get_correlated_expectation(state):
		probs=Probability.get_probabilities(state)
		return probs[0]+probs[3]-probs[1]-probs[2]
5-qubit симулация на квантов компютър