Compiler for a quantum computer

Abstract

A quantum computer is an array of quantum bits (qubits) together with some hardware for manipulating these qubits. quantum Bayesian (qb) nets are a method of modeling quantum systems graphically in terms of network diagrams. This invention comprises a classical computer running a computer program that expresses the information contained in a qb net as a sequence of elementary operations (SEO). One can then run these sequences on a quantum computer. We show how to reduce a qb net into a SEO by a two step process. First, express the information contained in the qb net as a sequence of unitary operators. Second, express each of those unitary operators as a SEO. An appendix to this document contains the C++ source code of a computer program called “Qubiter1.0”, which is a preferred embodiment of the invention.

Description

REFERENCE TO A MICROFICHE APPENDIX

The present application includes a microfiche appendix comprising the C++ source code of a fully functional computer program called Qubiter1.0. Qubiter1.0 is a possible embodiment of the software of the present invention. The microfiche appendix comprises 1 microfiche with a total of 52 frames. The first frame is a test pattern for focusing. The second frame, called Appendix A, is a list of files contained in a CodeWarrior™ project for Qubiter1.0. Subsequent frames are labelled Appendix B, and comprise source code contained within said files.

TECHNICAL FIELD

The invention relates to an array of quantum bits known as a quantum computer. More specifically, it relates to the generation of the instruction sequences that are used to manipulate such an array.

BACKGROUND OF THE INVENTION

This invention deals with quantum computers. A quantum computer is an array of quantum bits (qubits) together with some hardware for manipulating these qubits. Quantum computers with only a few bits have already been built. For a review of quantum computers, see DiV95: D. P. DiVincenzo, Science 270, 255 (1995). See also Ste97:A. M. Steane, Los Alamos eprint

The function A_j with values A_j[x_j|(x.)_Sj] determines a matrix that we will call the node matrix of node {circumflex over (x)}_j, and denote by Q_j. Thus, x_j is the matrix's row index and (x.)_Sj is its column index.

For example, FIG. 1 gives the four node matrices Q₁, Q₂, Q₃, Q₄ associated with the four nodes of the graph shown there.

One can translate a QB net into a SEO by performing 4 steps: (1) Find eras. (2) Insert delta functions. (3) Find unitary extensions of era matrices. (4) Decompose each unitary matrix into a SEO. Next we will discuss these 4 steps in detail. We will illustrate our discussion by using Teleportation as an example. Teleportation was first discussed in Tel93: C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, W. K. Wootters, Physical Review Letters 70, 1895 (1993). FIG. 2 shows a QB net for Teleportation. Reference Bra96: G. Brassard, Los Alamos eprint
B₂(x₂,x₃|x₁)=A₂(x₂|x₁)A₃(x₃|x₁),   (5b)
B₃(x₅|x₂,x₄)=A₅(x₅|x₂,x₄),   (5c)
B₄(x₆|x₃,x₅)=A₆(x₆|x₃,x₅),   (5d)
and
A(x.)=B₄B₃B₂B₁.   (6)
Step 2: Insert Delta Functions

The Feynman Integral FI for a QB net is defined by $\begin{array}{cc}\mathrm{FI}\left[{\left(x.\right)}_{Z_{\mathrm{ext}}}\right]=\sum _{\left(x,\right)z_{\mathrm{int}}}A\left(x.\right).& \left(7\right)\end{array}$
Note that we are summing over all stories x. that have (x.)_Zext as their external state. We want to express the right side of Eq. (7) as a product of matrices.

Consider how to do this for Teleportation. In that case one has $\begin{array}{cc}\mathrm{FI}\left(x_6\right)=\sum _{x_1,x_2,{\mathrm{\dots x}}_5}{B}_4{B}_3{B}_2{B}_1,& \left(8\right)\end{array}$
where the B_a are given by Eqs. (5). The right side of Eq. (8) is not ready to be expressed as a product of matrices because the column indices of B_a+1 and the row indices of B_a are not the same for all aεZ_1,|τ|−1. Furthermore, the variable x₃ occurs in B₄ and B₂ but not in B₃. Likewise, the variable x₄ occurs in B₃ and B₁ but not in B₂. Suppose we define B_a for aεZ_1,|τ| by $\begin{array}{cc}{\stackrel{\_}{B}}_1\left(x_1^1,x_4^1\right)=B_1\left(x_1^1,x_4^1\right),& \left(9a\right)\end{array}$ $\begin{array}{cc}{\stackrel{\_}{B}}_2\left(x_2^2,x_3^2,x_4^2❘x_1^1,x_4^1\right)=B_2\left(x_2^2,x_3^2❘x_1^1\right)\delta \left(x_4^2,x_4^1\right),& \left(9b\right)\\ {\stackrel{\_}{B}}_3\left(x_3^3,x_5^3❘x_2^2,x_3^2,x_4^2\right)=B_3\left(x_5^3❘x_2^2,x_4^2\right)\delta \left(x_3^3,x_3^2\right),& \left(9c\right)\\ {\stackrel{\_}{B}}_4\left(x_6❘x_3^3,x_5^3\right)=B_4\left(x_6❘x_3^3,x_5^3\right).\mathrm{Then}& \left(9d\right)\\ \mathrm{FI}\left(x_6\right)=\sum _{\mathrm{interm}}{\stackrel{\_}{B}}_4{\stackrel{\_}{B}}_3{\stackrel{\_}{B}}_2{\stackrel{\_}{B}}_1,& \left(10\right)\end{array}$
where we sum over all intermediate indices; i.e., all x_j^a except x₆. Contrary to Eq. (8), the right side of Eq. (10) can be expressed immediately as a product of matrices since now B_a+1 column indices and B_a row indices are the same. The purpose of inserting a delta function of x₃ into B₃ is to allow the system to “remember” the value of X₃ between non-consecutive eras T₄ and T₂. Inserting a delta function of x₄ into B₂ serves a similar purpose.

