Fast Fourier transforms and fast Wigner and Weyl functions in large quantum systems

Two methods for fast Fourier transforms are used in a quantum context. The first method is for systems with dimension of the Hilbert space $D=d^n$ with $d$ an odd integer, and is inspired by the Cooley-Tukey formalism. The `large Fourier transform' is expressed as a sequence of $n$ `small Fourier transforms' (together with some other transforms) in quantum systems with $d$-dimensional Hilbert space. Limitations of the method are discussed. In some special cases, the $n$ Fourier transforms can be performed in parallel. The second method is for systems with dimension of the Hilbert space $D=d_0...d_{n-1}$ with $d_0,...,d_{n-1}$ odd integers coprime to each other. It is inspired by the Good formalism, which in turn is based on the Chinese reminder theorem. In this case also the `large Fourier transform' is expressed as a sequence of $n$ `small Fourier transforms' (that involve some constants related to the number theory that describes the formalism). The `small Fourier transforms' can be performed in a classical computer or in a quantum computer (in which case we have the additional well known advantages of quantum Fourier transform circuits). In the case that the small Fourier transforms are performed with a classical computer, complexity arguments for both methods show the reduction in computational time from ${\cal O}(D^2)$ to ${\cal O}(D\log D)$. The second method is also used for the fast calculation of Wigner and Weyl functions, in quantum systems with large finite dimension of the Hilbert space.


I. INTRODUCTION
The fast implementation of large Fourier transforms is very important for many technological applications.Roughly speaking in this paper we express the Fourier transform in a Hilbert space of large dimension, as a combination of many Fourier transforms in Hilbert spaces of small dimension.This is a fast Fourier transform, because performing many 'small' Fourier transforms instead of one 'large' Fourier transform, is computationally beneficial.The 'small' Fourier transforms can be performed in a classical computer or as quantum Fourier transforms in a quantum computer.In the latter case, we will also have an additional well known reduction of the computational time by quantum Fourier transform circuits (e.g., [1,2]).Our methodology (and the associated reduction of computational time) is applicable to the calculation of other quantities also, like the Wigner and Weyl functions.
Two important approaches are the Cooley-Tukey formalism [3,4], and the Good formalism [5][6][7] which is based on the Chinese remainder theorem.There are also many variations of these schemes (reviewed in [8][9][10]).In this paper we study the implementation of fast Fourier transforms in quantum systems with large dimension of the Hilbert space.We also study the fast calculation of the Wigner and Weyl functions.This is an important application of the physics of quantum systems with finite -dimensional Hilbert space(e.g.[11]).
We consider a finite quantum system Σ(D) with variables in Z(D) (the ring of integers modulo d) where D is an odd integer.This system is described by the D-dimensional Hilbert space H(D).There are well known technical differences between quantum systems with odd dimension D and even dimension D (e.g., [12][13][14]).In this paper we consider systems with odd dimension D. We discuss the fast implementation of the Fourier transform F in Σ(D), using two methods described briefly below.
A. First method for the case D = d n with d an odd integer The fast implementation of the Fourier transform F in Σ(D), is using a sequence of n Fourier transforms (together with some other transforms) in a multipartite system Σ n (d) comprised of n components each of which is described with variables in Z(d).Positions and momenta in Σ n (d) take values in [Z(d)] n = Z(d) × ...Z(d) and the corresponding Hilbert space is H A = H(d) ⊗ ... ⊗ H(d).The Hilbert spaces H(D) and H A are isomorphic (they have the same dimension), and in this sense Σ(D) and Σ n (d) are two different descriptions of the same system.However, Fourier transforms and other phase space methods are different in these two cases [15].
Mathematically, this approach is inspired by the Cooley-Tukey formalism [3] for fast Fourier transforms (see also [8][9][10]), and is used here in a quantum context.But we note that the most popular Cooley-Tukey algorithm is for D = 2 n , whilst in our approach D is a power of an odd number.
A quantum circuit for the implementation of this fast Fourier transform is given in Fig. 1.In some special cases, the various operations can be performed in parallel (parallel computing).
We discuss the complexity of this method and show that the computational time is reduced from O(D 2 ) to O(D log D).We also present numerical work that supports this.
A limitation of the method is the fact that the ring Z(D) (with D = d n ) is not isomorphic to the ring [Z(d)] n .Although there is a bijective map between them, sum and products do not correspond to sums and products (section 2.A).The implications of this are discussed in section 4.C.For example, this method cannot be used in Eqs(63) below, for the fast calculation of the Wigner and Weyl functions.Mathematically, this approach is inspired by the Good formalism [5][6][7] for fast Fourier transforms (see also [8][9][10]), which in turn is based on the Chinese remainder theorem, and is used here in a quantum context.A quantum circuit for the implementation of this fast Fourier transform is given in Fig. 5.
The complexity of the method is discussed, and it is shown that the computational time is reduced from O(D 2 ) to O(D log D).This is supported with numerical work.
A strength of the method is the fact that the ring Z(D) is isomorphic to the ring Z(d 0 ) × ... × Z(d n−1 ) (section 2.B).Because of this the method is used for the fact calculation of Wigner and Weyl functions in section 6.

