Arithmetic Circuits for Multilevel Qudits Based on Quantum Fourier Transform
Archimedes Pavlidis Emmanuel Floratos
quant-ph cs.ET

    

Arithmetic Circuits for Multilevel Qudits  Based on Quantum Fourier Transform
A. Pavlidis and E. Floratos

empty










1 




ARITHMETIC CIRCUITS FOR MULTILEVEL QUDITS


BASED ON QUANTUM FOURIER TRANSFORM 





ARCHIMEDES PAVLIDIS


Department of Informatics, University of Piraeus
=10pt
80, Karaoli & Dimitriou str., GR 185 34, Piraeus, Greece
Department of Informatics and Telecommunications, National and Kapodistrian University of Athens
=10pt
Panepistimiopolis, Ilissia, GR 157 84, Athens, Greece
=10pt
e-mail: adp@unipi.gr


EMMANUEL FLORATOS


Department of Physics, National and Kapodistrian University of Athens
=10pt
Panepistimiopolis, Ilissia, GR 157 84, Athens, Greece
=10pt
Institute of Nuclear and Particle  Physics, N.C.S.R. Demokritos
=10pt
27, Neapoleos str., Agia Paraskevi, GR 153 41, Athens, Greece
=10pt
e-mail: mflorato@phys.uoa.gr














We present some basic integer arithmetic quantum circuits, such as adders and multipliers-accumulators of various forms, as well as  diagonal operators,  which operate on multilevel qudits. The integers to be processed are represented in an alternative basis after they have been Fourier transformed. Several arithmetic circuits operating on Fourier transformed integers have appeared in the literature for two level qubits. Here we extend these techniques on multilevel qudits, as they may offer some advantages relative to qubits implementations. The arithmetic circuits presented can be used as basic building blocks for higher level algorithms such as quantum phase estimation, quantum simulation, quantum optimization etc., but they can also be used in the implementation of a quantum fractional Fourier transform as it is shown in a companion work presented separately.  












   







§ INTRODUCTION
     
  
The common representation of the elementary quantum information is the qubit, where its state is a superposition a|0⟩+b|1⟩ which belongs to a two-dimensional Hilbert space with two basis states |0⟩ and |1⟩ known as the computational basis. A quantum computer is a finite dimensional quantum system composed of a  qubits collection, performing various unitary operations on the qubits (quantum gates) and quantum measurements. Accordingly, there is a correspondence between a qubit and a classical bit, in the sense  that the basis states of a qubit follow the binary logic. We can extend this correspondence to multivalued logic instead of two values only by enlarging the dimension of the elementary Hilbert space used. The qudit is a generalization of the qubit to a larger Hilbert space of dimension d>2. The state of a qudit is a superposition a_0|0⟩+a_1|1⟩+⋯ + a_d-1|d-1⟩, where |0⟩, |1⟩, …  |d-1⟩ are the computational basis states. Qutrit is a special name for the case  d=3, while ququart corresponds to d=4. In many cases, the employment of a multivalued quantum logic is more natural. E.g. in ion traps we could exploit more than two energy levels. Multiple laser beams could be used to manipulate the transitions between these levels <cit.>. 

Working with qudits instead of qubits may offer some advantages. The required number of qudits is smaller by a factor log _2d than the corresponding number of qubits for the same dimension a quantum computer has to explore. E.g. the dimensions of a composite system of n qubits  is 2^n, while the same dimension  can be reached with only log_d2^n=log_22^n / log_2d=n/log_2d qudits. Such as reduction of the required number of physical carriers of quantum information is advantageous, considering the difficulty of reliably controlling a large number of  carriers. Also, when fewer quantum information carriers are used, a decrease in the overall decoherence is expected and this fact favors the scalability issues <cit.>. 

Another advantage, which is also related to the adverse effect of decoherence, is that fewer multilevel qudit gates are required to construct a quantum circuit implementing a given unitary operation compared to the case of using two-level gates <cit.>.  
Fewer gates reduce the number of steps needed to complete the circuit operation (circuit depth), and consequently less errors are accumulated during the overall operation of the circuit. Even so, protection of quantum information against environmental interaction is inevitable. Quantum error correcting codes and fault tolerant gate constructions to combat decoherence on multilevel qudits have been proposed and they are similar to the ones used for the qubit case <cit.>. 

At a higher level, generalizations of known quantum algorithms and circuits using d-level qudits may offer improvements with respect to their qubits implementation counterparts. E.g., quantum phase estimation, which is the core part of Shor's algorithm <cit.>	 and also it is used in quantum simulation <cit.>, is improved in terms of success probability when multilevel qudits are incorporated <cit.>. Multiple-valued version of Deutsch-Josza algorithm  has been reported in <cit.> while an implementation proposal for five level superconducting qudit appeared in <cit.>. Qudits version for Grover's algorithm <cit.> has been reported in <cit.>. The high dimensional Deutsch-Josza algorithm may find applications in image processing, while the high dimensional Grover's algorithm offers a trade-off between space and time. 

An assortment of quantum gates operating on qudits have been proposed or experimentally realized on various technologies. Single and two qudit d-level gates proposed in <cit.> for the ion trap technology. Single qudit gates for d=5 implemented in superconducting technology and used to emulate spins of 1/2, 1 and 3/2 in <cit.>. Proposals for single and two qudit gates based on superconducting technology appeared in <cit.>. Single qudit gates based on optical technology reported in <cit.>. Three dimensional entanglement between photons observed in <cit.>.

In this work we present some quantum arithmetic circuits operating on d-level qudits by extending results given in prior works <cit.>. These circuits exploit the quantum Fourier transform and various single qudit and two qudit rotation gates to perform the desired calculations. Processing in the Fourier domain may offer some advantages related to speed <cit.> and robustness to decoherence <cit.>. Among the proposed circuits are various versions of adders (adder with constant, generic adder, adder with constant controlled by single qudit) and multipliers  (multiplier with constant and accumulator, multiplier with constant). Such circuits are useful in many quantum algorithms, e.g. quantum phase estimation, quantum simulation. 

The increased interest in quantum information processing exploiting d-level qudits, both in theoretical and experimental aspects, was one of the stimulation for this work. 
However, the main motivation was the particular application targeted by the quantum circuits presented in this manuscript, which is a new definition of the fractional Fourier transform. Unitary operations on high dimensional d-level qudits fit more naturally for this specific application because the proposed fractional Fourier transform operates on a Hilbert space of dimension d^n, where d ≠ 2 is a prime. The development of the quantum fractional Fourier transform and its implementation on qudits is presented in a separate work <cit.>.

The rest of the paper is organized as follows: A short background about design and synthesis of qudits quantum circuits is given in section <ref>. The elementary and basic qudit gates used in the proposed designs are given in section <ref>. The quantum Fourier transform definition and its circuit for q qudits of d levels is presented in section <ref>. Section <ref> introduces integer arithmetic circuits like adder with constant, adder of two integers, controlled adder with constant, multiplier with constant and accumulator, and multiplier with constant. All of the arithmetic units accept one of their operands after it has been Fourier transformed.  In section <ref> a method to implement a diagonal operator on q qubits is analyzed where the diagonal elements are some powers of roots of unity. A quantum multiplier of two integers is introduced  and then a quantum squarer is built upon this multiplier. It is demonstrated how to introduce relative phases between the basis states of superposition which depend quadratically on the index of the basis state. It can be generalized for a function that is polynomial in the states index. This operation is a necessary part of  the quantum fractional Fourier transform presented in the companion article and also it may find other applications, such as quantum simulation algorithms and Grover's search algorithm. Appendix A gives the decomposition of a three qudits rotation gate introduced in section <ref>. Complexity analysis in terms of quantum cost, depth and width is reported in section <ref>. In Appendix B we discuss how it is possible to use a discrete library of components to approximate the proposed designs and the impact to cost and depth. A discrete library of gates is necessary if fault tolerance is to be incorporated. Finally, we conclude in section <ref>.



§ BACKGROUND AND RELATED WORK


The construction of a complex quantum circuit operating on multilevel qudits is based on the selection of a set of elementary qudit gates and their interconnection so as to achieve the target operation. A multilevel gate operates on a single qudit, on two qudits or more qudits. A single d-level qudit gate is represented by a unitary matrix U of dimensions d× d. It transforms an initial qudit state |ψ⟩ _in to |ψ⟩ _out=U|ψ⟩ _in. As an example consider the  general superposition qutrit state |ψ⟩ _in=a|0⟩+b|1⟩+c|2⟩. The application of the gate 

    U=
    [ 0 0 1; 1 0 0; 0 1 0 ]
    =|1⟩⟨ 0| + |2⟩⟨ 1|  + |0⟩⟨ 2|


on this state results to the state |ψ⟩ _out= c|0⟩+a|1⟩+b|2⟩.

Two qudit gates operate on states of two qudits which are d^2 dimensional, so their representing unitary matrices have dimensions of d^2 × d^2. Two single qudit gates V_1 and V_2 operating on two different qudits can be seen as a two qudit gate which is their tensor product U=V_1⊗ V_2. However, not every two qudit gate can be decomposed as a tensor product of two single qudit gates, in which case we have an entangling gate.  

Consecutive application of single, two or more qudit gates to a collection of q qudits results in a quantum circuit which is represented by a unitary matrix of dimensions d^q × d^q. The design of a quantum circuit is the procedure of interconnecting various elementary gates so as to fulfill the given specifications. These specifications are given in the form of a unitary matrix or the relationship between the desired input-output state relation in the computational basis.

It is proven that single qudit gates  and a two qudit gate alone are adequate to form a universal set of gates, provided that the two qudit gate is an entangling gate <cit.>. A universal set of gates can be used to approximate any target  quantum circuit with arbitrary precision. Various sets of qudit gates (gate libraries) and methods to exploit them to build more complex unitaries have been introduced in the literature. The library used in  <cit.> consists of single and two qudit gates with continuous parameters and the synthesis method is based on spectral decomposition of the target unitary matrix. Cosine-Sine decomposition is another method used in  <cit.>. A discrete set of single qudit gates and a single two qudit gate is used in  <cit.>  to synthesize the large unitary matrix using QR decomposition. In <cit.> a different two qudit gate and a set of single qudit gates is used along with quantum Shannon decomposition to synthesize the target unitary. The previous methods and results are similar to the two-level qubits synthesis cases.  It is proven that the cost of the resulting circuit in terms of two qudit gates is  upper bounded by O(d^2n) where n is the qudits number <cit.>. Thus, these automated methods are suitable only for small quantum circuits due to the exponential cost increase. 

When the target circuit is an arithmetic or logic block where its unitary matrix is a permutation matrix consisting of 0 and 1 elements, then multiple-valued reversible synthesis methods could be applied. These methods are extension of the binary reversible logic case and may be applied to a specific value of d, e.g. <cit.> (d=3), or applied to any value of d <cit.>. Similarly to the quantum synthesis case, these algorithms are not suitable for large circuits.