In the Teleportation net of FIG. 2, the last era contains all the external nodes. However, for some QB nets like the one in FIG. 5, this is not the case. For the net of FIG. 5.
B₁(x₁)=A₁(x₂),   (11a)
B₂(x₂,x₃|x₁)=A₃(x₃|x₁)A₂(x₂|x₁),   (11b)
B₃(x₄|x₃)=A₄(x₄|x₃),   (11c)
B₄(x₅|x₄)=A₅(x₅|x₄).   (11d)
Even though node {circumflex over (x)}₂ is external, the variable x₂ does not appear as a row index in B₄. Suppose we set $\begin{array}{cc}{\stackrel{\_}{B}}_1\left(x_1^1\right)=B_1\left(x_1^1\right),& \left(12a\right)\\ {\stackrel{\_}{B}}_2\left(x_2^2,x_3^2❘x_1^1\right)=B_2\left(x_2^2,x_3^2❘x_1^1\right),& \left(12b\right)\\ {\stackrel{\_}{B}}_3\left(x_2^3,x_4^3❘x_2^2,x_3^2,\right)=B_3\left(x_4^3❘x_3^2\right)\delta \left(x_2^3,x_2^2\right),& \left(12c\right)\\ {\stackrel{\_}{B}}_4\left(x_2,x_5❘x_2^3,x_4^3\right)=B_4\left(x_5❘x_4^3\right)\delta \left(x_2,x_2^3\right).\mathrm{Then}& \left(12d\right)\\ \mathrm{FI}\left(x_2,x_5\right)=\sum _{\mathrm{interm}}{\stackrel{\_}{B}}_4{\stackrel{\_}{B}}_3{\stackrel{\_}{B}}_2{\stackrel{\_}{B}}_1,& \left(13\right)\end{array}$
where we sum over all intermediate indices; i.e., all x_j^a except x₂ and x₅. Contrary to B₄, the rows of B₄ are labelled by the indices of both external nodes {circumflex over (x)}₂ and {circumflex over (x)}₅.

This technique of inserting delta functions can be generalized as follows to deal with arbitrary QB nets. For jεZ_1,N, let a_min(j) be the smallest aεZ_1,|τ| such that x_j appears in B_a. Hence, a_min(j) is the first era in which x_j appears. If {circumflex over (x)}_j is an internal node, let a_max(j) be the largest a such that x_j appears in B_a (i.e., the last era in which x_j appears). If {circumflex over (x)}_j is an external node, let a_max(j)=|τ|+1. For aεZ_1,|τ|+1. For aεZ_1,|τ|, let
Δ_a={jεZ_1,N|a_min(j)<a<a_max(j)},   (14)
$\begin{array}{cc}{\stackrel{\_}{B}}_a=B_a\left[{\left(x\stackrel{a}{.}\right)}_{T_a}❘{\left(x\stackrel{a-1}{.}\right)}_{{\Gamma }_a}\right]\prod _{j\in \mathrm{\Delta a}}\delta \left(x_j^a,x_j^{a-1}\right).& \left(15\right)\end{array}$
In Eq. (15), x_j^|T| should be identified with x_j and x_j⁰ with no variable at all. Equation (7) for FI can be written in terms of the B_a functions: $\begin{array}{cc}\mathrm{FI}\left[{\left(x.\right)}_{Z_{\mathrm{ext}}}\right]=\sum _{\mathrm{interm}}{\stackrel{\_}{B}}_{T}\dots {\stackrel{\_}{B}}_2{\stackrel{\_}{B}}_1,& \left(16\right)\end{array}$
where the sum is over all intermediate indices (i.e., all x_j^a for which a≠|τ|). For all a, define matrix M_a so that the x,y entry of M_a is B_a(x|y). Define M to be a column vector whose components are the values of FI for each external state. Then Eq. (16) can be expressed as:
M=M_|τ| . . . M₂M₁.   (17)
The rows of the column vector M are labelled by the possible values of (x.)_Zext. The rows of the column vector M₁ are labelled by the possible values of (x.)_T1, where T₁ is the set of root nodes.
Step 3: Find Unitary Extensions of Era Matrices

So far, we have succeeded in expressing FI as a product of matrices M_a, but these matrices are not necessarily unitary. In this step, we will show how to extend each M_a matrix (by adding rows and columns) into a unitary matrix U_a.

By combining adjacent Ma's, one can produce a new, smaller set of matrices M_a. Suppose the union of two consecutive eras is also defined to be an era. Then combining adjacent M_a's is equivalent to combining consecutive eras to produce a new, smaller set of eras. We define a breakpoint as any position aεZ_1,|τ|−1 between two adjacent matrices M_a+1 and M_a. Combining two adjacent M_a's eliminates a breakpoint. Breakpoints are only necessary at positions where internal measurements are made. For example, in Teleportation experiments, one measures node {circumflex over (x)}₅, which is in era T₃. Hence, a breakpoint between M₄ and M₃ is necessary. If that is the only internal measurement to be made, all other breakpoints can be dispensed with. Then we will have M=M₂′M₁′ where M₂′=M₄, M₁′=M₃M₂M₁. If no internal measurements are made, then we can combine all matrices M_a into a single one, and eliminate all breakpoints.

We will henceforth assume that for all aεZ_1,|τ|, the columns of M_a are orthonormal. If for some a₀εZ_1,|τ|, M_a0 does not satisfy this condition, it may be possible to “repair” M_a0, so that it does. First: If a row β of M_a0−1 is zero, then eliminate the column β of M_a0, and the row β of M_a0−1. Next: If a row β of the column vector M_a0−1 . . . M₂M₁ is zero, then flag the column β of M_a0. The flagged columns of M_a0 can be changed without affecting the value of M. If the non-flagged columns of M_a0 are orthonormal, and the number of columns in M_a0 does not exceed the number of rows, then the Gram Schmidt method, to be discussed later, can be used to replace the flagged columns by new columns such that all the columns of the new matrix M_a0 are orthonormal. If it is not possible to repair M_a0 in any of the above ways (or in some other way that might become clear once we program this), one can always remove the breakpoint between M_a0+1 and M_a0.