C. Contents
The work is complementary to the work on quantum Fourier transforms.It reduces a large Fourier transform to many small Fourier transforms, and this reduces the computational time.The small Fourier transforms can be preformed with a classical computer, or (if available) with a quantum computer so that we have the additional (and well known) advantages of quantum Fourier transforms [1,2].
In section 2 we discuss the number theory related to the two methods.In section 3 we consider a quantum system Σ(D) with variables in Z(D) where D is an odd integer, described by the D-dimensional Hilbert space H(D).In section 4 we present the first method for the case where D = d n .In section 5 we present the second method for the case where D = d 0 ...d n−1 with d 0 , ..., d n−1 odd integers coprime to each other.In section 6 we use the second method for the fast calculation of the Wigner and Weyl functions.We conclude in section 7 with a discussion of our results.

II. NUMBER THEORY FOR THE TWO FAST FOURIER TRANSFORMS
A. A bijective map between the non-isomorphic rings [Z(d)] n and Z(D) when D = d n Z(D) is the ring of integers modulo D, where D is an odd integer.We take D = d n (where d is an odd integer) and consider a bijective map between [Z(d)] n = Z(d) × ... × Z(d) and Z(D).We use upper case letters for elements in Z(D), and lower case letters for elements in Z(d).We also take j r ∈ Z(d) and J ∈ Z(D) in the 'periods' correspondingly.
We introduce the following bijective map between the sets [Z(d)] n and Z(D) Given J, we can find the j 0 , ..., j n−1 as the remainders in the following sequence of divisions: • We divide J by d and we get j 1 + j 2 d + ...j n−1 d n−2 and remainder j 0 .
We note that the [Z(d)] n as a ring (with addition and multiplication componentwise), is not isomorphic to the ring Z(D) because addition and multiplication is different [15].Indeed does not correspond to J + K.The sum in Z(D) has the 'carry' rule and the r-component might be j r + k r + 1 rather than j r + k r .In contrast, there is no 'carry' rule in [Z(d)] n .Also the multiplication in Z(D) does not correspond to the componentwise multiplication in [Z(d)] n (j 0 , ..., j n−1 ) Due to the non-isomorphism of the rings Z(D) and [Z(d)] n , there is a limitation (see subsection 4.C) of the corresponding fast Fourier transform method in section 4. We use the notation For later use, we use Eq(4) and we get Example II.1.We consider the bijective map between the sets [Z(3)] 2 and Z( 9): Then (1, 1) corresponds to 4 ∈ Z (9).A different method for Fast Fourier transforms is the Good method [5][6][7] which is based on the Chinese remainder theorem.In a quantum context it has been used in [11,16].
If d 0 , ..., d n−1 are coprime, then the ring We first define the integers b ν is the inverse of a ν within Z(d ν ), and it exists because the a ν , d ν are coprime.We also define the c ν = a ν b ν as an element of Z(D), which is an integer multiple of d ν plus one (N d ν + 1).
Proof.In the first relation, for ν ̸ = µ we get a multiple of D, which is 0(mod D).
In the second relation, we get In the third relation, we get We define a bijective map between Z(d 1 ) × ... × Z(d n ) and Z(D) as follows: (j 0 , ..., j n−1 ) ↔ J; The Chinese remainder theorem ensures that this map is bijective.Using Eq.( 12), we prove that and therefore the ring We also define a different bijective map From Eqs( 13), (15) we find the relationship between j ν and j ν : Using Eqs( 12),( 13), (15) we prove that It then follows the important relation: Example II.3.Let d 0 = 3 and d 1 = 5.Then D = 15 and a 0 = 5; b 0 = 2; c 0 = 10 Then As an example, we take J = 11 and we find the corresponding (j 0 , j 1 ) = (2, 1) and ( j 0 , j 1 ) = (4, 2).We confirm Eq(16): We consider a quantum system Σ(D) with variables in the ring Z(D), where D is an odd integer.H(D) is the D-dimensional Hilbert space describing this system.
Let |X; J⟩ where J ∈ Z(D) be an orthonormal basis in H(D).The X in the notation is not a variable, it simply indicates 'position states'.The finite Fourier transform F is given by [17] We act with F † on position states and get the dual basis The P in the notation is not a variable, it simply indicates 'momentum states'.A state |s⟩ in H(D) can be written as Below we study the fast implementation of this Fourier transform.
We consider a multipartite system Σ n (d) comprised of n components each of which is described with variables in Z(d).Positions and momenta take values in [Z(d)] n .This system is described with the d n -dimensional Hilbert space H A = H(d) ⊗ ... ⊗ H(d).We consider the basis An arbitrary state is written as |s⟩ = s(j 0 , ..., j n−1 )|X; j 0 , ..., j n−1 ⟩; |s(j 0 , ..., j n−1 )| 2 = 1. (26) Fourier transforms are defined as: We assume that D = d n and compare and contrast the systems Σ n (d) and Σ(D).Then H A is isomorphic to H(D) (because they both have the same dimension), and therefore Σ(D) and Σ n (d) are two different descriptions of the same system.However as discussed in ref [15], Fourier transforms and phase space methods (displacement operators, Wigner and Weyl functions, etc) are different in Σ(D) and Σ n (d) (F is different from F A ).This is because in these techniques we use addition and multiplication and as we explained above, the rings [Z(d)] n and Z(D) are not isomorphic to each other.Furthermore (proposition 4.4 in ref [15]), depending on the d, n, the Fourier transforms in Σ n (d) and Σ(D) are unitarily inequivalent or unitarily equivalent .
Below we explain how the equivalent of the Fourier transform F in Σ(D), is a sequence of transformations in Σ n (d) that involve Fourier transforms in the various components together with some other transformations.The latter is a fast Fourier transform in a quantum context.
Example IV.1.For n = 2, Eq.(29) becomes Acting on a vector s(K) = s(k 0 , k 1 ) we get In this case the fast Fourier transform given above becomes with C. Limitation of the method We have calculated the Fourier transform of the function s(K) = s(k 0 , ..., k n−1 ).For other functions it is not easy to apply this method.For example in the Weyl function in Eq.( 63) below, we want to calculate the Fourier transform of the function s(K)s * (B + K).Because the rings [Z(d)] n and Z(D) (with D = d n ) are not isomorphic to each other, if It is then difficult to apply directly the above formalism to Eq.( 63) for the fast calculation of the Weyl and Wigner functions.
In general, this fast Fourier transform is not directly applicable to functions which involve various sums and products of the variables.The fact that the rings [Z(d)] n and Z(D) are not isomorphic to each other, limits the practical use of the method.