As many algorithms widely use quantum arithmetic blocks like adders or multipliers recurrently, it is crucial to have available efficient arithmetic and logic blocks. Ad hoc design of such blocks usually offers better results compared to the automated synthesis methods. One can exploit the iterative and regular structure of these arithmetic blocks or extend known classical designs to the quantum case. A diversity of ad hoc designed quantum arithmetic and logic circuits for two-level qubits can be found in the literature <cit.>, but few (usually adders) are known for multilevel qudits and also they are mostly designed for a specific value of d. In contrast, the proposed designs are parametrized for any value of d.

One of the first ternary (d=3) quantum adder for 3-inputs only appeared in <cit.> as an application example of the proposed synthesis method. Ternary quantum adders/subtractors ad hoc designed for any number of inputs is given in <cit.>. In <cit.> a ternary extension of the well known VBE ripple-carry adder <cit.> is reported. Quaternary (d=4) comparators proposed in <cit.>. Improved designs of ternary ripple carry and carry look-ahead adders along with modifications that lead to subtractors and comparators are given in  <cit.>. The previous ternary ripple carry adder is a modification of the CDKM binary quantum adder appeared in  <cit.>  and it has also a depth of O(n) using one ancilla qutrit. Similarly, the previous ternary carry look-ahead quantum adder is an extension of the DKRS binary quantum adder appeared in <cit.> and thus it offers a depth O(log(n)) using O(n) ancilla qutrits. 

Several of the previous multilevel qudits (usually qutrits) designs are modifications of binary quantum adders. The gates libraries used are differentiated among each design. This is justifiable as the implementation technology for multilevel qudits is far apart to be considered matured. However, gates of one library can be expressed as gates of another one, provided that the libraries are universal. 

Diagonal operator circuits on qubits or qudits don't change the absolute value of the amplitudes of a superposition, but rearrange their relative phases. Such circuits are useful in quantum algorithms <cit.> like quantum optimization, quantum simulation, Grover's search etc.  Synthesis methods for diagonal unitary matrices of two-level qubits  have been developed  for diagonals of special structure <cit.> or for any diagonal <cit.>. Recently, diagonal synthesis for multilevel qudits reported in <cit.>. In general, the synthesis cost is related to the dimensionality of the Hilbert space covered (that is exponential in the number of the qubits or qudits) and the number of the distinct phases of the diagonal. In this work, based on ideas of <cit.> and the arithmetic circuits designed,  we develop a diagonal operator circuit which has a special structure, that is the phases are quadratic functions of the coordinates, with polynomial cost and depth. Using same techniques, other powers, instead of quadratic, can be achieved.

The gates used in our proposed designs are the ones introduced in <cit.> where also physical implementation directions are given. The proposed designs use extensively various rotation gates. As the rotation angles of the gates vary with the size of the circuit and also small angles are required, implementation and fault tolerance issues are addressed in Appendix B, using results of the binary quantum case. 

Many of the arithmetic multilevel quantum designs of this manuscript are direct extension or modifications of binary quantum designs appeared in  <cit.> which use QFT before applying the rotation gates to one of the two integer operands and then applying the inverse QFT to bring back the result in the computational basis. The following arithmetic units are presented:



  * Adder of two integers of q qudits with depth of O(q) and width of 2q qudits  (Subsection <ref>).


  * Adder of an integer of q qudits with a constant integer. Its depth is O(q) including the direct and inverse QFT blocks or O(1) without the QFT blocks. Its  width is q qudits  (Subsection <ref>).


  * Single state controlled adder of an integer of q with a constant. The adder is enabled if the control qudit is in a particular basis state of the different d possible states, otherwise it acts as an identity. Its depth is O(q) and its width is q+1 qudits  (Subsection <ref>).


  * Generalized controlled adder of an integer of q qudits with a constant. It adds a multiple of the constant to the integer. The multiple depends on the state of control qudit, being between 0 and d-1. Its depth is O(q) and its width is q+1 qudits (Subsection <ref>).


  * Multiplier with constant and accumulator.  It multiplies an integer of q qudits with a constant and adds the product to a second integer of q qudits. Its depth is O(q) and its width is 2q qudits (Subsection <ref>).


  * Multiplier with constant. It multiplies an integer of q qudits with a constant provided that the constant is relative prime with d^q, which is always the case when p is prime.  Its depth is O(q) and its width is 2q qudits of which q qudits are ancilla initialized to the zero state and then are reset back to the zero state (Subsection <ref>).


  * Multiplier of two integers and accumulator. It multiplies two integer of q qudits and adds the product to a third integer of q qudits. Its depth is O(q^2)  and its width is 3q qudits (Subsection <ref>).


  * Squarer/Multiplier with constant/Accumulator. It performs the transform |x⟩ |z⟩→ |x⟩  |z + γ x^2⟩, where γ is the integer constant. Its depth is O(q^2) and its width is 4q qudits   (Subsection <ref>).


  * General diagonal operator. It operates diagonally on a general superposition state of  q qudits and changes the phases of the superposition amplitudes by applying the matrix  
∑_k=0^d^q-1 e^i2π/d^qf(k) |k⟩⟨ k|, where f(k) is a function of k. The specific diagonal operator presented here is based on the previous squarer and some other blocks as it applies the function f(k)=γ k^2. It can be generalized for other powers of k or even for polynomial functions on k by exploiting similar techniques. It has a depth of O(q^2) and its width is 4q qudits of which 3q qudits are ancilla (Section <ref>). 


Detailed complexity analysis in terms of quantum cost and depth is given in section <ref>, where the parameter d enters the previous rough approximations. This is because many gates like the basic rotation gates used and introduced in subsection <ref> are synthesized using more elementary gates with a cost (and consequently depth) which depends on the dimension d of the qudits.



§ ELEMENTARY AND BASIC GATES ON QUDITS

We followed a hierarchical bottom-up approach to design the arithmetic circuits. At the lowest level, elementary gates operating in a two dimensional subspace of the d-dimensional space of a qudit are used. Upon them, more complex gates (which are basic for the designs) operating in the whole d-dimensional space are built. Some of the basic gates are reported in <cit.>, while others like the generalized controlled and doubly controlled rotation gates are introduced here (subsection <ref>, subsection <ref>  and Appendix A ). 





 §.§ Generalized X gates

The X^(jk) gates <cit.> operate on a two-dimensional subspace of a d-level qudit by exchanging the basis states |j⟩,|k⟩,   and leaving intact the other basis states, thus they are a generalization of the well  known X gate for qubits which exchanges the basis states   |0⟩ and |1⟩. They are defined by the d× d matrix 


    [ X^(jk)=|j ⟩⟨ k| + |k ⟩⟨ j| +
    ∑_l
    n=0 
    
    n ≠ j 
    
    n ≠ k
    ^d-1 |n ⟩⟨ n|                                                        j,k=0 … d-1 ]


  
It holds that X^(jk)= X^(kj), so there are d(d-1)/2 different such gates in this family.



 §.§ Rotation gates of two levels

These gates perform a rotation on a two dimensional subspace <cit.> of a d-level qudit and are defined as 


    [ R_a^jk(θ)=exp(-iθσ _a^(jk)/2   ),                    0 ≤ j,k ≤ d-1,                        a ∈{x,y,z} ]


where σ_x^(jk)=|j⟩⟨ k| +  |k⟩⟨ j|,
    σ_y^(jk)=-i|j⟩⟨ k| +  i|k⟩⟨ j| and
    σ_z^(jk)=|j⟩⟨ j| - |k⟩⟨ k| for j,k=0… d-1 are 
matrices of dimensions d × d. Parameter θ is the rotation angle, while i=√(-1). 



 §.§ Generalized Controlled X gates

The GCX_(m)^(jk) gates are generalization in the qudits of the CNOT gates acting on qubits <cit.>. Thus, they are gates which operate on a control and a target qudit. A GCX gate has three parameters, m, j and k, which define its operation. A GCX_(m)^(jk) acts like a
X^(jk) on the target qudit iff the control qudit is on the basis state |m ⟩. Consequently, the definition matrix of such a gate is block diagonal with dimension d^2× d^2 consisting of d  blocks of d× d dimensions each and it is given by


    [ GCX_(m)^(jk)= |m ⟩⟨ m| ⊗( |j ⟩⟨ k| + |k ⟩⟨ j| +
     ∑_l
     n=0 
    
     n ≠ j 
    
     n ≠ k
     ^d-1 |n ⟩⟨ n| ) +
     ∑_l
     n=0 
    
     n ≠ m 
     ^d-1 |n ⟩⟨ n| ⊗ I_d                                                                                                                           j,k,m=0 … d-1 ]


where I_d is the identity matrix of dimensions d× d. Equation (<ref>) can be equivalently written as


    GCX_m^(jk)= diag(I_d,I_d,…,  m-th blockX^(jk) ,…,I_d  )




 §.§ Hadamard gate


The Hadamard gate H^(d)  on d-level qudits is defined by the matrix


    H^(d)=    1/√(d)[                1                1                1                1                 ;                1         e^i2π1/d         e^i2π2/d                ⋯       e^i2πd-1/d;                1      e^i2π 2 1/d      e^i2π 2 2/d                ⋯    e^i2π 2 d-1/d;                ⋮                ⋮                ⋮                ⋱                ⋮;                1    e^i2π(d-1)1/d   e^i2π (d-1)2/d                ⋯ e^i2π (d-1)d-1/d ] =  
       1/√(d)[                  1                  1                  1                  1                   ;                  1         e^i2π(0.1)         e^i2π(0.2)                  ⋯       e^i2π(0.d-1);                  1       e^i2π 2(0.1)      e^i2π 2 (0.2)                  ⋯    e^i2π 2 (0.d-1);                  ⋮                  ⋮                  ⋮                  ⋱                  ⋮;                  1    e^i2π(d-1)(0.1)   e^i2π (d-1)(0.2)                  ⋯ e^i2π (d-1)(0.d-1) ]



In the above equation the notation (0.n) is the fractional representation of n/d in the base-d arithmetic system. The application of the H^(d) gate  to a basis state |j ⟩  is shown below


    H^(d)|j ⟩   =    1/√(d)[               1      e^i2π(0.j)    e^i2π 2(0.j)               … e^i2π(d-1)(0.j) ]
    ^T = 
       1/√(d) (   |0⟩+e^i2π(0.j)|1⟩+⋯ +e^i2π(d-1)(0.j) |d-1⟩   )



The Hadamard gate for qudits essentially performs the order-d Fourier transform, likewise the Hadamard gate for qubits performs the order-2 Fourier transform. Methods for implementation of the H^(d)  gate are proposed in <cit.>.

The symbols that will be used throughout the text for the three families of elementary gates defined and the H^(d)  gate are shown in Figure <ref>. 





 §.§ Diagonal Gates of one and two qudits