We will call d_a the number of rows of matrix M_a and d_a−1 its number of columns. We define D and N_S by
D=max{d_a|1≦a≦|τ|}  (18)
N_S=min{2ⁱ|iεZ_1,∞,D≦2ⁱ}.   (19)
Let d_a=N_S−d_a for all a. For each a≠1. we define U_a to be the matrix that one obtains by extending M_a as follows. We append an d_a×d_a−1 block of zeros beneath M_a and an N_S× d_a−1 block of gray entries to the right of M_a. By gray entries we mean entries whose value is yet to be specified. When a=1, M₁ be can extended in two ways. One can append a column of d₁ zeros beneath it and call the resulting N_S dimensional column vector v₁. Alternatively, one can append a column of d₁ zeros beneath M₁ and an N_S×(N_S−1 block of gray entries to the right of M₁, and call the resulting N_S×N_S matrix U₁. In this second case, one must also insert e₁ to the right of U₁. By e₁ we mean the N_S dimensional column vector whose first entry equals one and all others equal zero. Which extension of M₁ is used, whether the one that requires e₁ or the one that doesn't, should be left as a choice of the user. Henceforth, for the sake of definiteness, we will assume that the user has chosen the extension without the e₁. The other case can be treated similarly. Equation (17) then becomes
v=U_|τ| . . . U₃U₂v₁,   (20)
where v is just the column vector M with d_|τ| zeros attached to the end.

To determine suitable values for the gray entries of the U_a matrices, one can use the Gram-Schmidt (G.S.) method. (See Nob88: B. Noble and J. W. Daniels, Applied Linear Algebra, Third Edition (Prentice Hall, 1988)). This method takes as input an ordered set S=(v₁, v₂, . . . , v_N) of vectors, not necessarily independent ones. It yields as output another ordered set of vectors S′=(u₁, u₂, . . . , u_N), such that S′ spans the same vector space as S. Some vectors in S′ may be zero. Those vectors of S′ which aren't zero will be orthonormal. For rεZ_1,N, if the first r vectors of S are already orthonormal, then the first r vectors of S′ will be the same as the first r vectors of S. Let e_j for jεZ_1,NS be the j'th standard unit vector (i.e., the vector whose j'th entry is one and all other entries are zero). For each aεZ_1,|τ|, to determine the gray entries of U_a one can use the G.S. method on the set S consisting of the non-gray columns of U_a together with the vectors e₁, e₂, . . . e_Ns.

Step 4: Decompose Each Unitary Matrix into a SEO

In this section we present a CS method for decomposing an arbitrary unitary matrix into a SEO. By following the previous 3 steps, one can reduce a QB net to a product of unitary operators U_a. By applying the CS method of this section to each of the matrices U_a, one can reduce the QB net to a SEO.

We will use the symbol N_B for the number (≦1) of bits and N_S=2^NB for the number of states with N_B bits. We define Bool={0,1}. We will use lower case Latin letters a,b,c . . . εBool to represent bit values and lower case Greek letters α,β,δ, . . . εZ_0,Nb−1 to represent bit positions. A vector such as {right arrow over (a)}=a_Nb−1 . . . a₂a₁a₀ will represent a string of bit values, a_μ being the value of the μ'th bit for μεZ_0,Nb−1. A bit string {right arrow over (a)} has a decimal representation $\begin{array}{cc}d\left(\ddot{a}\right)=\sum _{\mu =0}^{N_M-1}{2}^{\mu }{a}_{\mu }.& \left(21\right)\end{array}$
For βεZ_0,NB−1, we will use {right arrow over (u)}(β) to denote the β'th standard unit vector, i.e, the vector with bit value of 1 at bit position β and bit value of zero at all other bit positions.

I_r will represent the r dimensional unit matrix. Suppose βεZ_0,NB−1 and M is any 2×2 matrix. We define M(β) by
M(β)=I₂(×) . . . (×)I₂(×)M(×)I₂(×) . . . (×)I₂,   (22)
where the matrix M on the right side is located at bit position β in the tensor product of N_B2×2 matrices. The numbers that label bit positions in the tensor product increase from right to left (←), and the rightmost bit is taken to be at position 0.

For any two square matrices A and B of the same dimension, we define the (·) product by A(·)B=ABA^†, where A^† is the Hermitian conjugate of A.

{right arrow over (σ)}=(σ_x,σ_y,σ_z) will represent the vector of Pauli matrices, where $\begin{array}{cc}{\sigma }_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_y=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).& \left(23\right)\end{array}$

The Sylvester-Hadamard matrices H_r are 2^r×2^r matrices whose entry at row {right arrow over (a)} and column {right arrow over (b)} is given by
(H_r)_{{right arrow over (a)}, {right arrow over (b)}}=(−1)^{{right arrow over (a)}·{right arrow over (b)}},   (24)
where $\begin{array}{cc}\stackrel{\rightharpoonup }{a}·\stackrel{\rightharpoonup }{b}=\sum _{\mu =0}^{r-1}{a}_{\mu }{b}_{\mu }.& \left(25\right)\end{array}$

The qubit's basis states |0< and |1> will be represented by $\begin{array}{cc}0\rangle =\left[\begin{array}{c}1\\ 0\end{array}\right],1\rangle =\left[\begin{array}{c}0\\ 1\end{array}\right].& \left(26\right)\end{array}$

The number operator n of the qubit is defined by $\begin{array}{cc}n=\left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right]=\frac{1-{\sigma }_x}{2}.& \left(27\right)\end{array}$
Note that
n|0>=0,n|1>=|1>.   (28)

We will often use n as shorthand for $\begin{array}{cc}\stackrel{\_}{n}=1-n=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]=\frac{1+{\sigma }_x}{2}.& \left(29\right)\end{array}$
We define P₀ and P₁ by $\begin{array}{cc}{P}_0=\stackrel{\_}{n}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right],P_1=n=\left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right].& \left(30\right)\end{array}$
For βεZ_0,NB−1, we define P₀(β),P₁(β), n(β) and n(β) by means of Eq. (22). For {right arrow over (a)}εBool^BB,let
P_{{right arrow over (a)}}=P_aNB−1(×) . . . (×)P_a2(×)P_a1(×)P_a0.   (31)

As mentioned earlier, we utilize a mathematical technique called the CS Decomposition. In this name, the letters C and S stand for “cosine” and “sine”. Next we will state the special case of the CS Decomposition Theorem that arises in a preferred embodiment of the invention.