D. Parallelism in the special case of factorisable states
We consider the factorisable state In this case Also Also etc.The last one is We note that for factorisable functions, we can calculate independently each of the n factors G 0 (j 0 , ..., j n−1 ), G 1 (j 0 , ..., j n−2 ), ..., g n−1 (j 0 ) and multiply them at the end.Therefore this scheme is suitable for parallel computation.The calculation in the previous subsection for general functions, needs to be done sequentially.
Remark IV.3.The 'parallel formalism' of this section is limited to special cases where we know that the factorisation in Eq.( 42) holds.Given s(K) = s(k 0 , ...k n−1 ), we give a necessary (but not sufficient) condition for the factorisation to hold.We define the Here we have a summation over all indices, except one.A necessary (but not sufficient) condition for Eq.( 42) to hold, is that E. Time complexity of the Fourier transform: counting the number of multiplications The estimate of the computational time is usually based on the number of multiplications, because they require more computational time than additions.It is easily seen that the number of multiplications for 'normal' Fourier transform is O(D 2 ) (it is a multiplication of a D × D matrix with a D-dimensional vector).For the fast Fourier transform it is known that a lower bound for the computational time is O(D log D), and we now give an approximate estimate for this.
In the fast transform in section IV B, the first step in Eq.( 31) is a Fourier transform in a d-dimensional space and it requires d 2 multiplications.This needs to be repeated for all values of the n − 1 variables k 0 , .., k n−2 which take d values each, therefore the number of multiplications is d 2 d n−1 = Dd.The second step in Eq.(32) involves another Dd multiplications (plus some extra multiplications which we ignore because we are interested in a lower limit).In this way we find that a lower bound for the number of multiplications is Many authors pointed out that this is a lower bound and that 'real' numerical fast Fourier transforms take a bit more time than that.We consider a Hilbert space H(D) with D = d 2 , where d that takes all the odd values 51, ...., 101.Using a random vector s(K) = s(k 0 , k 1 ) (produced by qiskit [19]) , we calculated s(J) = s(j 0 , j 1 ) using both Eq(38) (that involves the multiplication of a D × D matrix times a D-dimensional vector) and also the fast Fourier transform in Eqs(39), (40).We call T (D) the computational time for the calculation of all components s(j 0 , j 1 ) with the 'normal' Fourier transform in Eq.(38), and T f (D) the computational time for the calculation with the fast Fourier transform in Eqs(39), (40).In Figs.2,3 we plot • A Fourier transform of s(K) = s(k 0 , ..., k n−1 ) with ω dn−1 (j n−1 b n−1 k n−1 ) (we note here the constant b n−1 ) and summation over k n−1 : • A Fourier transform of s 1 (j n−1 |k 0 , ..., k n−2 ) with ω dn−2 (j n−2 b n−2 k n−2 ) (we note here the constant b n−2 ) and summation over k n−2 : etc.The last step is Similarly to the previous method, for factorisable functions these n steps can be done in parallel.But for general functions, they need to be done sequentially.