The qudit elementary gates of the previous section affect a 2-dimensional subspace of the whole d-dimensional Hilbert space of a single qudit. In this section single and two qudit diagonal basic gates affecting the whole d-dimensional space of one of the qudits are described and synthesized using elementary gates of the previous section.



  §.§.§ Diagonal D^'(a_1,a_2,…,a_d-1) and D(φ_1,φ_2,…,φ_d-1) gates 



The diagonal D^'(a_1,a_2,…,a_d-1)  gate <cit.> is defined by the equation 


    D^'(a_1,a_2,…,a_d-1)=e^iφ
    diag( e^-i(a_1+a_2+… +a_d-1) , e^ia_1, e^ia_2, … , e^ia_d-1 )



It can be easily proved that such a gate can be constructed by sequentially applying d-1    R_z^(jk)(θ) gates as shown in the following equation 


    D^'(a_1,a_2,…,a_d-1)=e^iφ
    R_z^(01)(a_1) R_z^(02)(a_2) ⋯  R_z^(0(d-1))(a_d-1)



A related gate is the D(φ_1,φ_2,…,φ_d-1) defined as 


    D(φ_1,φ_2,…,φ_d-1) =
    diag(1,e^iφ_1,e^iφ_2,…, e^iφ_d-1 )



The D(φ_1,φ_2,…,φ_d-1)  gate is identical with the  D^'(a_1,a_2,…,a_d-1) gate if we set 


    [ a_j= φ_j-1/d∑_k=1^d-1φ_k                j=1 … d-1 ]


and add a global phase of angle φ=1/d∑_k=1^d-1φ_k to every diagonal element of 
D^'(a_1,a_2,…,a_d-1).



  §.§.§ Controlled Diagonal CD^'(a_1,a_2,…,a_d-1) and CD(φ_1,φ_2,…,φ_d-1) gates 