Suppose that U is an N×N unitary matrix, where N is an even number. Then the CS Decomposition Theorem states that one can always express U in the form $\begin{array}{cc}U=\left[\begin{array}{cc}{L}_0& 0\\ 0& L_1\end{array}\right]D\left[\begin{array}{cc}{R}_0& 0\\ 0& R_1\end{array}\right],& \left(32\right)\end{array}$
where L₀,L₁,R₀,R₁ are N/2×N/2 unitary matrices and $\begin{array}{cc}D\left[\begin{array}{cc}{D}_{00}& D_{01}\\ D_{10}& D_{11}\end{array}\right],& \left(33a\right)\\ D_{00}=D_{11}=\mathrm{diag}\left(C_1,C_2,\dots ,C_{\frac{N}{2}}\right),& \left(33b\right)\\ D_{01}=\mathrm{diag}\left(S_1,S_2,\dots ,S_{\frac{N}{2}}\right),& \left(33c\right)\\ D_{10}=-D_{01}.& \left(33d\right)\end{array}$
For $i\in \left\{1,2,\dots ,\frac{N}{2}\right\},$
C_i and S_i are real numbers that satisfy
C_i²+S_i²=1.   (33e)
Henceforth, we will use the term D matrix to refer to any matrix that satisfies Eqs. (33). If one partitions U into four blocks U_ij of size $\frac{N}{2}\times \frac{N}{2},$
then
U_ij=L_iD_ijR_j,   (34)
for i,jε{0,1}. Thus, D_ij gives the singular values of U_ij.

More general versions of the CS Decomposition Theorem allow for the possibility that we partition U into 4 blocks of unequal size.

Note that if U were a general (not necessarily unitary) matrix, then the four blocks U_ij would be unrelated. Then to find the singular values of the four blocks U_ij would require eight unitary matrices (two for each block), instead of the four L_i,R_j. This double use of the L_i,R_j is a key property of the CS decomposition.