C. Time complexity of the Fourier transform: counting the number of multiplications
We first give an approximate estimate that a lower bound for the computational time in the present scheme, is O(D log D).
In the fast transform in section V B, the first step in Eq.( 58) is a Fourier transform in a d n−1 -dimensional space and it requires d We consider Hilbert spaces H(d 1 d 2 ) where d 1 = 53 and d 2 takes the odd values 55, 57, ...., 101.Since 53 is a prime number the d 1 , d 2 are coprime.As in section IV E we used a random vector s(K) = s(k 0 , k 1 ) (produced by qiskit [19]) , we calculated s(J) = s(j 0 , j 1 ).In Figs.6,7 we plot The result for T (D) D 2 in Fig. 6 is a horizontal line and this confirms that the computational time for the normal Fourier transform is O(D 2 ).The result for 7 is also a horizontal line and this confirms that a good lower bound for the computational time of the fast Fourier transform is approximately O(D log D).
In Fig. 8 we compare T (D) with T f (D).It is seen that T f (D) is much smaller than T (D).We checked that the Fourier transform of different random vectors give similar results.

VI. FAST WIGNER AND WEYL FUNCTIONS USING THE SECOND METHOD
Phase space methods for the system Σ(D) (Wigner and Weyl functions, etc) rely heavily on Fourier transforms.Therefore fast Fourier transforms can be used for the fast calculation of various quantities within the phase space formalism.
As an example, we consider the Weyl function W (A, B) and the Wigner function W (A, B) for the state |s⟩ = K s(K)|X; K⟩ of the system Σ(D).They are given by the following Fourier transforms(e.g., [20]): The 2 −1 = D+1 2 (mod D) for odd D. We explained in subsection IV C that the first method for fast Fourier transforms (in the case D = d n ) is not directly applicable to Eqs(63), for the fast calculation of these functions.This is related to the fact that the rings [Z(d)] n and Z(D) (with D = d n ) are not isomorphic to each other.
The second method for fast Fourier transforms (in the case D = d 0 ...d n−1 with coprime d 0 , ..., d n−1 ) is directly applicable in the fast calculation of the Weyl and Wigner functions.We present in detail the fast Weyl function.We use the bijective map in Eq.( 13 Similarly to the previous methods, for factorisable functions s(k 0 , ..., k n−1 ) these n steps can be done in parallel.But for general functions, they need to be done sequentially.

B
. Fast Fourier transform F in Σ(D) as a sequence of transformations in Σn(d) with D = d nWe use the following dual notation for functions and states in Σ(D), based on the bijective map in Eq.(2):

2 n− 1 1 =
multiplications.This needs to be repeated for all values of the n − 1 variables k 0 , .., k n−2 , therefore the number of multiplications is d 1 ...d n−2 d 2 n−Dd n−1 .The second step in Eq.(59) involves another Dd n−2 multiplications.In this way we find that a lower bound for the number of multiplications is D(d 0 + ... + d n−1 ) ≥ D(log d 0 + ... + log d n−1 ) = D log D.