The diagonal gates of the previous subsection can be extended to operate on two qudits, where the first is the control qudit and the second is the target qudit, in the following manner: A diagonal gate  D^'(a_1,a_2,…,a_d-1)   or   D(φ_1,φ_2,…,φ_d-1)  is applied on the target qudit iff the control qudit is in state |m⟩  , otherwise no operation is effective on the target. Thus, the d^2× d^2 matrices representing such gates have the following block diagonal form 


    CD_m^'(a_1,a_2,…,a_d-1)=
    diag(I_d, … ,I_d,  m-th  blockD^'(a_1,a_2 ,…,a_d-1) ,I_d, … ,I_d  )


and


    CD_m(φ_1,φ_2,…,φ_d-1)=
    diag(I_d, … ,I_d,  m-th  blockD(φ_1,φ_2 ,…,φ_d-1) ,I_d, … ,I_d  )



A construction of a CD_m^'(a_1,a_2,…,a_d-1) gate using 4(d-1) elementary GCX_(m)^(jk) and  R_z^(jk)(θ) gates is shown in 
Figure <ref>. Single qudit gate  S_m=diag(1, … ,1, m-th pose^iφ,1,… , 1) is a phase gate which is identical to a D' gate up to a global phase.







 §.§ Generalized Controlled Rotation gate R_k^(d)  


The controlled diagonal gates CD_m^' and CD_m of the previous subsection are activated whenever the control state is equal to one of the d possible basis states, e.g. |m⟩. We define a basic controlled diagonal gate, R_k^(d), such that each one of the d possible control states have a different effect on the target qudit. Such gates will be useful in the QFT and arithmetic circuits presented in the following sections.  
The R_k^(d)  gate is parametrized by the integer k. The matrix defining this gate is block diagonal of the form 


    R_k^(d)=diag ( 
    ( Φ_k^(d)) ^0, 
    ( Φ_k^(d)) ^1,  … ,
    ( Φ_k^(d)) ^d-1)


where the matrix Φ_k^(d) is diagonal too, and defined with


    Φ_k^(d)=diag ( 
    1, 
    e^iφ_1,   
    e^iφ_2, … ,
    e^iφ_(d-1))



The angles φ_1,φ_2,…,φ_(d-1) depend on the parameter k as follows 


    [ φ_m=2π/d^km,    m=1,…,d-1 ]



The R_k^(d)  gates can be equivalently written in a more detailed form consisting of a sum of tensor products of the basis states of the two qudits as 


    R_k^(d)=
    ∑_j=0^d-1∑_m=0^d-1 
    e^i2π/d^kjm
    |j⟩⟨ j| ⊗ |m⟩⟨ m| =
    ∑_j=0^d-1∑_m=0^d-1 
    e^i2π (0.00… 0_k-1j)m
    |j⟩⟨ j| ⊗ |m⟩⟨ m|



We can see by inspecting Eq. (<ref>)  that an R_k^(d)  gate is a generalization on qudits of the controlled rotation gates R_k=R_z(2π/2^k)=diag(1,1,1,e^i2π/2^k)   for the qubit case (where d=2) and this generalization will be exploited when constructing the QFT and various arithmetic circuits based on the QFT. To understand this, it is useful to see what is the effect of an R_k^(d) gate when the control qudit is on a basis state |j_1⟩ (j_1=0,1,…,d-1) and the target qudit is in a superposition of equal amplitudes, but with different phases, such as  
|b⟩=1/√(d)∑_l=0^d-1e^iφ_l|l⟩. The joint state of the two qudits after the application of the R_k^(d) gate is 


    R_k^(d)( |j_1⟩ |b ⟩)=    1/√(d)∑_j=0^d-1∑_m=0^d-1 
    e^i2π (0.00… 0_k-1j)m
    |j⟩⟨ j|j_1⟩_=δ_jj_1⊗|m⟩⟨ m| ∑_l=0^d-1e^iφ_l|l⟩
    _=e^iφ_m|m⟩
     =  
       1/√(d)∑_m=0^d-1 
    e^i2π (0.00… 0_k-1j_1)m
    |j_1⟩ e^iφ_m |m⟩
     =  
       1/√(d) 
     |j_1⟩∑_m=0^d-1 
    e^i2π[ (0.00… 0_k-1j_1)m +φ_m] 
    |m⟩






Thus, an angle 2π  (0.00… 0_k-1j_1)m=2π/d^k j_1  m is added to every component |m⟩ of the target qudit superposition and this angle is proportional to the value |j_1⟩ of the control qudit and also proportional to the |m⟩ component of target qudit superposition.

The implementation of an R_k^(d) can be achieved by sequentially combining d-1 controlled diagonal gates CD_m(φ_1,φ_2,…,φ_d-1) for m=1 … d-1 and different angles for each case of m as shown below (see also Eqs. (<ref>),(<ref>) and (<ref>) ) 


    R_k^(d) =    
    CD_(1)( 
    2π/d^k,2π/d^k2, … , 2π/d^k(d-1)
    ) ·
    CD_(2)( 
    2π/d^k2,2π/d^k4, … , 2π/d^k2(d-1)
    ) ⋯
       
    CD_(d-1)( 
    2π/d^k(d-1),2π/d^k(d-1)2, … , 2π/d^k(d-1)(d-1)
    )




Taking into account that a  CD_(m)(φ_1,φ_2,…,φ_d-1) gate is composed by 4(d-1) elementary GCX_(m)^(jk) and R_z^(jk)(θ) gates, then we conclude that an R_k^(d) gate requires 4(d-1)^2 elementary gates.  The symbol used for the R_k^(d)  gate in this text is shown in Figure <ref>.









§ QUANTUM FOURIER TRANSFORM


The Quantum Fourier Transform on the N-dimensional computational basis  { |0⟩ ,|0⟩  , … ,|N-1 ⟩}  is defined by


    |j⟩1/√(N)∑_k=0^N-1 
    e^i2π/Njk|k⟩


Using q qudits of d levels <cit.>,<cit.>, and setting N=d^q, the q qudits basis consists of |j⟩ = |j_1… j_q⟩ = |j_1⟩… |j_q⟩ where for l-th qudit it holds |j_l⟩∈{ |0⟩ , … , |d-1 ⟩}. Then, the QFT action on a basis state |j⟩ (j=0 … d^q-1) is


    |j⟩ = |j_1… j_q⟩   1/√(N)∑_k_1=0^d-1⋯∑_k_q=0^d-1 
    e^i2π/d^qj ∑_l=1^qk_ld^q-l
    |k_1… k_q⟩ 
    = 
       ∑_k_1=0^d-1⋯∑_k_q=0^d-1⊗_l=1^q 
    e^i2π jk_ld^-l |k_l⟩
    = 
       ⊗_l=1^q∑_k_l=0^d-1  
    e^i2π jk_ld^-l |k_l⟩
    = 
    
    ( 
    ∑_m=0^d-1
    e^i2π (0.j_q) m |m⟩)
    ( 
    ∑_m=0^d-1
    e^i2π (0.j_q-1j_q) m |m⟩)
    ⋯( 
    ∑_m=0^d-1
    e^i2π (0.j_1j_2… j_q-1j_q) m |m⟩)



The d-ary representation (j_1j_2… j_q) of the integer j=j_1d^q +j_2d^q-1+ ⋯ +j_q as well as the fractional d-ary representation (0.j_1j_2…j_q)= j_1/d +j_2/d^2+ … +j_q/d^q are used in the above definition. This tensor product form is similar to the form of the QFT of order 2^n implemented using n qubits of two levels. Thus, the structure of a QFT circuit implemented with qudits is similar to the binary QFT case as depicted in Figure <ref>.


Indeed, comparing the state ∑_m=0^d-1
e^i2π (0.j_lj_l+1… j_q-1j_q) m |m⟩ 
of the l-th qudit after the transformation of Eq.(<ref>) with Eq.(<ref>) and (<ref>) we can conclude that this state can be generated by applying at the basis state |j_l⟩  of the l-th qudit a Hadamard gate H^(d)  and a sequence of q-l generalized rotation gates R_k^(d) , with k=2… q-l+1,  controlled by the qudits l+1 … q, respectively. At the end, the order of the qudits must be reversed with swap gates as in the case of the QFT operated on qubits. This swapping of the qudits is not shown in Figure <ref>. The inverse QFT circuit is derived by reversing horizontally  the direct QFT circuit of Figure <ref> (including the SWAP gates not shown) with  opposite signs in the angles of the rotation gates.








§ ARITHMETIC CIRCUITS


The integer arithmetic circuits presented in this section are developed in a bottom up succession, starting from the simpler ones and proceeding gradually to more complex ones. The arithmetic operation are assumed to be modulo d^q where d are the qudit levels and q is the number of qudits used to represent the integers. All the adders can be easily converted to subtractors  by using opposite sign in the angles of the rotation gates while retaining the same circuit structure.



 §.§ Adder of two integers (ADD)


A basic arithmetic operation block is an adder of two integers of q d-ary digits each, e.g a=(a_1a_2… a_q) and b=(b_1b_2… b_q) or two superpositions of integers. Following  the previous sections, the most significant d-ary digit of an integer is indexed with 1 while the least significant digit is indexed with q.
The circuit in Figure <ref>  operates on 2q qudits, the state |b_1… b_q⟩  of the q upper qudits (upper register) represents integer b while the state of the lower q qudits (lower register) represents the Fourier transformed state of the other integer a, that is |φ_1(a)⟩ |φ_2(a)⟩⋯ |φ_q(a)⟩, where |φ_l(a)⟩= 
∑_m=0^d-1e^i2π ( 0.a_la_l+1… a_q)m (see Eq. (<ref>)). It is a generalization on qudits of the circuit proposed in <cit.>.
 



The first qudit of the lower register is initially in the state |φ_1(a)⟩. The effect of the first rotation gate R_1^(d) controlled by state |b_1⟩ to this qudit (step[1,1]), taking into account Eq. (<ref>), is to evolve it  in the state 


    |φ_1(a)⟩ |φ_1(a)⟩ _1,1 =
    1/√(d)∑_m=0^d-1
    e^i2π[ ( 0.a_1a_2… a_q) +( 0.b_1) ]m   |m⟩



The effect of the second gate R_2^(d)  controlled by |b_2⟩   is to further evolve it (step[1,2]) in the state


    |φ_1(a)⟩ _1,1 |φ_1(a)⟩ _1,2 =
    1/√(d)∑_m=0^d-1
    e^i2π[ ( 0.a_1a_2… a_q) +( 0.b_1) + ( 0.0b_2) ]m   |m⟩



Proceeding in a similar way up to gate R_q^(d)  controlled by |b_q⟩, we find the final state (step[1,q]) of the first qudit which becomes


    |φ_1(a)⟩ _1,q-1 |φ_1(a)⟩ _1,q =
    1/√(d)∑_m=0^d-1
    e^i2π[ ( 0.a_1a_2… a_q) +( 0.b_1b_2… b_q)   ]m   |m⟩



In general, the final state of the l-th qudit of the lower register is found to be


    |φ_l(a)⟩ _l,q-l+1 =
    1/√(d)∑_m=0^d-1
    e^i2π[ ( 0.a_la_l+1… a_q) +( 0.b_lb_l+1… b_q)   ]m   |m⟩



Applying Eq. (<ref>) to each lower register qudits we can find that the  lower register has the final joint state 


    |φ(a)⟩ _1 |φ(a)⟩ _2⋯ |φ(a)⟩ _q⊗_l=1^q1/√(d)∑_d=0^d-1
    e^i2π[ ( 0.a_la_l+1… a_q) +( 0.b_lb_l+1… b_q) +  ]m   |m⟩ =
    |φ (a+b) ⟩



This is the quantum Fourier transform of the sum state |a+b  d^q⟩. By applying the inverse QFT at the lower register we can get the desired sum in the computational basis, while the upper register remains in the initial state |b ⟩. The required direct and inverse QFT blocks are not shown in Figure <ref>. 



 §.§ Adder of an integer with constant (ADDC_b)


Whenever one of the integers is constant, e.g. b=(b_1b_2… b_q), then the upper register in Figure <ref> is not necessary and all the controlled rotation gates become single qudit rotation gates  with their angles defined by the constant integer b. Thus (see Eqs. (<ref>) and (<ref>)), we must apply on the l-th qudit of the lower register  a sequence of q-l+1 rotation gates ( Φ_k^(d)) ^b_k+l-1 = ∑_m=0^d-1 e^i2π/d^kmb_k+l-1 |m⟩⟨ m|, for k=1… q-l+1. This product of gates can be merged into one gate of the form



    B_l(b) =    ∏_k=1^q-l+1( Φ_k^(d)) ^b_k+l-1 
    =
       ∏_k=1^q-l+1(  
    ∑_m=0^d-1 e^i2π m/d^k |m⟩⟨ m|
    ) ^b_k+l-1
    =
       ∑_m=0^d-1(  
    ∏_k=1^q-l+1
    e^i2π m/d^kb_k+l-1) |m⟩⟨ m|



These are diagonal gates of the form of Eq. (<ref>), and their angles depend on the constant b by the relation 


    φ_l,m(b)=
    ∑_k=1^q-l+12π/d^kmb_k+l-1


so they can be constructed with elementary R_z^(jk)(θ)  gates  using the procedure described in  subsection <ref>. Figure <ref> shows the constant b adder (direct and inverse QFT blocks not included in the diagram). Likewise the general adder ADD, this adder performs the addition modulo d^q.






 §.§ Single State Controlled Adder of an integer with constant (C_cADDC_b) 


The constant adder ADDC_b can be easily converted to a constant adder controlled by the state  of an additional control qudit so as to perform the transformation


    C_cADDC_b( |e⟩ |a⟩)=
    |e⟩ |a+bδ_ce⟩




where δ_ce is the Kronecker delta function. Consequently, the addition is performed iff the control state  equals |c⟩, otherwise the target state   |a⟩ remains unaltered. The one state controlled constant adder C_cADDC_b can be constructed as shown in Figure <ref> if the one qudit rotation gates B_l(b) of Figure <ref> are converted to the respective two qudits diagonal gates controlled by state |c⟩. These gates are exactly the CD_(c) gates of  subsection <ref>.






 §.§ Generalized Controlled Adder of an integer with constant (GCADDC_b) 


A useful generalization of the previous C_cADDC_b circuit can be achieved if we permit all the basis states of the control qudit to have an influence on the result of the addition. Such a circuit will be named Generalized Controlled Adder with constant b and is defined by  the relation


    GCADDC_b( |e⟩ |a⟩)=
    |e⟩ |a+be⟩





The above equation can be rewritten as 


    GCADDC_b( |e⟩ |a⟩)=
    |e⟩ |a+bδ _1e+2bδ _2e+ ⋯  +(d-1)bδ _(d-1)e⟩




Equation (<ref>) directly leads to the implementation of Figure <ref> where d-1 consecutive applications of C_cADDC_bc (c=1… d-1) adders are employed. 






 §.§ Multiplier with constant and Accumulator (MAC_b) 


A Multiplier with constant and Accumulator MAC_b multiplies a q qudits integer x with a constant b of q d-ary digits, and accumulates the product bx to a q qudits integer a (modulo d^q). Namely, the MAC_b circuit consists of two q qudits registers holding initially the states |x⟩  and |a⟩  and transforms them as 


    MAC_b( |x⟩ |a⟩)=
    |x⟩ |a+bx⟩



Taking into account that x can be written as (x_1x_2… x_q)=∑_l=1^qx_ld^q-l then Eq. (<ref>) can be written as 


    MAC_b( |x⟩ |a⟩)=    
    |x⟩ |a+b∑_l=1^qx_ld^q-l⟩ = 
        |x⟩ |a + x_qb + x_q-1db + ⋯ + x_1(d^q-1b) ⟩



This means that the above transformation can be implemented by applying q GCADDC circuits, where  the control is done consecutively by the qudits x_q,x_q-1,…,x_1 and the constant parameter for each one GCADDC block is b,db,…,d^q-1b (modulo d^q), respectively, as shown in Figure <ref>.






 §.§ Multiplier with constant  (MULC_b) 


A multiplier (modulo d^q) with constant b implements the function f:0… d^q-1→0… d^q-1 with y=f(x)=bx q^q. When cosntant b is relative prime to d^q then there exists the inverse b^-1d^q and consequently there exists the inverse function f^-1(y)=b^-1y d^q = b^-1bx d^q=x. This always happens when d is a prime number. 
Figure <ref> shows how to construct a Multiplier with constant b using two MAC_b blocks and the necessary direct and inverse QFT blocks. It requires a q qudits register initially holding the integer x and another q qubits ancilla register initially in zero state. At the end, one register is set to the state |bx d^q⟩ while the other register is set to state zero, so effectively the ancilla register is reset back and can be reused.





In the diagram of Figure <ref>, the boxes with the black strip at their right side are the "direct" blocks while these with the black strip at their left side are the respective inverses. The operation of the inverse MAC with parameter b^-1 is to perform substraction instead of accumulation, that is referring to Figure <ref>,  we have the operation MAC_b^-1^-1 | bx⟩ |φ(x) ⟩ = | bx⟩ |φ (x-b^-1(bx)) ⟩ = | bx⟩ |0⟩. 
The inverse MAC^-1 has the same internal topology as the direct MAC of Figure <ref> (of course with parameter b^-1 instead of b) with the only difference that the angles of its rotation gates have a minus sign. By inspecting the labels at the qudit buses of Figure <ref> describing the respective states we can conclude that the circuit implements the multiplication 


    MULC_b( |x⟩ |0⟩)= |bx⟩ |0⟩



Excluding the ancilla register, which is in the zero state before and after the operation and thus it remains unentangled, we conclude that this circuit performs the desired multiplication operation.





§ DIAGONAL OPERATORS ON Q QUDITS


The diagonal operator on q qudits of d levels, as its name implies, is a circuit whose unitary matrix of  dimensions d^q× d^q
has a diagonal form. The circuit developed in this section is such that the diagonal elements of the matrix are integer powers of the principal root of unity e^i2π/d^q and the integer powers are a function f(k) of the coordinate k=1… d^q-1 of the elements. In what follows, the diagonal operator circuit developed is for the function f(k)=γ k^2, where γ is an integer constant. Thus, the definition of our diagonal operator on q qudits is  


    Δ_γ^(q)=∑_k=0^Q-1 e^i2π/Qf(k) |k⟩⟨ k|


where Q=d^q and f(k)=γ k^2d^q. All the diagonal entries of the above matrix are integer powers of the basic phase ω =e^i2π/Q. The effect of this matrix upon a general superposition state of q qubits will be 


    ∑_k=0^Q-1 c_k |k ⟩∑_k=0^Q-1 c_k e^i2π/Qf(k)  |k ⟩



The circuit that implements the operator of Eq. (<ref>) will be derived by exploiting results of <cit.> which are given for the case of binary quantum circuits. A prerequisite for this construction is a Squarer/Multiplier with constant/Accumulator circuit (SMAC) that computes the function f  involving two q qudits registers as in 


    SMAC_γ(|k⟩ |z⟩ ) = |k⟩ |z + f(k) ⟩  =
    |k⟩ |z + γ k^2Q⟩



Such an SMAC circuit will be described in  subsection <ref>. 

Figure <ref> shows the diagonal operator circuit with entries dependable on the function f(k)=γ k^2Q. Two  quantum registers are used, each q qudits wide, namely Reg1 and Reg2. The upper register Reg1 is assumed to be in a general superposition state prior the operator Δ_γ^(q) is applied as described in Eq. (<ref>), while the lower register Reg2 is an ancilla register with zero initial and final state.

The first step is to form in the ancilla register Reg2 the uniform superposition state |R⟩ = 1/√(Q)∑_h=0^Q-1 |h⟩. This is accomplished with the application of q Hadamard gates H^(d) on each qudit of the register. 
Then, we apply on each qudit the diagonal gates D_Q^d^m† for m=0… q-1 . The matrix representing these gates is 





    D_Q^d^m†=
    diag(1,ω _Q ^-1· d^m,ω _Q ^-2· d^m,… ,ω _Q ^(d-1)· d^m ),      m=0… q-1


and it has exactly the same form of the diagonal gates of Eq. (<ref>). The joint affect of these gates at Reg2 is given by their tensor product which is a diagonal matrix too, of dimensions d^q× d^q 


    D_QQ^†=
    diag(1,ω _Q  ^-1,ω _Q ^-2, … ,ω _Q ^Q-2, ω _Q ^Q-1 )



Then, the state of Reg2  becomes


    |R⟩ _1= 1/√(Q) = ∑_h=0^Q-1ω ^(Q-h) |h ⟩



The initial state of Reg1 is assumed to be a general superposition of basis states and can be expressed as


    |S ⟩ _1= ∑_k=0^Q-1 c_k |k ⟩ = 
    ∑_n=0^Q-1∑_k ∈ K_n c_k |k ⟩ , 
         K_n={k: f(k)=n   }



In the right hand side of Eq. (<ref>) we have grouped all the basis states with value k such that f(k)=n in a set K_n and then sum over all the states belonging to sets K_n. The expediency of this grouping will be clear later.
Combining Eq. (<ref>) and (<ref>) we find the joint state of Reg1 and Reg2 just before the application of the SMAC block, which is given by the tensor product 


    |S ⟩ _1⊗ |R ⟩ _1 = 
    1/√(Q)∑_n=0^Q-1∑_k ∈ K_n∑_h=0^Q-1
    c_kω ^Q-h |k ⟩ |h ⟩



Taking into account the effect of the SMAC block given by Eq. (<ref>) we get the state of the two registers after the application of the SMAC 


    |SR ⟩ _2  =     SMAC(|S ⟩ _1⊗ |R ⟩ _1) = 
       1/√(Q)∑_n=0^Q-1∑_k ∈ K_n∑_h=0^Q-1
    c_kω ^Q-h |k ⟩ |h + f(k) Q ⟩
       1/√(Q)∑_n=0^Q-1∑_k ∈ K_n∑_h=0^Q-1
    c_kω ^Q-h |k ⟩ |h + n Q ⟩



We are going to use m=h+n Q as the index of the inner summation in place of h. We observe that for a particular n, as h takes the values from 0 to Q-1, then m=h+n Q takes one value a time ("1-1" mapping), that is the new index m will be in the same range from 0 to Q-1. Thus the lower and upper limits of the new index m remain the same and we have h=m-n Q and Q-h=n+(Q-m) Q. Also, it holds ω ^Q=1. Then Eq. (<ref>) becomes 


    |SR ⟩ _2  =    1/√(Q)∑_m=0^Q-1∑_n=0^Q-1∑_k ∈ K_n 
    c_kω ^nω ^Q-m |k ⟩ |m⟩
       1/√(Q)∑_n=0^Q-1∑_k ∈ K_n c_kω ^n |k ⟩⊗∑_m=0^Q-1ω ^Q-m   |m ⟩
       ∑_k=0^Q-1 c_kω ^f(k) |k ⟩⊗1/√(Q)∑_m=0^Q-1ω ^Q-m   |m ⟩
       ( Δ_γ^(q) |S ⟩) ⊗ |R ⟩



This shows that Reg1 has the desired state of Eq. (<ref>) and it is disentangled with respect to Reg2 which remains in state of Eq. (<ref>).  Thus, the ancilla Reg2 can be reset without any effect on the Reg1. The resetting can be accomplished as shown in Figure <ref> by applying in the reverse sequence (a) the inverse of the gates H^(d) and (b) the inverse of D_Q^d^m†,  which are the H^(d)†=H^(d)* (conjugate Hadamard) and D_Q^d^m, respectively. An alternative method would be to measure Reg2 and depending on the measurement result to apply GCX gates controlled by the measurement classical result. This measurement would not affect Reg1 as it is disentangled with respect to Reg2.



 §.§ Multiplier of two integers / Accumulator (MMAC)


The construction of the SMAC block requires a  multiplier of two integers and accumulator block (MMAC) whose operation is to multiply  integer x with integer y and accumulate the product xy to  integer z (modulo d^q). This means that the MMAC block is applied on three q qudits registers and performs the transformation


    MMAC ( |x⟩ |y⟩ |z⟩) = 
    |x⟩ |y⟩ |z+xy⟩




If x=(x_1x_2… x_q)=∑_t=1^q x_td^q-t and y=(y_1y_2… y_q)=∑_s=1^q y_sd^q-s are the d-base representations of the two integers , then their product (modulo d^q) is given by 


    xy=∑_s=0^q-1 d^s∑_t=0^q-1x_q-ty_q-s+t



In Eq. (<ref>) the full product terms corresponding to powers d^s with s ≥ q have not been included, because the product is to be calculated modulo d^q. Also, digits with negative index (e.g. x_-1), as well as with index greater than q (e.g. x_q+1), are assumed zero.
The calculation of the product and the accumulation can be performed in an similar way as in the MAC circuit given in subsection <ref>. We assume that the state corresponding to the accumulation register integer |z⟩ is already Fourier transformed and taking into account Eq. (<ref>) which expresses the QFT we expect that the l-th qudit of the accumulation result |z+xy⟩ prior the inverse QFT is 


    |φ_l(z+xy)⟩ = 
    |0⟩ +
    e^i2π/d^l (z+xy)   |1⟩ + 
    e^i2π/d^l (z+xy)2   |2⟩ + ⋯ +
    e^i2π/d^l (z+xy)(d-1)   |d-1⟩



Thus, to bring and initial state |φ_l(z)⟩  of the l-th qudit to the state of Eq. (<ref>)   we must  add various integer multiples of the basic angle (2π/d^l). Namely, taking into account Eq. (<ref>), the angles that must be added to the amplitude phases of a basis state |r⟩ (r=0… d-1) in the superposition  |φ _l(z)⟩ of Eq. (<ref>) are



    Φ_l,r = 2π/d^lxyr =
    2π r  ∑_s=0^q-1
    d^s-l∑_t=0^sx_q-ty_q-s+t=
    2π r ∑_s=0^l-1 d^s-l∑_t=0^s x_q-ty_q-s+t



The restriction s<l at the upper limit of the first sum of Eq. (<ref>) comes due to the periodicity exp(φ +2π d^n )= exp(φ ) that holds for any integer d and any non negative integer n. The restriction t≤ s at the upper limit of the second sum results because y_q-s+t=0 for t>s. Replacing with k=l-s, Eq. (<ref>) becomes 


    Φ_l,r = 
    2π r ∑_k=1^l d^-k∑_t=0^l-k x_q-ty_q+k-l+t



Consequently, the angles that must be added to the phase amplitude of the |r⟩ component of the superposition are (2π/d^k)x_my_nr and depend on indices m=q-t and n=q+k-l+t. This can be attained if we introduce the notion of a double controlled generalized rotation gate applied to three qudits, two controls and one target.
Similarly to Eq. (<ref>) which is the definition of the simply controlled generalized rotation gate, we define the double controlled generalized rotation gate R _k^(d) with the d^3× d^3 matrix


    R _k^(d) =
    ∑_m=0^d-1 ∑_n=0^d-1 ∑_r=0^d-1 
    e^i2π/d^kmnr( |m⟩⟨ m| ) ⊗( |n⟩⟨ n| ) ⊗( |r⟩⟨ r| )


Figure <ref> depicts the symbol for this double controlled rotation gate. In Appendix A a construction of R _k^(d)   will be presented using some of the elementary and basic gates introduced in Section 3.



The topology of the MMAC circuit can be directly concluded from Eq.(<ref>), as this equation describes  which  gates have to applied and which are their control connections to the qudits carrying |x⟩ and |y⟩. Figure <ref>  shows an example MMAC for the case of q=4. In this figure the R _k^(d)  gates are represented with the value k inside the circle. Generalization for any value of q is obvious. 




We observe in Figure <ref> and in Eq. (<ref>) that l-k+1 R _k^(d) gates are sequentially applied on the l-th target qudit for a specific k (l=1 … q,  k=1 … l). In total, C_l=∑_k=1^l (l-k+1) =l(l+1)/2 R _k^(d)  gates are applied on the l-th target qudit. Summing over all target qudits we find the total number of gates used  C_MMAC(q)=∑_l=1^q C_l = ∑_l=1^ql(l+1)/2 =1/6q^3+1/2q^2 + 1/3q. The same value gives the depth of the circuit as arranged in Figure <ref>. Indeed, for the example q=4 we find C_MMAC(4)=20.  
We can exploit the fact that gates R _k^(d)  mutually commute as they are diagonals and rearrange them so as to achieve a parallelization in their execution. The gates that can be executed simultaneously are those that operate on different qudits. An example of the proposed parallelization for the case q=4 is shown Figure <ref>, where below each gate is shown the soonest timestep in which it can be executed. E.g. at the first timestep three gates can be executed in parallel as none of these gates operate on the same qudit as the other two. We can generalize this parallelization scheme and conclude that we can achieve a depth of about q(q+1)/2 which is quadratic instead of cubic without the proposed rearrangement.




 §.§ Squarer/Multiplier/Accumulator  (SMAC)


The MMAC circuit allows the  construction of the SMAC_γ circuit described by Eq. (<ref>) and required for the q qudits diagonal operator Δ_γ^(q). The Squarer/Multiplier with constant γ /Accumulator modulo d^q is presented in block diagram in Figure <ref>. It uses 4q qudits, 2q of which are ancilla qudits with zero initial and final state, where q is the number of qudits used to represent the argument x. The 4q qudits are grouped into four registers of q qudits each. The second register from top holds the argument |x⟩ while the bottom register holds the accumulation value |z+γ x^2 d^q⟩. 



The first step is to set the state of the top register into the same state as the second one, which is |x⟩. This is accomplished with the adder block sandwiched between two QFT blocks, direct and inverse (This operation could be also achieved using a sequence of GCX gates to "copy" the second's register state to the first). In a second step, the two states |x⟩ and |x⟩  of the two top registers are multiplied together by the MMAC and  the product is accumulated to the third register from top, which was initially in the zero state. At this stage, the joint state of the three top registers is |x⟩ |x⟩ |x^2d^q⟩. Next, the MAC_γ block follows to multiply the constant γ with the |x^2d^q⟩ state of the third register. The result is accumulated to the bottom register, which was initially in state |z⟩. At this point the joint state of the four registers can be described by |x⟩ |x⟩ |x^2d^q⟩ |z+γ x^2d^q⟩.
What remains is to reset the first and the third ancilla registers. The inverse MMAC resets the third register by performing substraction instead of accumulation of the product |x^2⟩. The inverse MMAC is constructed like the direct MMAC with opposite angles in its rotation gates. Last, the inverse adder resets the top ancilla register. Consequently, the circuit of Figure <ref> implements the transformation 


    SMAC ( |0⟩ |x⟩ |0⟩ |z⟩) =
    |0⟩ |x⟩ |0⟩ |z + γ x^2d^q⟩

which is exactly the transformation of Eq. (<ref>) if the ancilla registers are ignored.




§ COMPLEXITY ANALYSIS 




The arithmetic quantum circuits proposed in the previous sections are broken down to the level of elementary gates H^(d), R_z^(jk)(θ), R_x^(jk)(θ) and GCX_m^(jk) introduced in Section <ref>. This decomposition is depicted in Figure <ref> in a tree structure, where the root of of each tree is some of the complete circuits proposed and the leaves of the tree (trapezoids) represent the elementary gates. The edges of each tree are labeled with the number of components needed by each level from one level below (no label stands for 1). The SMAC and the Diagonal operator are not included in this Figure, but their costs and depths can be easily calculated after the calculations of the blocks shown in this figure.


A rough complexity analysis in terms of quantum cost (number of elementary gates used) and depth (execution time) can be done with the help of Figure <ref>. The analysis assumes that single and two qudits gates are equivalent in terms of costs and execution time. Exact costs and depths depend on the particular implementations. The total gates count for each block can be found by traversing the tree emerging from the inspected block down to each leaf of the subtree. The labels of the edges for each path are multiplied and then the products of each path used are summed together. E.g. the QFT circuit needs q Hadamard gates, (q^2-q/2)(d-1)(2d-1) GCX_m^(jk)gates and (q^2-q/2)(d-1)2d R_z^(jk)(θ)  gates. Similar calculations provide us with the quantum costs shown in Table <ref>, which shows only the highest order terms. 



The depth calculation will be done in more detail by finding first the depths of QFT, ADD, MAC and MMAC blocks and then the depths of MULC and SMAC. 



  QFT At first glance Figure <ref> exhibits a quadratic depth O(q^2), but it can be easily shown that we can parallelize the execution with an appropriate reordering of the gates and thus achieve a linear depth, namely depth(QFT)=8d^2q.


  ADD Similarly as in the QFT case, a reordering of gates in Figure <ref> offers a linear depth too, that is depth(ADD)=4d^2q.


  MAC Concurrent execution of gates is possible in this case, too. It can be easily seen that by flattening the hierarchy MAC-GCADDC-CADDC, q different controlled gates B_l(b) (Eq. (<ref>)) belonging in different GCADDC blocks can be executed concurrently. Thus, the depth of the MAC is of the order O(4d^2q) instead of O(4d^2q^2) as directly calculated by the number of elementary gates.


  MULC Observing Figure <ref> we find depth(MULC)=3depth(QFT)+2depth(MAC), as the two middle QFT blocks (direct and inverse) can be executed simultaneously. Thus, we derive depth(MULC)= 32d^2q.


  MMAC The reordering of gates achieves q(q+1)/2 execution steps of double controlled rotation gates. Taking into account the decomposition of these three qudit gates into single and two qudits gates (see Appendix A) we end up in depth(MMAC)=21d^3q^2. 


  SMAC From the previous calculations and Figure <ref> we find that the dominant depth of SMAC in leading order is twice the depth of the MMAC block.


  Δ _γ^q The depth of the diagonal circuit is essentially the depth of the SMAC.








§ CONCLUSIONS AND FUTURE WORK

In this paper we presented an assortment of quantum circuits for multilevel qudits. These are basic integer arithmetic operations circuits (like addition, multiplication/accumulation and multiplication) as well as more complex circuits such as squarers. Additional extensions can be applied. E.g., the ADD, ADDC, MAC and MULC circuits can be converted to single qudit controlled versions. Such controlled versions could be useful to the multilevel qudits quantum phase estimation algorithms and quantum simulations. 

The general diagonal operator has been developed for the special case of a quadratic function f(x)=γ k^2, where k is the coordinate of the diagonal element, however using the same techniques we can easily construct diagonal operators for any power of k and even for a polynomial function on k. E.g. the Squarer/Multiplier/Accumulator can be converted to a circuit that accumulates  the third power by inserting additional MMAC units in Figure <ref>.

The designs  are based on the alternative representation of an integer after QFT transformed instead of the usual computational basis representation, a method which has been already exploited in the binary qubits case. QFT based arithmetic circuits design is a versatile method to develop many arithmetic circuits. E.g. there is no need to handle carries which leads to space reduction. Moreover, if it is suitably used, it can offer advantages in terms of speed. This is possible when similar blocks are iterated to act on a datapath whose state follows the QFT representation. The extensive usage of rotation gates (which mutually commute) on such a datapath permits their rearrangement so as they execute concurrently. This capability is observed in the MAC block, where the application of a suitable reordering of gates led to depth reduction from O(q^2) to O(q). Similarly, the depth of the MMAC block reduced from O(q^3) to O(q^2).

Another advantage that has been observed in designs adopting the QFT method is their robustness to various kinds of deviations from the ideal operation. E.g. approximate QFT <cit.> or QFT banding is the design procedure of eliminating small angle rotation gates. Studies of the Shor's algorithm which uses the QFT showed that the algorithm still works sufficiently even when a large proportion of the QFT rotation gates are eliminated <cit.>. Recent studies extended to circuits beyond QFT. In <cit.> the simultaneous gate pruning of rotation gates of the QFT circuit and the QFT based modular exponentiator of Beauregard's circuit <cit.> were simulated. The simulation results showed similar robustness of Shor's algorithm to these gates eliminations. This robustness is sustained even if the parameters of the remaining rotation gates are randomly selected <cit.>. The above results suggest that a similar robustness is expected in the multidimensional qudits case and further investigation to be carried.

On the other side, there is a drawback related to the requirement of reliable implementing high accuracy small angles rotation gates. Moreover, these gates must belong to a set of fault tolerant gates if large scale quantum computation is considered. Fortunately, as shown in Appendix B, approximation of these gates is possible, albeit with a cost. However, the remarks of the previous paragraph suggest that this cost may be much lower if approximate computation is adopted.

For the above reasons and also because the exact cost depends on the exact technology used, which for qudits is at an early stage, the complexity analysis of section <ref> is to be considered as a crude indicator of performance. Despite that, we think that the proposed designs enrich the toolkit of the future quantum computing.




References
000
Mut:2000 A.Muthukrishnan and C.R.Stroud, Multivalued logic gates for quantum computation, Phys. Rev. A, Vol.62, Iss.5, pp.052309, (2000).  

Mut:2002 A.Muthukrishnan and C.R.Stroud, Quantum fast Fourier transform using multilevel atoms, J. Modern Optics, vol.49, Iss.13, pp.2115-2127, (2002).

Lan:2009 B.P.Lanyon, M.Barbieri, M.P.Almeida, T.Jennewein, T.C.Ralph, K.J.Resch, J.G.Pryde, J.L.O'Brien, A.Gilchrist and A.G.White, Simplifying quantum logic using higher-dimensional Hilbert spaces, Nature Phys., Vol.5, Iss.2, pp. 134-140 (2009).

Got:1999 D.Gottesman, Fault-Tolerant Quantum Computation Higher-Dimensional Systems, Chaos, Solitons & Fractals, Vol.10, Iss.10, pp.1749-1758, (1999).

Ket:2006  A. Ketkar, A. Klappenecker, S.Kumar and P.K. Sarvepalli, Nonbinary Stabilizer Codes Over Finite Fields, IEEE Trans. Info. Theory, Vol.52, Iss.11, pp.4892-4914, (2006).

Cam:2014 E.T.Campbell, Enhanced Fault-Tolerant Quantum Computing in d-Level Systems, Phys. Rev. Lett., Vol.113, Iss.23, pp.230501, (2014).

Shor:1994	
P.W.Shor, Algorithms for quantum computation: Discrete log and factoring, Proc. 35th IEEE Ann. Symp. Foundations  Comp. Sci. (FOCS), pp. 124-134, (1994).

Asp:2005 A.Aspuru-Guzik, A.D.Dutoi, P.J.Love and M.Head-Gordon, Simulated Quantum Computation of Molecular Energies, Science, Vol. 309, Iss.5741, pp.1704-1707, (2005). 

Par:2011 V.Parasa and M.Perkowski, Quantum Phase Estimation Using Multivalued Logic, Proc. 41st IEEE Int. Symp. Multiple-Valued Logic (ISMVL), pp.224-229, (2011).

Fan:2007 Y.Fan A Generalization of the Deutsch-Jozsa Algorithm to Multi-Valued Quantum Logic, Proc. 37th IEEE Int. Symp.  Multiple-Valued Logic (ISMVL), pp.12, (2007).

Kik:2015
E.O.Kiktenko, A.K.Fedorov, A.A.Strakhov and V.I.Man'ko, Single qudit realization of the Deutsch algorithm using superconducting many-level quantum circuits, Phys. Lett. A, Vol.379, Iss.22-23, pp.1409-1413, (2015).

Grover:1996 L.K.Grover, A Fast Quantum Mechanical Algorithm for Database Search, Proc. 28th ACM Symp. Theory Comp. (STOC), pp. 212-219, (1996).

Iva:2012
S.S. Ivanov, H.S.Tonchev, and N.V.Vitanov, Time-efficient implementation of quantum search with qudits, Phys. Rev. A, Vol.85, Iss.6, pp.062321, (2012).

Di:2013 Y.-M. Di and H.-R.Wei, Synthesis of Multivalued Quantum Logic Circuits by Elementary Gates, Phys. Rev. A,  Vol.87, Iss.1, pp.012325, (2013).

Nee:2009 M.Neeley, M.Ansmann, R.C.Bialczak, M.Hofheinz, E.Lucero, A.D.O'Connell, D.Sank, H.Wang, J.Wenner, A.N.Cleland, M.R.Geller, J.M.Martinis, Emulation of a Quantum Spin with a Superconducting Phase Qudit, Science, Vol.325, Iss.5941, pp. 722-725, (2009).

Str:2011 F.W.Strauch, Quantum logic gates for superconducting resonator qudits, Phys. Rev. A, Vol. 84, Iss. 5, pp. 052313, (2011).

Bab:2017 A.Babazadeh, M.Erhard, F.Wang, M.Malik, R.Nouroozi, M.Krenn and A.ZeilingerHigh-Dimensional Single-Photon Quantum Gates: Concepts and Experiments, arXiv e-print quant-ph/1702.07299, (2017).

Mal:2016 M.Malik, M.Erhard, M.Huber, M.Krenn, R.Fickler and A.Zeilinger, Multi-photon entanglement in high dimensions, Nature Photonics, Vol.10, Iss.4, pp.248–252, (2016).

Dra:1998
T.G.Draper, Addition on a Quantum Computer, arxiv e-print quant-ph/0008033, (1998).

Bea:2003
S.Beauregard, Circuit for Shor's algorithm using 2n+3 qubits, Quantum Info. & Comput., Vol.3, Iss.2 , pp 175-185, (2003).

Pav:2014	A.Pavlidis and D.Gizopoulos, Fast quantum modular exponentiation architecture for Shor’s factoring algorithm, Quantum Info. & Comput., Vol.14, no.7&8, pp.0649-0682, (2014).

Bar:1996 A.Barenco, A.Ekert, K.-A.Suominen and P.Törmä, Approximate quantum Fourier transform and decoherence, Phys. Rev. A, Vol.54, Iss.1, pp.139-146, (1996).

Nam:2015 Y.S.Nam and R.Blümel, Analytical formulas for the performance scaling of quantum processors with a large number of defective gates, Phys. Rev. A, Vol.92, Iss. 4, pp. 042301, (2015).

Flo:2017 E.Floratos and A.Pavlidis, Quantum Fractional Fourier Transform Implementation on d-level qudits, to appear in arxiv/quant-ph.

Bry:2002 J.-L. Brylinski and R. Brylinski, Universal Quantum Gates, Mathematics of Quantum
Computation, edited by R. Brylinski and G. Chen, CRC
Press, (2002).

Kha:2006 F.S.Khan and M.Perkowski, Synthesis of multi-qudit hybrid and d-valued quantum logic circuits by decomposition, Theor. Comp. Sci., Vol.367, Iss.3, pp.336-346, (2006).

Bre:2005 G.K.Brennen, D.P.O'Leary  and S.S.Bullock, Criteria for exact qudit universality, Phys. Rev. A, Vol.71, Iss.5, pp.052318, (2005).

Bul:2005 S.S.Bullock, D.P.O'Leary  and G.K. Brennen
Asymptotically Optimal Quantum Circuits for d-Level Systems, Phys. Rev. Lett., Vol.94, Iss.23, pp.230502, (2005).

Per:2004 M.Perkowski and E.Curtis, A Transformation Based Algorithm for Ternary Reversible Logic Synthesis using Universally Controlled Ternary Gates, Proc. 30th Int. Workshop Logic Synth. (IWLS), pp. 345-352, (2004).

Mil:2004 D.M.Miller, G.W.Dueck and  D.Maslov, A Synthesis Method for MVL Reversible Logic, Proc. 34th Int. Symp. Multiple-Valued Logic (ISMVL), pp.74-80, (2004).

Den:2004 N.Denler, B.Yen, M.Perkowski and P.Kerntopf, Synthesis of reversible circuits from a subset of Muthukrishnan–Stroud quantum realizable multi-valued gates, Proc. 30th Int. Workshop Logic Synth. (IWLS), pp.321–328, (2004).




Ved:1996
V.Vedral, A.Barenco, A.Ekert, Quantum networks for elementary arithmetic operations, Phys. Rev. A, Vol. 54, Iss. 1, pp. 147-153,(1996).

Cuc:2005
S.A.Cuccaro, T.G.Draper, S.A.Kutin, D.P.Moulton, A new quantum ripple carry addition circuit, QIP workshop MIT, arxiv e-print quant-ph/0410184v1 (2005).

VanM:2005
R.Van Meter, K.M.Itoh, Fast quantum modular exponentiation, Phys. Rev. A, Vol. 71, Iss.5, 052320, (2005).

Dra:2006
T.G.Draper, S.A.Kutin, E.M.Rains, K.M.Svore, A logarithmic-depth quantum carry-lookahead adder, Quantum Info. & Comput., Vol. 6, Iss.4, pp.351-369, (2006).

Kho:2011
A.Khosropour, H.Aghababa, B.Forouzandeh, Quantum Division Circuit Based on Restoring Division Algorithm, Proc. 8th Int. Conf. Info. Tech. New Generations (ITNG), pp. 1037-1040, (2011).

Cho:2012
B.-S.Choi, R.Van Meter, Θ(√(n))-depth Quantum Adder on a 2D NTC Quantum Computer Architecture, ACM J. Emerg. Tech. Comput. Syst., Vol.8, Iss.3, Art. 24 (2012).


Kha:2007 M.H.Khan and M.A.Perkowski, Quantum Ternary Parallel Adder/Subtractor with Partially-look-ahead Carry, J. Syst. Archit., Vol.53, Iss.7, pp.454-464, (2007).

Sat:2007 T.Satoh, S.Nagayama and  R. Van Meter, A Reversible Ternary Adder for Quantum Computation, 17th Asian Conf. Quantum Info. Sci. (AQIS), (2007).

Kha:2008 M.H.A.Khan, Synthesis of quaternary reversible/quantum comparators, J. Syst. Archit., Vol.54, Iss.10, pp.977-982, (2008).

Boch:2016  A.Bocharov, S.X.Cui, M.Roetteler and K.M.Svore, Improved Quantum Ternary Arithmetics, Quantum Info. & Comput., Vol.16, Iss.9&10, pp.862-884, (2016).
 
 

Zal:1998 C.Zalka,  Simulating quantum systems on a quantum computer, Proc. Roy. Soc., Vol. 454, Iss.1969, pp.313-322, (1998).

Kas:2008 I.Kassal,  S.P.Jordan,  P.J.Love, M.Mohseni and A.Aspuru-Guzik  Polynomial-time quantum algorithm for the simulation of chemical dynamics, Proc. Natl. Academy Sci., Vol.105, Iss.48, pp.18681-18686, (2008). 

Wel:2014 J.Welch, D.Greenbaum, S.Mostame and A.Aspuru-Guzik, Efficient Quantum Circuits for Diagonal Unitaries Without Ancillas, New J. Physics, Vol.16, Iss.3, pp.033-040, (2014).

Hog:1999 T.Hogg, C.Mochon, W.Polak and E.Rieffel, Tools for Quantum Algorithms, Int. J. Mod. Phys. C, Vol.10, Iss.7, pp.1347-1362, (1999).

Wel:2016 J.Welch, A.Bocharov,and K.M.Svore, Efficient Approximation of Diagonal Unitaries over the Clifford+T Basis, Quantum Info. & Comput., Vol.16, Iss. 1&2, pp.87-104, (2016).


Bul:2004 S.S.Bullock and I.L.Markov, Asymptotically optimal circuits for arbitrary n-qubit diagonal computations, Quantum Info. & Comput. Vol.4, Iss. 1, pp.27-47, (2004).


Bee:2016 K.Beer and F.A.Dziemba, Phase-context decomposition of diagonal unitaries for higher-dimensional systems, Phys. Rev. A, Vol.93, Iss.5, pp.052333, (2016).


Erm:2007 A.S.Ermilov and V.E.Zobov, Representation of the Quantum Fourier Transform on Multilevel Basic Elements by a Sequence of Selective Rotation Operators, Optics and Spectroscopy, Vol. 103, Iss.6, pp.969-975, (2007).



Coppersmith:1994 D.Coppersmith, An approximate Fourier transform useful in quantum factoring, IBM Research Div., T.J.Watson Research Center, Tech. Rep. RC19642, (1994).


Fowler:2004 A.G.Fowler and L.C.L.Hollenberg, Scalability of Shor's algorithm with a limited set of rotation gates, Phys. Rev. A, Vol. 70, Iss.3 pp. 032329, (2004).

Nam:2012 Y.S.Nam and R.Blümel, Performance scaling of Shor's algorithm with a banded quantum Fourier transform, Phys. Rev A, Vol.86, Iss.4, pp.044303, (2012).


Nam:2013a Y.S.Nam and R.Blümel, Robustness and performance scaling of a quantum computer with respect to a class of static defects, Phys. Rev. A, Vol. 88, Iss. 6, pp.062310, (2013).

Nam:2013b Y.S.Nam and R.Blümel, Streamlining Shor's algorithm for potential hardware savings, Phys. Rev. A, Vol.87, Iss.6, pp. 060304, (2013).



Nie:2010	M.A.Nielsen and I.L.Chuang, Quantum Computation and Quantum Information, Cambridge Press, (2010).


Har:2002 A.W.Harrow, B.Recht and I.L.Chuang, Efficient Discrete Approximations of Quantum Gates, J. Math. Phys., Vol.43, Iss.9, pp. 4445, (2002).

Fow:2011 A.G.Fowler, Constructing arbitrary Steane code single logical qubit fault-tolerant gates, Quantum Info. & Comput., Vol.11, Iss.9&10, pp.867-873, (2011).

Cod:2012	N.Cody Jones, J.D.Whitfield, P.L.McMahon, M.-H. Yung, R.Van Meter and A.Aspuru-Guzik and Y.Yamamoto, Simulating chemistry efficiently on fault-tolerant quantum computers, New J. of Physics, Vol. 14, Iss.11, pp.115023, (2012). 

Boc:2012	A.Bocharov and K.M.Svore, A Depth-Optimal Canonical Form for Single-qubit Quantum Circuits, Phys. Rev. Lett., Vol. 109, Iss.19, pp.190501, (2012).

Pha:2013	T.T.Pham, R.Van Meter and C.Horsman, Optimization of the Solovay-Kitaev algorithm, Phys. Rev. A, Vol.87, Iss.5, pp.052332, (2013).

Duc:2013	G.Duclos-Cianci and K.M.Svore, A State Distillation Protocol to Implement Arbitrary Single-qubit Rotations, Phys. Rev. A, Vol. 88, Iss.4, pp.042325, (2013).

Kli:2013a	V.Kliuchnikov, D.Maslov and M.Mosca, Asymptotically optimal approximation of single qubit unitaries by Clifford and T circuits using a constant number of ancillary qubits, Phys. Rev. Lett., Vol. 110, Iss.19, pp.190502, (2013).

Boc:2013	A.Bocharov, Y.Gurevich and K.M.Svore, Efficient Decomposition of Single-Qubit Gates into V Basis Circuits, Phys. Rev. A, Vol. 88, Iss.1, pp.012313, (2013). 

Kli:2013b	V.Kliuchnikov, D.Maslov and M.Mosca, Fast and efficient exact synthesis of single qubit unitaries generated by Clifford and T gates, Quantum Info. & Comput., Vol.13, No7&8, pp. 607-630, (2013).

Sel:2015	P.Selinger, Efficient Clifford+T approximation of single-qubit operators, Quantum Info. & Comput., Vol.15, Iss.1&2, pp. 159-180, (2015). 

Ros:2015	N.J.Ross and P.Selinger, Optimal ancilla-free Clifford+T approximation of z-rotations, Quantum Info. & Comput., Vol.15, Iss.11&12, pp.932-950, (2015).

Kli:2016	V.Kliuchnikov, D.Maslov and M.Mosca, Practical approximation of single-qubit unitaries by single-qubit quantum Clifford and T circuits, IEEE Trans. On Computers, Vol.65, Iss.1, pp.161-172, (2016). 

Kit:1997 A.Y.Kitaev, Quantum computations: Algorithms and error correction, Russ. Math. Surv., Vol. 52, Iss.6, pp.1191-1249, (1997).

Kit:1999 A.Y.Kitaev, A.H.Shen, and M.N.Vyalyi,  Classical and Quantum Computation, AMS Grad. Stud. Math., Vol.47, (1999).

Daw:2006 C.M.Dawson and M.A.Nielsen, The Solovay-Kitaev algorithm, Quantum Info. & Comput., Vol.6, Iss.1, pp. 81-95, (2006). 






The double controlled gates R _k^(d) which are used in the MMAC block can be decomposed to single and two qudits elementary gates. Equation (<ref>) states that the effect of the R _k^(d)  to the target qudit state is  multiplication by a diagonal matrix d × d of the form diag(1,e^i2π/d^kmn,e^i2π/d^k2mn, … ,    e^i2π/d^k(d-1)mn) iff the two control states are  |m⟩ and |n⟩. Consequently, to implement this gate we need (as was the case of the  gate CD_m(φ _1, … , φ _d-1 )), the construction of double controlled diagonal gates of the form 


    CCD_(m,n)( φ_1,φ_2,…,φ_d-1)=
    diag(I_d, … ,I_d,  mn-th    blockD(φ_1,φ_2 ,…,φ_d-1) ,I_d, … ,I_d  )


where the diagonal D(φ_1,φ_2 ,…,φ_d-1) is applied to the target qudit iff the two control qudits states are |m⟩ and |n⟩. The angles are φ _l =(2π / d^kmnl) for l=1 … d-1. Thus, the gate R _k^(d)   is constructed by successively using (d-1)× (d-1) three qudit gates CCD_(m,n)(φ_1,φ_2,…,φ_d-1) as follows 


    R _k^(d) =
    ∏_m=0^d-1∏_n=0^d-1 CCD_(m,n) (1,i2π/d^kmn,i2π/d^k2mn, … ,    i2π/d^k(d-1)mn )





The above decomposition is depicted in Figure <ref>. The parameter k inside the rectangular symbol of the CCD gates corresponds to the parameter k of R _k^(d) gate, while the values m and n inside the small circles of the same gate signify that the (mn)-th block of the diagonal matrix CCD is of the form diag(1,e^i2π/d^kmn,e^i2π/d^k2mn, … ,    e^i2π/d^k(d-1)mn) while the rest of the blocks are identity matrices (see Eq. (<ref>)). That is the gate transforms the target qudit with the matrix diag(1,e^i2π/d^kmn,e^i2π/d^k2mn, … ,    e^i2π/d^k(d-1)mn) iff the control states are |x⟩ = |m⟩ and |y⟩ = |n⟩.


  
The way to construct a double controlled rotation gate CCD_(m,n)(φ_1,φ_2,…,φ_d-1) is analogous to the one for the  simply controlled gate  CD_(m)(φ_1,φ_2,…,φ_d-1) which is equivalent to the CD^'_m(a_1,a_2,…,a_d-1). The difference in this case is that  we need  double controlled generalized NOT gates, which will be called GCCX_(m,n)^(jk). They can be thought as an extension of Toffoli gates to the qudit case and their operation is analogous to that of  the GCX_(m)^(jk), but in three qudits. They are defined by the equation 


    [                                                                     GCCX_(m,n)^(jk)= |m ⟩⟨ m| ⊗ 
    |n ⟩⟨ n| ⊗( |j ⟩⟨ k| + |k ⟩⟨ j| +
     ∑_l
     k=0 
    
     k ≠ j 
     ^d-1 |k ⟩⟨ k| ) +                                                                                     ;                                                                                                                      ∑_l
     l=0 
    
     l ≠ m 
    
     l ≠ n
     ^d-1 |l ⟩⟨ l| ⊗ I_d⊗ I_d                                                                      j,k,m,n=0 … d-1 ]



This description means that they interchange the two target qudit basis states |j⟩ and |k⟩ iff the two control states are |m⟩ and |n⟩. Having available the GCCX_(m,n)^(jk) gates we can construct a CCD^'_(m,n)(a_1,a_2,…,a_d-1) as shown in Figure <ref> (the example depicts the qutrit case of d=3, generalization for other values of d is obvious). The controlled gate S_n is the analogous of the single qudit gate S_m of Figure <ref>. Namely, it is a CD_m(1,…,1, n-th pose^iφ,1,…,1) gate, with angle φ=1/d∑_k=1^d-1a_k. The CCD_(m,n)^' gate is identical with the desired CCD_(m,n) gate if angles redefinition similar to the ones of Eq. (<ref>) are applied.





What remains is the GCCX_(m,n)^(jk) gate  construction. This gate operates in a two dimensional subspace of the target qudit. Thus, extension to the qudit case of a Toffoli gate decomposition into single and two qubits <cit.> can be exploited and this is shown in Figure <ref>. The single and two qudit gates used will be the generalization of the S,T and H qubit gates to the d dimension of the qudits but operating only on a 2-dimensional subspace. Concretely, we define the gates 


    [ S^(jk)=|j ⟩⟨ j| + i|k ⟩⟨ k| +
    ∑_l
    n=0 
    
    n ≠ j 
    
    n ≠ k
    ^d-1 |n ⟩⟨ n|                                                         j,k=0 … d-1 ]




    [ T^(jk)=|j ⟩⟨ j| + e^iπ /4|k ⟩⟨ k| +
    ∑_l
    n=0 
    
    n ≠ j 
    
    n ≠ k
    ^d-1 |n ⟩⟨ n|                                                               j,k=0 … d-1 ]




    [ H^(jk)=1/√(2)( |j ⟩⟨ j| + |j ⟩⟨ k| + 
    |k ⟩⟨ j| - |k ⟩⟨ k| ) +
    ∑_l
    n=0 
    
    n ≠ j 
    
    n ≠ k
    ^d-1 |n ⟩⟨ n|                                                                                         j,k=0 … d-1 ]






Gates S^(jk) and T^(jk) are effectively R_z^(jk)(θ) with θ equal to π /2 and π /4, respectively, ignoring some global phase. Also, we can build the H^(jk) gate using the easily proved identity H^(jk)=e^iπ /2 R_z^(jk)(π /2) R_x^(jk)(π /2) R_z^(jk)(π /2), which is similar to the one for the qubit Hadamard case.
We have finally achieved to synthesize a three qudits rotation gate R _k^(d) with elementary single qudit rotation R_z^(jk)(θ), R_x ^(jk)(θ) gates and two qudits GCX_(m)^(jk) gates.



The design of arithmetic quantum circuits based on the QFT involves a library of elementary gates  { H^(d),  R_z^(jk)(θ), R_x^(jk)(θ), GCX_m^(jk)}. The size of this library is not constant as the parameter θ of the rotation gates depends on the size q of the circuit (number of qudits used). However, it is possible to approximate these rotation gates with arbitrary precision using a constant set of gates.  Much research has been done recently in this area focused on gates operating on qubits <cit.>. An extension of some of these results can be easily applied to the case of the qudits for the specific elementary gates library and thus we can use a constant library for the synthesis. This is important both for physical implementation reasons and for the fault tolerance aspect of the circuit, as fault tolerance techniques have developed for a restricted set of gates (e.g Clifford + T gates for the case of binary qubits). 
Before proceeding to the extension of the well established qubit gates approximation methods to the qudit case, some definitions are necessary. A unitary matrix U of dimensions d × d is called a two level matrix if it has the form  <cit.>


    U=[    1                                   ;         ⋱                              ;              1                         ;                u_jj      u_jk          ;                        1               ;                u_kj      u_kk          ;                                  ⋱     ;                                       1;      ]



This kind of matrix leaves invariant a subspace of d-2 dimensions and operates only on the two dimensions corresponding to coordinates j and k. A more compact notation for the above matrix is 


    [                U_[ jk ](V), V=
    [ v_11 v_12; v_21 v_22 ] ]



In this notation, we only need to define a unitary 2 × 2 dimensional matrix V and also define which coordinates j and k this matrix operates on. A multiplication homomorphism is valid since


    U_[ jk ](V_1)· 
    U_[ jk ](V_2) =
    U_[ jk ](V_1· V_2)



The R_z^(jk)(θ) and R_x^(jk)(θ) elementary gates have exactly the form of Eq. (<ref>) with V_z(θ)=[ e^- iθ/2        0;        0  e^ iθ/2 ]
and 
V_x(θ)=[   cosθ /2 -isinθ /2; -isinθ /2   cosθ /2 ]
, respectively.  The V_z(θ) and V_x(θ) gates are in fact qubits rotation gates  and thus we can  exploit the known approximation results of the literature for the qubit gates. These results state that an arbitrary rotation gate like V_z(θ) can be approximated by a finite sequence of gates belonging to a discrete set, e.g.  V̂_z(θ)=(HT⋯ T)(HT⋯ T)⋯ (HT⋯ T), where H is the Hadamard gate and T is the π/8 gate, such as the approximation error ϵ = V̂_z(θ) - V_z(θ) can be arbitrary small (Solovay-Kitaev Theorem and improvements <cit.>). Using this fact and Eq. (<ref>) we find that every rotation gate R_z^(jk)(θ) can be approximated by another one R̂_z^(jk)(θ) with arbitrary precision as 


    R̂_z^(jk)(θ)=
    U_[jk ](V̂_z(θ) )=
    U_[jk ](H) U_[jk ](T) ⋯ 
    U_[jk ](T) 
    ⋯
    U_[jk ](H) U_[jk ](T) ⋯ 
    U_[jk ](T)



In essence, the U_[jk ](H) and U_[jk ](T) gates are the H^(jk) and T^(jk) gates of Eq. (<ref>) and (<ref>), respectively.
On the other side, the R_x^(jk)(θ) gates can be decomposed using the identity R_x^(jk)(θ)=H^(jk) R_z^(jk)(θ) H^(jk) and thus the proposed circuits can be synthesized using the discrete library of constant number of components shown in Table <ref>. The second column shows the number of different gates of the same family, which depends on the family parameters (none, j, k and m). The constant library consists of a total of (3+d)d(d-1)/2+1 gates.



The first Solovay-Kitaev algorithms <cit.> generate a sequence of such gates of length O(log^3.97(1/ϵ)) and synthesis time in order of O(log^2.71(1/ϵ)). In the last few years extensive research resulted in great improvements both in terms of the sequence length and synthesis time. They used a diverse set of techniques (usage of ancilla or not, different libraries, approximate or exact synthesis etc). Some of the best results in terms of the generated sequence length can be found in <cit.>. These works offer a length of less than 10log(1/ϵ)  T gates (T gates are considered more costly if they are built fault-tolerantly). In the presented circuits, the worst case angle of a R_z^jk(θ)  gate is θ=2π / d^q, so the desired approximation error should be of the same order ϵ≈ 2π/d^q. Consequently, each R_z^(jk)(θ) gate can be adequately approximated by a sequence of H^(jk) and T^(jk) gates of length of the order 10log(d^q/2π)≈ 10qlogd. Thus if we have to use a constant library of components due to implementation and/or fault tolerance reasons we have a linear  in q multiplicative overhead in the quantum costs and depths calculated in Section <ref>.  