Consider FIG. 6. We start at 61 with a unitary matrix U. Without loss of generality, we can assume that the dimension of U is 2^NB for some N_B≧1. (If initially U's dimension is not a power of 2, we replace it by a direct sum U(+)I_r whose dimension is a power of two.) We apply the CS Decomposition method to U. This yields node 62 comprising matrix D(0,U) of singular values, two unitary matrices L(0,U) and L(1,U) on the left and two unitary matrices R(0,U) and R(1,U) on the right. Then we apply the CS Decomposition method to each of the 4 matrices L(0.U),L(1.U),R(0,U) and R(1,U) and obtain nodes 63 and 64. Then we apply the CS Decomposition method to each of the 16 R and L matrices in nodes 63 and 64. And so on. This process ends when the current row of nodes in the pyramid of FIG. 6 has L's and R's that are 1×1 dimensional, i.e., just complex numbers.

Call a central matrix either (1) a single D matrix, or (2) a direct sum D₁(+)D₂(+) . . . (+)D_r of D matrices, or (3) a diagonal unitary matrix. From FIG. 6 it is clear that the initial matrix U can be expressed as a product of central matrices, with each node of the tree providing one of the central matrices in the product. Next, we show how to decompose each of the 3 possible kinds of central matrices into a SEO.

Case 1: Central Matrix is a Single D Matrix

Consider how to decompose a central matrix when it is a single D matrix. Before dealing with arbitrary N_B consider N_B=3. Then the central matrix D can be expressed as: $\begin{array}{cc}D=\sum _{a,b\in \mathrm{Bool}}\exp \left({i\varphi }_{\mathrm{ab}}{\sigma }_y\right)(\times )P_a(\times )P_b.& \left(35\right)\end{array}$
Suppose {right arrow over (φ)} (ditto, {right arrow over (θ)}) is a column vector whose components are the numbers φ_{{right arrow over (a)}} (ditto, θ_{{right arrow over (a)}}) arranged in order of increasing {right arrow over (a)}. We define new angles {right arrow over (θ)} in terms of the angles {right arrow over (φ)} by $\begin{array}{cc}\stackrel{\rightharpoonup }{\theta }=\frac{1}{4}{H}_2\stackrel{\rightharpoonup }{\varphi }.& \left(36\right)\end{array}$
Then one can show that
D=A₀₀A₀₁A₁₀A₁₁,   (37)
where
A₀₀=exp(iθ₀₀σ_y)(×)I₂(×)I₂,   (38a)
A₀₁=σ_x(2)ⁿ⁽⁰⁾(·)[exp(iθ₀₁σ_y)(×)I₂(×)I₂],   (38b)
A₁₀=σ_x(2)ⁿ⁽¹⁾(·)[exp(iθ₁₀σ_y)(×)I₂(×)I₂],   (38c)
A₁₁=[σ_x(2)ⁿ⁽¹⁾σ_x(2)ⁿ⁽⁰⁾](·)[exp(iθ₁₁σ_y)(×) I₂(×)I₂].   (38d)
Eqs. (37)-(38) achieve our goal of decomposing D into a SEO. Now consider an arbitrary N_B. D can be written as $\begin{array}{cc}D=\sum _{\stackrel{\rightharpoonup }{a}\in {\mathrm{Bool}}^{N_B-1}}\exp \left({i\varphi }_{\stackrel{\rightharpoonup }{a}}{\sigma }_y\right)(\times )P_{\stackrel{\rightharpoonup }{a}},& \left(39\right)\end{array}$
where the φ_{{right arrow over (a)}} are real numbers. We define new angles {right arrow over (θ)} in terms of the angles {right arrow over (φ)} by $\begin{array}{cc}\stackrel{\rightharpoonup }{\theta }=\frac{1}{2^{N_B-1}}{H}_{N_B-1}\stackrel{\rightharpoonup }{\varphi }.\mathrm{Then}\mathrm{one}\mathrm{can}\mathrm{show}\mathrm{that}& \left(40\right)\\ D=\prod _{\stackrel{\rightharpoonup }{b}\in {\mathrm{Bool}}^{N_B-1}}{A}_{\stackrel{\rightharpoonup }{b}},& \left(41\right)\end{array}$
where the operators A_{{right arrow over (b)}} on the right side commute and will be defined next. For any {right arrow over (b)}εbool^ND−1 we can write $\begin{array}{cc}\stackrel{\rightharpoonup }{b}=\sum _{j=0}^{r-1}\stackrel{\rightharpoonup }{u}\left({\beta }_j\right),& \left(42\right)\end{array}$
where
N_B−2≧β_r−1>. . . β₁>B₀≧o.   (43)
In other words, {right arrow over (b)} has bit value of 1 at bit positions β_j. At all other bit positions, {right arrow over (b)} has bit value of 0. r is the number of bits in {right arrow over (b)} whose value is 1. When {right arrow over (b)}=0, r=0. One can show that
A_{{right arrow over (b)}}=[σ_x(N_B−1)^n(βr−1). . . σ_x(N_B−1)^n(β1)σ_x(N_N−1)^n(β0)](·)exp[iθ_{{right arrow over (b)}σy}(N_B−1]  (44)
There are other ways of decomposing A_{{right arrow over (b)}} into a SEO.

Case 2: Central Matrix is Direct Sum of D Matrices

Next, consider how to decompose a central matrix when it is a direct sum of D matrices. Consider first the case N_B=3. Let R(φ)=exp(iσ_yφ). Previously we mentioned the fact that any D matrix D can be expressed as $\begin{array}{cc}D=\sum _{a,b\in \mathrm{Bool}}R\left({\varphi }_{\mathrm{ab}}^{″}\right)(\times )P_a(\times )P_b.\mathrm{One}\mathrm{can}\mathrm{also}\mathrm{show}\mathrm{that}& \left(46\right)\\ D_0(+)D_1=\sum _{a,b\in \mathrm{Bool}}{P}_a(\times )R\left({\varphi }_{\mathrm{ab}}^{\prime }\right)(\times )P_b,& \left(47\right)\\ D_{00}(+)D_{01}(+)D_{10}(+)D_{11}=\sum _{a,b\in \mathrm{Bool}}{P}_a(\times )P_b(\times )R\left({\varphi }_{\mathrm{ab}}\right),& \left(48\right)\end{array}$
where the D_j and D_ij are D matrices. It follows that by permuting the bit positions, we can change such a direct sum of D matrices into a single D matrix. The latter can then be decomposed into a SEO by the method already discussed.

Case 3: Central Matrix is a Diagonal Unitary Matrix

Finally, consider how to decompose a central matrix when it is a diagonal unitary matrix. Before dealing with arbitrary N_B, consider N_B=2. Then the central matrix C can be expressed as
C=diag(e^iφoo, e^iφ01, e^iφ10, e^iφ11).   (49)
We define {right arrow over (φ)} by $\begin{array}{cc}\stackrel{\rightharpoonup }{\theta }=\frac{1}{4}{H}_2\stackrel{\rightharpoonup }{\varphi }.& \left(50\right)\end{array}$
One can show that
C=A₀₀A₀₁A₁₀A₁₁,   (51)
where
A₀₀=exp(iθ₀₀),   (52a)
A₀₁=I₂(×)exp(iθ₀₁σ_z), (52b)
A₁₀=exp(iθ₁₀σ_z)(×)I₂,   (52c)
A₁₁=σ_x(0)ⁿ⁽¹⁾(·)[I₂(×)exp(iθ₁₁σ_z)].   (52d)
Now consider an arbitrary N_B. Any diagonal unitary matrix C can be expressed as $\begin{array}{cc}C=\sum _{\stackrel{\rightharpoonup }{a}\in {\mathrm{Bool}}^{N_B}}\exp \left({i\varphi }_{\stackrel{\rightharpoonup }{a}}\right)P_{\stackrel{\rightharpoonup }{a}},& \left(53\right)\end{array}$
where the φ_{{right arrow over (a)}} are real numbers. We define {right arrow over (θ)} by $\begin{array}{cc}\stackrel{\rightharpoonup }{\theta }=\frac{1}{2^{N_B}}{H}_{N_B}\stackrel{\rightharpoonup }{\varphi }.\mathrm{Then}\mathrm{one}\mathrm{can}\mathrm{show}\mathrm{that}& \left(54\right)\\ C=\prod _{\stackrel{\rightharpoonup }{b}\in {\mathrm{Bool}}^{N_B}}{A}_{\stackrel{\rightharpoonup }{b}},& \left(55\right)\end{array}$
where the A_{{right arrow over (b)}} operators commute and will be defined next. For any _{{right arrow over (b)}}εBool^NB, we can write $\begin{array}{cc}\stackrel{\rightharpoonup }{b}=\sum _{j=0}^{r-1}\stackrel{\rightharpoonup }{u}\left({\beta }_j\right),& \left(56\right)\end{array}$
where
N_B−1≧β_r−1>. . . >β₁>β₀≧0.  (57)
One can show that $\begin{array}{cc}{A}_{\stackrel{\rightharpoonup }{b}}=\{\begin{array}{c}\exp \left[{i\theta }_0\right]\mathrm{if}r=0\\ \exp \left[{i\theta }_{\stackrel{\rightharpoonup }{b}}{\sigma }_z\left({\beta }_0\right)\right]\mathrm{if}r=1\\ \begin{array}{c}\left[{{\sigma }_x\left({\beta }_0\right)}^{n\left({\beta }_{r-1}\right)}{{\mathrm{\cdots \sigma }}_x\left({\beta }_0\right)}^{n\left({\beta }_2\right)}{{\sigma }_x\left({\beta }_0\right)}^{n\left({\beta }_1\right)}\right]\\ (·)\exp \left[{i\theta }_{\stackrel{\rightharpoonup }{b}}{\sigma }_z\left({\beta }_0\right)\right]\mathrm{if}r\ge 2\end{array}\end{array}& \left(58\right)\end{array}$
There are other ways of decomposing A_{{right arrow over (b)}} into a SEO.
(B) Implementation of New Method on Classical Computer

So far in Section (A), we have described a mathematical algorithm for decomposing any QB net into a SEO. Next we describe a particular implementation of the algorithm, a computer program called “Qubiter” that can be run on a classical computer.

The use of a computer is practically indispensable for obtaining useful numerical answers through the method of Section (A). In all but the simplest of cases, vast amounts of data storage and processing are necessary to obtain final numerical answers from the method of Section (A). The necessary book-keeping and number crunching are prohibitively error prone and time consuming to a human, but not to a computer.

A preferred embodiment of the invention is a classical computer that feeds data (a SEO) to a quantum computer. By a classical computer, we mean a device that makes a desired calculation using digital circuits which implement deterministic (classical, non-quantum) logic. By a quantum computer we mean a device that makes a desired calculation using an array of quantum bits (qubits). Besides their calculational circuits, classical and quantum computers may comprise input, output and memory devices. The important difference is that an array of quantum bits may be put in an entangled quantum state, whereas a digital deterministic logic circuit cannot be put in such a state (in practice, for useful periods of time). The classical computer of our preferred embodiment is a Mac computer, produced by Apple Computer Inc. of Cupertino, California. The Mac is running a program written in the computer language C++. Of course, this invention could just as easily be implemented in a language other than C++, and on a platform other than a Mac. FIG. 7 is a block diagram of a classical computer feeding data to a quantum computer. 70 represents a classical computer. It comprises 71, 72, 73. 71 represents input devices, such as a mouse or a keyboard. 72 represents the CPU, internal and external memory units. 72 does calculations and stores information. 73 represents output devices, such as a printer or a display screen. An image of a QB net (for example, FIG. 2) can be rendered on the display screen. 75 represents a quantum computer, comprising an array of quantum bits and some hardware for manipulating the state of those bits.

Software for a preferred embodiment of the invention was written using Code Warrior™. Code Warrior is a C++ Integrated Development Environment produced by MetroWerks Inc. of Austin, Texas. C++ source code for a computer program called “Qubiter1.0” is included as a Microfiche Appendix to this document. The Microfiche Appendix has two parts: Appendix A and Appendix B.

Appendix A is a listing of the names of all the files in the “CodeWarrior Project” for Qubiter. Appendix B is a listing of Qubiter source code. Apart from libraries provided with CodeWarrior, Qubiter requires parts of the C library called Clapack. This library is freeware and can be downloaded via the Internet from the website “www.netlib.org”.

Qubiter uses files that list the entries of a unitary matrix U. FIG. 8 shows an example of such files. Qubiter ignores the lines labelled 81 and 83. The line labelled 82 tells how many bits are necessary to label U's columns (or rows, since U must be square). In the example of FIG. 8, according to 82, two bits are necessary to label U's column, so one expects U to be a 4×4 matrix with 16 elements. 84 labels all the lines after the line labelled 83. Each line in 84 represents a complex number z, the real part of z is listed first, then some white space to the right of it, and then the imaginary part of z to the right of that. 84 contains the entries of U, arranged so that the entries in the first column of U are listed first, then those of the second column of U, etc.

Qubiter also uses files that list SEOs. FIG. 9 is an example of such a file. The line labelled 91 tells how many bits the gates listed in the file are to operate on, (i.e., what we called N_B in Section (A)). 92 labels all the other lines in the file. Each line in 92 represents a gate. There are four types of gates, and they are specified as follows:

(a) PHAS ang

where ang is a real number. This signifies a phase factor $\exp \left(i\left(\mathrm{ang}\right)\frac{\pi }{180}\right).$
(b) CNOT α a char β

• • where α,βεZ_0,NB−1 and char ε{T,F}. T and F stand for true and false. If char is the character T, this signifies σ_x(β)^n(α). Read it as “c-not: if α is true, then flip β.” If char is the character F, this signifies σ_x(β) ^n(α)n. Read it as “c-not: if α is false, then flip β.”
    (c) ROTY α ang
  • where αεZ_0,NB−1 and ang is a real number. This signifies the rotation of qubit α about the Y axis by an angle ang in degrees. In other words, $\exp \left({i\sigma }_y\left(\alpha \right)\mathrm{ang}\frac{\pi }{180}\right).$
    (d) ROTZ α ang
  • This is the same as (c) except that the rotation is about the Z axis instead of the Y one.

The matrix given in FIG. 8 can be decomposed into the SEO given in FIG. 9. Such decompositions are not unique. After doing the trivial optimization of removing all factors A_{{right arrow over (b)}} for which the rotation angle is zero, the 33 operations in FIG. 9 reduce to 25 operations in FIG. 10.

Qubiter starts by looking for a parameter file entitled “qbtr-params.in”. FIG. 11 is an example of such a file. Lines labelled 111, 113, 115 and 117 are ignored by the program. Lines labelled 112, 114, 116, 118 are not.

The user should enter into line 112 the name of a matrix. In FIG. 11, we have used “DiscFourier2bits”.

If the user enters “yes” into line 114 as an answer to “Do compilation?”, then Qubiter will look for a file named “DiscFourier2bits.in”. In other words, it will look for a file whose name is the string in line 112 plus the suffix “.in”. Qubiter expects to find in this file the entries of the unitary matrix U_initial to be decomposed. The file should be of the form exemplified by FIG. 8.

If the user enters “yes” into line 116 as an answer to “Do decompilation?”, then Qubiter will look for a file named “DiscFourier2bits-cmnd.out”. In other words, it will look for a file whose name is the string in line 112 plus the suffix “-cmnd.out”. Qubiter expects to find in this file a SEO. The file should be of the form exemplified by FIGS. 9 and 10.

If the user enters “yes” into line 118 as an answer to “Do zero angle optimization?”, then Qubiter will produce a file of the type exemplified by FIG. 10 instead of the type exemplified by FIG. 9. Thus, it will omit those gates arising from a rotation by a zero angle.

Qubiter has 2 main modes of operation.

The first mode of operation is when the user enters “yes” in line 114 as an answer to “Do compilation?”, and “yes” in line 116 as an answer to “Do decompilation?” In this mode, the user must provide 2 input files entitled “qbtr-params.in” and “mat.in”, where “mat” is the string in line 112. Qubiter will output a file called “mat-cmnd.out”. Then it will use mat-cmnd.out as input, multiply the SEO listed in this file, arrive at a unitary matrix U_final, and output a file called “mat-check.out” which lists the entries of U_final. If everything goes well, the matrix U_initial specified by file “mat.in” and the matrix U_final specified by file “mat-check.out” will be the same matrix, within machine precision.

The second mode of operation is when the user enters “no” in line 114 as an answer to “Do compilation?”, and “yes” in line 116 as an answer to “Do decompilation?” In this mode, the user must provide 2 input files entitled “qbtr-params.in” and “mat-cmnd.out”. Qubiter will multiply the SEO listed in “mat-cmnd.out”, arrive at a unitary matrix U_final, and output a file called “mat-check.out” which lists the entries of U_final.

We should also mention a small frill to the first mode of operation. At the same time that Qubiter outputs the file “mat-cmnd.out”, it also outputs a file called “mat-pict.out” which is a translation of “mat-cmnd.out” to a pictorial language. Each elementary gate is represented by a line in “mat-pict.out”. Consider a single line of the file. There is a 1 to 1 correspondence between the characters in the line and the qubits of an array of qubits. The rightmost character represents bit 0. The next to rightmost character in the line represents bit 1. And so on. A bit that is not operated on is represented by a “|” character. A bit that is rotated about the Z axis (ditto, Y axis) is represented by a “Z” character (ditto, “Y” character). If the gate is a c-not that flips a bit when the control bit is true (ditto, false), then the control bit is represented by a “@” character (ditto, “0” character). The bit to be flipped is represented by an “X” character. If the gate is a pure phase acting on all bits, all bits are represented by the ˜ character.

The CS Decomposition is intimately related to the Generalized Singular Value Decomposition (GSVD). In fact, Qubiter1.0 does CS decompositions by means of a Clapack subroutine for doing GSVD. For more information about the GSVD and its connection to the CS Decomposition, see Pai94 and references therein. See also the Clapack documentation that comes with the subroutine zggsvd.c

So far, we have described version 1.0 of Qubiter. Future versions of Qubiter are planned that will: (1)Take as input an arbitrary QB net (not just 2 connected nodes), and return as output a SEO. (2) Add quantum error correction code to the input QB net. (3)Include optimizations enabling it to produce SEOs with fewer steps.

In classical computation, the basic set of elementary operators is not unique. For example, instead of using AND, NOT and OR gates, one can just use NAND gates. The same is true in quantum computation: the basic set of elementary operators is not unique. In this preferred embodiment of the invention, we use the set {CONTROLLED-NOT, QUBIT ROTATION} of elementary operators, but the invention also applies to other sets of elementary operators.

So far, we have described what are at present the preferred embodiments of this invention. Those skilled in the art will be able to come up with many modifications to the given embodiments without departing from the present invention. It is therefore desired that the scope of the invention be determined by the appended claims rather than by the given embodiments.

Claims

1. A method of operating a classical computer, wherein said method must be stored in a computer readable medium which said classical computer can read, , wherein said method must be stored in the external or internal memory units of said classical computer, to calculate a sequence of operations with the purpose of applying said sequence of operations to a quantum computer comprising an array of qubits, to induce said quantum computer to execute a desired calculation, said method comprising the steps of:

storing in said classical computer a qb net (quantum Bayesian net) data-set comprising:

(a) graph information comprising a node label for each node of a plurality of n nodes, and also comprising a plurality of directed lines, wherein a directed line comprises an ordered pair of said node labels, wherein one member of the label pair labels the source node and the other member labels the destination node of the directed line,

(b) state information comprising, for each jε{1, 2, . . . n}, a finite set Σ_j containing labels for the states that the j'th node {circumflex over (x)}_j may assume assumes, and

(c) amplitude information comprising, for each jε{1, 2, . . . n}, a representation of a complex number $A_j\left[x_j❘x_{k1},x_{k2},\dots ,x_{kS_j}\right]$

2. The method of claim 1, wherein said sequence of operations comprises c-not operations elementary operations on said array of qubits.

3. The method of claim 1, wherein said sequence of operations comprises qubit rotation operations is a sequence of elementary operations on said array of qubits.

4. The method of claim 1, further also utilizing a quantum computer, comprising the additional step of:

manipulating said quantum computer largely according to said sequence of operations.

5. The method of claim 1, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequence of operations, wherein if the decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies specifies a represetation of three unitary matrices: L, D, R, such that the following matrix equation holds: U=LDR.

6. The method of claim 1, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequence of operations, wherein if the decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies eight unitary matrices: specifies a represetation of four matrices D₀₀, D₀₁, D₁₀, D₁₁, and four unitary matrices L₀, L₁, R₀, R₁, such that if one partitions said U into four blocks U₀₀, U₀₁, U₁₀, U₁₁, then the following four matrix equations hold: U_ij=L_iD_ijR_j for iεBool and jεBool.

7. The method of claim 1, comprising the additional step of:

calculating with said classical computer and using said qb net data-set, a tree data-set that comprises data that can be represented as a tree-node matrix for each node contained in a subset of the nodes of a tree, wherein the product, in some order defined by the determined in accordance with said tree, of all said tree-node matrices is equivalent to said v_in.

8. A method of operating a classical computer having display, storage and calculation means, wherein said method must be stored in a computer readable medium which said classical computer can read, , wherein said method must be stored in the external or internal memory units of said classical computer, to analyze a quantum physical system that exhibits quantum mechanical behavior comprising an array of qubits, said method comprising the steps of:

displaying on said display means a graph comprising a plurality of n nodes, and a plurality of directed lines connecting certain pairs of said nodes,

storing in said storage means a qb net (quantum Bayesian net) data-set comprising:

(a) graph information comprising a node label for each of said n nodes, and also comprising, for each said directed line, said node label for the source node and for the destination node of the directed line,

(b) state information comprising, for each jε{1, 2, . . . n}, a finite set Σ_j containing labels for the states that the j'th node {circumflex over (x)}_j may assume assumes, and

(c) amplitude information comprising, for each jε{1, 2, . . . n}, a representation of a complex number $A_j\left[x_j❘x_{k1},x_{k2},\dots ,x_{kS_j}\right]$

9. The method of claim 8, wherein said sequence of operations comprises c-not operations elementary operations on said array of qubits.

10. The method of claim 8, wherein said sequence of operations comprises qubit rotation operations is a sequence of elementary operations on said array of qubits.

11. The method of claim 8, further also utilizing a quantum computer, comprising the additional step of:

manipulating said quantum computer largely according to said sequence of operations.

12. The method of claim 8, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequence of operations, wherein if the decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies specifies a represetation of three unitary matrices: L, D, R, such that the following matrix equation holds: U=LDR.

13. The method of claim 8, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequences sequence of operations, wherein if the decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies eight unitary matrices; specifies a represetation of four matrices D₀₀, D₀₁, D₁₀, D₁₁, and four unitary matrices L₀, L₁, R₀, R₁, such that if one partitions said U into four blocks U₀₀, U₀₁, U₁₀, U₁₁, then the following four matrix equations hold: U_ij=L_iD_ijR_j for iεBool and jεBool.

14. The method of claim 8, comprising the additional step of:

calculating with said classical computer and using said qb net data-set, a tree data-set that comprises data that can be represented as a tree-node matrix for each node contained in a subset of the nodes of a tree, wherein the product, in some order defined by the determined in accordance with said tree, of all said tree-node matrices is equivalent to said v_in.

15. A method of operating a classical computer, wherein said method must be stored in a computer readable medium which said classical computer can read, , wherein said method must be stored in the external or internal memory units of said classical computer, to calculate a sequence of operations with the purpose of applying said sequence of operations to a quantum computer comprising an array of qubits, to induce said quantum computer to execute a desired calculation, wherein said classical computer comprises a unitary matrix decomposer, wherein if said unitary matrix decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U of dimension greater than 2, then the decomposer will calculate a data-set that fairly directly specifies specifies a represetation of three unitary matrices: L, D, R, such that the following matrix equation holds: U=LDR, wherein L and R each yields unitary matrices whose dimension is smaller than that of U, said method comprising the steps of:

storing in said classical computer an input data-set that fairly directly specifies specifies a represetation of a unitary matrix v_in, wherein at least one row of v_in has 3 or more non-zero entries and at least one column of v_in has 3 or more non-zero entries,

applying said unitary matrix decomposer to decompose the initial unitary matrix U=V_in,

applying said unitary matrix decomposer repeatedly to decompose initial unitary matrices obtained from the output of a previous application of said unitary matrix decomposer,

calculating with said classical computer, using said input data-set and data obtained by applying said unitary matrix decomposer, said sequence of operations on said array of qubits, wherein said sequence of operations and said v_in both would, if applied to an said array of qubits, produce equivalent transformations of the array of qubits.

16. The method of claim 15, wherein if said unitary matrix decomposer is given data that fairly directly specifies a unitary matrix U of dimension greater than 2 species a represetation of said initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies eight unitary matrices: specifies a represetation of four matrices D₀₀, D₀₁, D₁₀, D₁₁, and four unitary matrices L₀, L₁, R₀, R₁, such that if one partitions said U into four blocks U₀₀, U₀₁, U₁₀, U₁₁, then the following four matrix equations hold: U_ij=L_iD_ijR_j for iεBool and jεBool.

17. The method of claim 15, wherein said sequence of operations comprises c-not operations.

18. The method of claim 15, wherein said sequence of operations comprises qubit rotation operations elementary operations on said array of qubits.

19. The method of claim 15, wherein some of said sequences of operations comprise qubit rotations and c-nots said sequence of operations is a sequence of elementary operations on said array of qubits.

20. The method of claim 15, further also utilizing a quantum computer, comprising the additional step of:

manipulating said quantum computer largely according to said sequence of operations.

21. A method of operating a classical computer, wherein said method must be stored in a computer readable medium which said classical computer can read, , wherein said method must be stored in the external or internal memory units of said classical computer, to calculate a sequence of operations with the purpose of applying said sequence of operations to a quantum computer comprising an array of qubits, to induce said quantum computer to execute a desired calculation, said method comprising the steps of:

storing in said classical computer an input data-set that fairly directly specifies specifies a represetation of a unitary matrix v_in, wherein at least one row of v_in has 3 or more non-zero entries and at least one column of v_in has 3 or more non-zero entries,

calculating with said classical computer and using said input data-set, a tree data-set that comprises data that can be represented as a node matrix M_j for each node j contained in a subset J of the nodes of a tree, wherein the product, in some order defined by the determined in accordance with said tree, of all said node matrices is equivalent to said v_in,

calculating with said classical computer and using said tree data-set, for each of the node matrices M, a string of operations on qubits, wherein said string of operations and said M both would, for each jεJ, a product π_j of operations on said array of qubits, wherein M_j and π_j both would, if applied to an said array of qubits, produce equivalent transformations of the array of qubits.

22. The method of claim 21, wherein some of said sequences of operations comprise for some jεJ, π_j comprises c-not operations.

23. The method of claim 21, wherein some of said sequences of operations comprise for some jεJ, π_j comprises qubit rotation operations.

24. The method of claim 21, wherein some of said sequences of operations comprise qubit rotations and c-nots for some jεJ, π_j comprises elementary operations on said array of qubits.

25. The method of claim 21, wherein said sequences of operations are sequences of elementary operations for all jεJ, π_j is a sequence of elementary operations on said array of qubits.

26. The method of claim 21, further also utilizing a quantum computer, comprising the additional step of:

manipulating said quantum computer largely according to said sequence of operations.

27. The method of claim 21, wherein if said v_in is a square matrix with n rows, then said node matrices are also square matrices with n rows and with at most two non-zero entries in each row and in each column.

28. The method of claim 21, wherein if said v_in is a square matrix with n rows, then said node matrices are also square matrices with n rows and with about two non-zero entries in each row and in each column.

29. The method of claim 21, wherein if said v_in is a square matrix with n rows, then said node matrices have theoretically a number of nonzero real parameters which grows linearly with n or slower, for large n grows, for large n₁ linearly with n or slower.

30. The method of claim 21, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequences sequence of operations, wherein if said unitary matrix decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies specifies a represetation of three unitary matrices: L, D, R, such that the following matrix equation holds: U=LDR.

31. The method of claim 21, wherein said tree is a binary tree in which each node branches out into two nodes.

32. The method of claim 31 21, wherein said classical computer includes a unitary matrix decomposer which is used in the calculation of said sequences sequence of operations, wherein if the decomposer is given data that fairly directly specifies a specifies a represetation of an initial unitary matrix U, then the decomposer will calculate a data-set that fairly directly specifies eight unitary matrices: specifies a represetation of four matrices D₀₀, D₀₁, D₁₀, D₁₁, and four unitary matrices L₀, L₁, R₀, R₁, such that if one partitions said U into four blocks U₀₀, U₀₁, U₁₀, U₁₁, then the following four matrix equations hold: U_ij=L_iD_ijR_j for iεBool and jεBool.

33. The method of claim 32, wherein some of said sequences of operations comprise for some jεJ, π_j comprises c-not operations.

34. The method of claim 32, wherein some of said sequences of operations comprise for some jεJ, π_j comprises qubit rotation operations.

35. The method of claim 32, wherein some of said sequences of operations comprise qubit rotations and c-nots for some jεJ, π_j comprises elementary operations on said array of qubits.

36. The method of claim 32, wherein said sequences of operations are sequences of elementary operations for all jεJ, π_j is a sequence of elementary operations on said array of qubits.

37. The method of claim 32, further also utilizing a quantum computer, comprising the additional step of:

manipulating said quantum computer largely according to said sequence of operations.