Along with developing quantum algorithms for Computational Fluid Dynamics applications, an important part of the current work involves novel approaches for simulating quantum algorithms on classical hardware. Specifically, a detailed investigation into the potential of using FPGA-based hardware acceleration.

Simulation on classical hardware is essential to allow quantum algorithm researchers to test and validate new algorithms on a simulated “ideal” quantum computer (i.e., neglecting noise and decoherence errors for example) as well as in a simulated-error environment, allowing them to gain further insight into what is possible using this model. However, quantum computing simulation on classical hardware is extremely challenging in terms of computational time and memory requirement, particularly when the most general approach, i.e., storing the full quantum state vector amplitudes as used in the Schrödinger model, is employed. In general, computational time and memory requirements scale exponentially with the number of qubits, quickly reaching the limits of classical hardware.

Computing systems are becoming increasingly heterogeneous, using a variety of hardware accelerators to speed up computational tasks. One such type of accelerators, FPGAs (Field Programmable Gate Arrays), are reconfigurable devices that can be programmed after manufacturing to implement a user-specified digital circuit. They offer in particular very fine grained task parallelism and very high internal memory bandwidth, as well as excellent performance-per-Watt.

For the quantum algorithms for Computational Fluid Dynamics applications considered here, an important research question is if and how FGPAs can be effectively used in the evaluation of the quantum circuit implementations of the algorithms. To facilitate the simulation on FPGAs, this work details how quantum circuit transformations can be introduced to create multiple, smaller computational kernels by “specializing” the quantum state of one or more the most significant qubits in the quantum state vector. The underlying idea is to create a representation that will exploit the potential for fine-grained parallelism.

The main focus of the current work is the efficient simulation of large and complex quantum circuits representing lattice-based models of fluid mechanics. This way, the work also contributes to the ongoing quest to develop efficient quantum algorithms for LBM. The focus here is on a reduced lattice-based model for fluid mechanics, specifically the D1Q3 model in a modified form to facilitate quantum-computer implementation. A key feature in the present research work is the evaluation of the non-linear equilibrium distribution function using quantum computational basis encoding. For the quantum-circuit implementation of the D1Q3 model, a series of quantum-circuit transformation techniques are proposed to facilitate the efficient evaluation of the quantum algorithm while benefiting from the fine-grained parallelism offered by FPGAs. The main contributions can be summarized as follows:

• Demonstration of how quantum algorithms including non-linear terms can be derived by employing quantum computational basis encoding. Although the use of this type of encoding in general precludes exponential speed-ups relative to classical algorithms, it is expected that this approach has significant potential as part of larger quantum algorithms and when achieving polynomial speed-up is sufficient;

The structure and content of this work can be summarized as follows. Section 2 briefly reviews essentials aspects of quantum computing of particular relevance to the present work. Section 3 presents a detailed literature review focusing on quantum computer simulation techniques as well as previous work in developing quantum algorithms employing computational-basis encoding. Section 4 details the work flow used in the present work to go from quantum circuits to FPGA implementation. Section 5 briefly reviews the Lattice Boltzmann Method, introducing the concepts of “streaming” and “collision” operations. Since the most widely used LBM models, D2Q9 and D3Q27, were deemed too complex for the current proof-of-concept work, the reduced D1Q3 model is described in Section 6, followed by its modified formulation in Section 7. Section 7 also introduces the quantum floating-point approach used in the present work, introduced and detailed in previous work (Steijl, 2022). For the non-linear collision term of the D1Q3 model, Section 8 then details the design of the quantum-circuit implementation. The quantum-circuit transformation steps used to facilitate efficient simulation on FPGAs is detailed in Section 9. The computational complexity in terms of required number of qubits for ‘full’ and ‘reduced’ circuits is detailed in Section 10. Section 11 details two examples of reduced quantum circuits where the introduced transformation steps reduced the quantum circuit from 37 to 25 qubits. In Section 12, results from simulating the discussed example quantum circuits are presented and analyzed, and suggestions for improvements to the FPGA architecture are given. Finally, Section 13 presents concluding remarks and outlines steps for future research.

Typically, there are two methods of encoding the result of a quantum algorithm: encoding within the amplitudes of the quantum state and encoding within the computational basis of the quantum state. Most existing quantum algorithms employ the first method and encode the input information (as well as the output information) in the amplitudes of the quantum states. A key feature of this encoding method is that it enables the simultaneous manipulation of 2 ⁿ amplitudes by performing a single unitary operator on a n _q -qubit coherent quantum register. It is the most commonly used encoding method to take advantage of quantum parallelism in quantum computing. As an example, the widely-used Quantum Fourier Transform (QFT) uses amplitude-based encoding, where the N = 2 ^nq input defines as complex amplitude x _j with j ∈ [0, N − 1] gets transformed to an output state defined by N complex amplitudes y _j (j ∈ [0, N − 1]): ∑ j = 0 N − 1 x j j → ∑ k = 0 N − 1 y k k (1)

Encoding within the computational basis of the quantum state is not as widely used. For some tasks performed by quantum algorithms, this computational basis encoding is more efficient than the amplitude encoding. Also, for some computational tasks, such as the addition of two vectors, amplitude encoding creates problems in case the two vectors cancel each other and create a 0 vector that cannot be defined in terms of quantum amplitude encoding. As, a result, quantum computing in the computational basis (QCCB) is widely-used in quantum algorithms performing arithmetic operations. Also, in some quantum algorithms involving the Fourier transformation, the Fourier coefficients may be needed in the computational basis (Zhou et al., 2017). The quantum algorithm for computing the Fourier transform in the computational basis (termed QFTC) by Zhou et al. (2017) enables the Fourier-transformed coefficient to be encoded in the computational basis as follows, k 0 → k y k (2) where y k corresponds to the fixed-point binary representation of y _k ∈ (−1, 1) using two’s complement format. In the algorithm proposed by Zhou et al. (2017), the input vector x ⃗ is provided by an oracle O _x such that, O x 0 = ∑ j = 0 N − 1 x j j (3) which can be efficiently implemented if x ⃗ is efficiently computable or by using updatable quantum memory (qRAM) that takes complexity log(N) under certain conditions (Zhou et al., 2017). Comparing Eqs 1, 2, it is clear that encoding in the computational basis requires a number of additional qubits depending on the required fixed-point representation.

In terms of matrix computations, Ma et al. (2020) introduced a quantum algorithm for the computation of the QR matrix decomposition using computational basis encoding. The simulation time of the algorithm shows a polynomial speed-up relative to classical algorithms. In the QCCB-based QR algorithm, real numbers are encoded using a fixed-point representation. This means that any real number with its amplitude bounded by positive integer M can be approximate to precision O(M/2 ^m ) using m bits. A key factor in achieving the polynomial time complexity gain relative to classical algorithms, is that this representation accuracy is independent of the matrix size considered. In addition to the quantum QR decomposition algorithm, Ma et al. (2020) also proposed a general approach to simulate any quantum algorithm based on amplitude encoding by using the QCCB. The authors show that for this QCCB simulation, the simulation time does not exceed O(N ² polylog(N)) times the cost of the original quantum algorithm based on amplitude encoding. It should be noted that to achieve this performance, the QCCB-based algorithms introduced by Ma et al. require an updatable quantum memory (qRAM) where the cost of updating n _up entries is O(n _up log(n _up )). Such an updatable quantum memory model was previously used and investigated by Kerenidis and Prakash (2020).

Based on the published works using the quantum computational basis encoding, it is clear that this approach has not been as widely explored as amplitude encoding. The main reason is that the quantum parallelism and exponential saving in storage offered by quantum amplitude-encoding is lost. However, it is clear that important and meaningful quantum algorithms using computational-basis encoded can still be obtained when used as part of a larger quantum algorithm or in cases where a suitable-designed updatable quantum memory is available.

The most general approach to simulating the quantum state vector in a quantum computer employs the Schrödinger wave function approach where the coherent quantum state of an n _q -qubit register is defined by N = 2 ^nq complex amplitudes. In the following, this approach is termed “full state-vector approach.” Alternatively, the exponential scaling can be addressed by avoiding the Schrödinger wavefunction representation as used in the full-state vector approach. Aaronson and Gottesman (2004), Garcia and Markov (2015) and more recently, Lee et al. (2018), used the Heisenberg representation of the quantum circuit via Stabilized Frames. This approach is limited to specific quantum circuits, i.e., termed stabilizer circuits, and for these cases has a polynomial scaling of memory and computational cost with increasing qubits. However, this approach cannot be used for general quantum circuits. Therefore, in the present work, the focus is on simulation techniques based on Schrödinger wavefunction representation. This approach can be used for quantum algorithms using the quantum-amplitude encoding as well as for algorithms using computational-basis encoding.

Early works related to the development of massively parallel quantum computer simulators include De Raedt et al. (2007) and Tabakin and Julia-Diaz (2009). Both works involve highly optimized simulators based on the full state vector approach, implemented in Fortran using the MPI library to exploit parallelism on distributed-memory architectures. For a range of textbook examples, very good parallel scaling was observed for tests with upto 4,096 processors documented in De Raedt et al. (2007). The work on QCMPI reported by Tabakin and Julia-Diaz (2009), was mainly motivated to facilitate numerical examination of not only how QC algorithms work, but also to include noise, decoherence, and attenuation effects and to evaluate the efficacy of error correction schemes. Both works focus on evaluating quantum circuits, while earlier work by De Raedt et al. (2000) involved a Quantum Computer Emulator designed to simulate the physical realization of the quantum computer and a graphical user interface to program and control the simulator. The potential to speed-up large-scale parallel simulations such as those described by De Raedt et al. (2007)and Tabakin and Julia-Diaz (2009) using Graphics Processing Units (GPUs) has been thoroughly investigated in the past decade. Early work includes Gutiérrez et al. (2010) and Lu et al. (2013). In the work of Gutierrez and co-workers, the simulation of an ideal quantum computer using NVIDIA’s CUDA GPU programming model is discussed in the context of how such a problem can benefit from the high computational capacities of GPUs. The simulator discussed in their work takes into account memory reference locality issues. The work showed that the challenge of achieving a high performance depends strongly on the explicit exploitation of memory hierarchy.

Highly-optimized approaches of the full state vector approach aimed at avoiding the exponential scaling of memory and computational cost with increasing number of qubits have been investigated for more than 2 decades. One possibility (e.g., Viamontes et al., 2003; Rosenbaum, 2010) involves employing the Schrödinger wavefunction representation (as used in the full-state vector approach) along with compact representation of amplitudes using tree-based or decision-diagram based data structures. For a range of practically relevant quantum algorithms, significant memory and time savings were documented relative to the full-state vector approach. However, worst-case situations often occur where memory and time complexity are still exponential with number of qubits.

The earliest work in simulating quantum circuits on FPGAs dates back to Khalid et al. (2004). A compiler which produces VHDL, a hardware description language, code that emulates the quantum circuit on the FPGA was developed. It emulates quantum parallelism by constructing parallel data paths for the state vector amplitudes representing the qubits, i.e. implementing the whole quantum circuit in the FPGA fabric. State vector amplitudes are implemented by fixed point numbers to keep the size of the circuits manageable. Fixed point was also chosen since the probability amplitudes can only have a decimal part of 0 or 1. This approach emulates a full quantum circuit on the FPGA, requiring a full synthesis when changing the circuit.

The approach used in Aminian et al. (2008) divides quantum circuit simulation into two circuit types based on gates used in the circuit. For circuits involving only X, Y, Z, and CNOT, they reduce the Logic Cell (LC) usage required for each type of gate to a handful (X: 2, Y: 6, Z: 2, CNOT: 4). They do this by adding extra information bits (basis, complexity, sign) and simply operating on them when applying any of these gates (however there is more basis bits in the case of CNOT). The second group is H, PS (phase shift), and CR (controlled rotation), which are implemented directly as adders and multipliers, requiring resources which increase with the number of mantissa bits. For circuits involving both groups of gates, they apply a different simulation policy than when just XYZC gates are used.

Conceicao and Reis (2015) address the issue of re-synthesis present in prior works. They present a reusable architecture for which synthesis is only rerun when the number of qubits or mantissa bits of the fixed point representation is required to be changed. In their design, a control unit holds an address of a an instruction in some instruction memory (list of gates) and a quantum Arithmetic Logic Unit (ALU) is fed a gate operational code (opcode), target qubit, and two control qubits at each gate and then communicates with a quantum register to perform the gate. They report their LC usage for a number of algorithms and benchmarks. In terms of LC usage, they are outperformed by Aminian et al. (2008), which they point out, but also observe that the ratio between their average usage of logic cells decreases in comparison when increasing the number of qubits from 3 to 8, leading them to believe their system would be better when scaled up.

Lee et al. (2016) developed a serial-parallel architecture-based FPGA emulation framework for quantum computing and, for small numbers of qubits, demonstrated significant speed-ups relative to CPU-based emulations.

Pilch and Dlugopolski (2019) proposed, designed and implemented an easily scalable universal quantum computer emulator, focused on reflecting natural quantum processes in hardware, while maintaining the time complexity of quantum algorithms. The underlying idea is to move the weight of complexity from time to hardware resources by using the inherent parallelism of FPGAs. As proof-of-concept, the authors created a hardware-software system capable of running and correctly interpreting results of the Deutsch quantum algorithm. So far, only small circuits where considered, e.g., the Deutsch algorithm was emulated for a 2-qubit quantum computer.

The works discussed so far (Khalid et al., 2004; Aminian et al., 2008; Lee et al., 2016; Pilch and Dlugopolski, 2019) demonstrate the promise for emulating quantum circuits on FPGAs, albeit for low number of emulated qubits. Mahmud and El-Araby (2019) focus on scalability, presenting two architectures for emulation. The first is a CMAC (complex multiply-and-accumulate) unit-based system, which for a given quantum circuit, relies on having the full algorithm matrix precomputed. An optimization to this is to have a kernel which dynamically generates the values of the algorithm matrix, massively reducing the memory requirement. Using this architecture, Mahmud and El-Araby (2019) emulated a 20-qubit QFT, an increase in qubits compared to previous works in FPGA emulation of quantum circuits. This required the creation of a custom hardware architecture for generating the values of the QFT matrix. The second recognizes that there may be algorithms which have sparse algorithm matrices which may not be suited for the CMAC-based architecture. Instead, this architecture requires a custom acceleration kernel to be developed from the quantum algorithm, which is then applied to the input state vector. Using this architecture, they emulated a 30-qubit Quantum Haar Transform. This required the extraction of a simplified kernel from the mathematical description of the QHT rather than from its quantum circuit description. This is a considerably higher number of qubits than those achieved in previous works. However, no method of automating the generation of the kernels from a quantum circuit model description is discussed. The authors extend this method to Grover’s database search algorithm in Mahmud et al. (2020).

Khalid et al. (2021) describe a proposal for a resource-efficient FPGA-based abstraction of quantum circuits. A non-programmable embedded system capable of storing, measuring, and introducing a phase shift in qubits is implemented. The proposed single-input single-output architecture implements single-input gates with corresponding memory and measurement blocks. A fixed-point quantum gate representation is used, using 8 bits (2-bit integer, and 6-bit fraction). By increasing the number of bits used for qubit representation, the quantized superposition states of the of qubit increase, leading to enhanced accuracy of the emulation results. The main objective of the proposed abstraction was to provide an FPGA-based platform as the fundamental sub-block for the design of quantum circuits. The quantum key distribution algorithm BB84 was implemented using the proposed platform as a proof-of-concept. The proposed design exhibits two principal properties of quantum computing, i.e., parallelism and probabilistic measurement.

The concept of using exponentially-increasing resources (with problem size) on an FPGA to maintain the exponential time-complexity gain of quantum algorithms relative to their classical counterparts was also investigated by Bonny and Haq (2020) who implemented the quantum k-means clustering algorithm on an FPGA emulator. Clustering is a technique involving the classification of unlabeled data into a number of categories, and is widely used in machine learning and data mining. The main computational work in the k-means clustering algorithm is the computation of the distance between points. Bonny and Haq model the points as n-dimensional vectors on the Bloch sphere and then use the inner product as an estimation of the distance between two vectors. In the example implementation 2D data was used. The work presented forms an example of a quantum-inspired algorithm allowing (exponential) speed-up relative to classical algorithms even when running on classical hardware, by trading time-complexity and resource use on the FPGA. Clearly, the rapid (exponential) increase of logic-gate resources with problems size limits this approach to relatively small problems. For quantum circuit emulation, this means that only a few qubits can be used when trying to main exponential time-complexity gain. This work follows on from previous works into such quantum-inspired algorithms including work on modified Grover’s search algorithm (Aïmeur et al., 2007) as well as quantum-inspired evolutionary algorithms for optimization problems (Han and Kim, 2002). More recently, Fujitsu’s quantum-annealing inspired optimizer as described by Aramon et al. (2019) uses extensive hardware acceleration techniques to achieve time-complexity improvements.

This summary of works involving hardware acceleration of quantum computing simulation shows that there is a growing interest in simulating quantum circuits on FPGAs and their results show that there can be a considerable computational advantage to using an FPGA to simulate a quantum computer. However, most of the research so far only considered circuits with few ( < 10 ) qubits and also did not consider circuit transformation techniques for reducing the number of qubits. The work which demonstrated the most promise for scaling to a high number of qubits recently is Mahmud and El-Araby (2019) and Mahmud et al. (2020) specialized kernel-based approach, using which they ran a 30-qubit QHT, and a 32-qubit Grover’s search circuits on an FPGA. While this approach reaches the highest number of qubits simulated on an FPGA in the literature, of which we are aware, it is not very easily reusable. Using their more circuit-independent CMAC-based approach, they were able to simulate a 20-qubit QFT circuit. Our work is in line with their more generic approach because, as we discuss in the next section, reusability is one of our primary goals.

We developed a toolchain in the functional programming language Haskell to allow for efficient specification and testing of quantum circuits. The toolchain includes an embedded Domain-Specific Language (eDSL) that allows for specification and static analysis of quantum circuits. This embedded language allows users to label qubits and use high-level circuit constructs, like looping, tiling, etc. The functionality offered by this toolchain is similar to other functional language-based quantum computing tools like Quipper (Green et al., 2013). It was implemented from scratch, however, to have maximum control over the circuit reductions and any optimizations necessary for an FPGA-based architecture.

Since anything the Haskell tool would do on any form of quantum circuit input would be static compile-time checks/processing (nothing close to exponential or related to the number of qubits involved), then the host code controlling the FPGA environment should be capable of doing the same in negligible time compared to the execution time of the quantum circuit on the simulator, which will necessarily be exponential to maintain generality. However, writing the compiler step in a separate toolchain allows for more modularity of our designs and separation of concerns. While currently the toolchain has no functionality for adding architecture-specific optimizations (such as gate fusion or other system-level changes requiring a higher level instruction set), it acts as a starting point for future research. In addition, usually (and in this case) the host code is written in C++; writing and maintaining an eDSL in Haskell will be considerably easier and less time consuming.

The toolchain offers several constructs for looping gates over sets of qubits, managing complex controls, and defining intermediate circuit components. Some of these constructs are demonstrated in Listing 1.

fullAdd :: QReg -> QReg -> Qu -> Qu -> Circ

Listing 1

Generic input size full Cuccaro adder and square circuit example implementation in the presented Haskell eDSL.

After an FPGA board configuration is compiled, the quantum circuit simulation flow starts with a user-specified circuit description. We can describe the circuit using an eDSL in Haskell, allowing for complex higher-level circuit operators to be used; alternatively, the toolchain can read a QASM-like file specifying the circuit. This process is demonstrated in Figure 1. The toolchain first verifies the specified circuit, ensuring all qubits used are valid (have an index in the register) and no gates are specified with invalid target/controls. Then the named qubit identifiers are parsed away and the qubits are mapped to an index in the quantum register. At this point in the compilation process, some constructs are still available to the toolchain which would not necessarily be available to the FPGA (like direct calls to a SWAP gate, or a high number of controls), which need to be reduced away. SWAP gates are replaced with their equivalent CNOT specifications, negative controls are reduced by negating the control qubits before and after the gate, and gates with a higher number of controls than is supported are expanded to several gates with fewer controls. This results in a circuit which is ready to be converted to a “QP” file which is simply a list of integers specifying the circuit. The first element in this list is always the number of qubits required for the circuit. Taking into account the maximum number of controls allowed by an architecture (C), each emitted gate consists of its opcode, target qubit, followed by a constant number of controls. The resulting list is then written to disk, ready to be read by the simulator host.

HLS design tools for FPGAs (like Intel’s Quartus/Intel’s FPGA OpenCL Compiler and Xilinx’s SDAccel toolchain) have been growing in popularity and effectiveness over the past decade; making it possible to develop customized FPGA architectures for domain-specific problems in a relatively high-level language like OpenCL, without requiring knowledge of an HDL like Verilog.

The presented simulator is written in OpenCL and built using Intel/Altera SDK for OpenCL. We used Altera Offline Compiler (aoc) version 17.1. The architecture we have implemented is a baseline implementation of a full state vector-based quantum circuit simulator. The simulator is universal, with direct on board support for the H, X, Y, Z, and R _m (integer-parametrised phase gate) gates and any controlled variations of them, up to a static maximum number of controls per gate. We currently make use of the dynamic scheduling of the NDRange OpenCL construct to submit gate applications to the kernel architecture on the board. Figure 2 show the steps involved in a circuit simulation. The QP instructions generated from the compilation toolchain are read from disk on the host CPU. These instructions contain the qubit count of the circuit, and based on this sufficient memory is then allocated on the host and device. The current representation used for the state vector storage in memory is using 32-bit floating-point complex numbers, meaning we need to allocate 2 ⁿ × (2 × 4) bytes to simulate the given circuit. Based on user input at runtime, the initial quantum state of the system is set and then copied to the device DRAMs.

Currently all FPGA architectures we have explored are parametrised by some values at compile time; primarily these are the number of compute units (NCU), which are the duplicated kernels that perform the CMACs required for the gate iteration computations, and the maximum number of controls per gate (NCONTROLS). Other more customized architectures may have more than these parameters at compile time, particularly for caching and manual memory buffering. Note that the NCONTROLS parameter determines the format of the input QP instructions; which the Haskell tool facilitates by generating QP instructions with a fixed total number of controls per gate. For non-pipelined architectures, these duplicated kernels may each include their own memory interfaces. However, for pipelined architectures (generally seen as more efficient), the compute units will generally receive their inputs over channels or pipes.

As discussed above, our focus is optimising the QWM simulation method with the full state vector being stored in memory. The main difficulty in performing a gate application in parallel arises from having an exponential memory space which has to be fully iterated through for a single gate computation (each control reduces the required accesses by half, but in general we still at least access an exponential subset of the memory space).

In the current implementation of the D1Q3 model, the idea is to represent the (scaled) distribution function values g _i , i ∈ [0, 3] and u/c) with a specialized (biased) quantum floating-point format. The motivation for choosing this format is the wider range of numbers that can be represented than in a fixed-point representation for the same number of qubits used. To minimize quantum-circuit width, a reduced number of mantissa and exponent bits is used as compared to IEEE-754 single-precision format. However, for the considered problems a scaling was used so that the ranges of numbers to be represented is both limited and predictable. For the D1Q3 model, where the numbers will be ≪ 1, this choice of parameters is further detailed later this section.

The quantum floating-point representation used here builds on earlier work by Steijl (2022) and involves reduced-bit representations of exponent and mantissa relative to single-precision in the IEEE-754 standard to facilitate quantum circuit implementations on current and near-future quantum hardware with relatively small number of qubits ( < 100 ) . A key feature of the used quantum floating-point representation is that following the IEEE-754 standard, it employs sub-normal numbers and consistent rounding (here, rounding-down to nearest). The number of mantissa qubits is defined by N _M , where only N _M − 1 mantissa qubits are stored following the “hidden-qubit” approach from IEEE-754. Then, the number of qubits storing the exponent is defined by N _E . In the present work, N _E = 3, and exponent 0 (exponent qubits in state 000 ) represent sub-normal numbers and zero as in the IEEE-754 standard. The maximum value for exponent is 7 (exponent qubits in state 111 ) refers to “overflow” conditions, as used in the IEEE-754 standard. For N _E = 3, an “unbiased” exponent formulation would be equivalent to an exponent bias of 3. To optimize for the small numbers occurring in considered problems, a bias toward smaller numbers is used here. Clearly, the choice of N _M is crucial in achieving the required accuracy. For equilibrium distribution functions, terms linear and quadratic in u/c are combined, and N _M = 3 was found to lead to excessive rounding or truncation for most values of u/c. Clearly, N _M ≥ 4 should be used. However, since the lattice-based models considered here represent conservation of mass, momentum and energy in the fluid, rounding of numbers will have a significant effect on accuracy, particularly in multiple time-step simulations. Therefore, realistically it can be expected that N _M ∈ (Fillion-Gourdeau and Lorin, 2019; Bonny and Haq, 2020) is needed for realistic engineering applications. In this work, circuits for N _M = 4 are shown as illustration, and the complexity analysis shown in Section 10 addresses in detail how the circuit width increases with N _M .

Following IEEE-754, signed floating-point numbers are stored as “ sign | exponent | mantissa ,” while the sign qubit is omitted where unsigned numbers are used in the quantum-circuit implementation. For “signed” additions or subtractions of two numbers, two (N _M + 1)-qubit registers as input to a modulo adder are created as follows. First, the hidden qubits are added to the mantissa, followed by re-normalization (to account for possible difference in exponent) to create two (N _M + 1)-qubit inputs with 0 as most-significant qubit. Where required a conversion to 2’s complement is performed, so that negative numbers have the most-significant qubit in state 1 . After addition/subtraction, the outcome is converted back to the quantum floating-point format, where a sign-qubit defines the sign and with mantissa no longer in 2’s complement, i.e., back into (N _M − 1)-qubit “hidden-qubit” representation.

For N _M = 4 and N _E = 3, the sub-normal numbers with the corresponding “negative” (N _M + 1)-qubit mantissa representation using 2’s complement method are shown in Table 1 for “bias = 3” and “bias = 8.” Here, “bias = 3” corresponds to the “standard” floating-point format with a symmetric bias. As can be seen from Table 1, with symmetric bias (“bias = 3”), terms involving u ² (O(10^–2) will always be sub-normal and often truncated to 0. For “bias = 8,” it can expected that u ² can be represented either as non-zero sub-normal or normalized numbers.

For N _M = 4 and N _E = 3, the normalized numbers for e 2 | e 1 | e 0 = 110 with the corresponding “negative” using 2’s complement method are shown in Table 2 for “bias = 3” and “bias = 8.”

For the maximum exponent ( e 2 | e 1 | e 0 = 110 ) , the values for “bias = 8” shown in Table 2 indicate that despite the greatly reduced value of the maximum number that can be represented relative to the symmetric bias case (“bias = 3”), the components of distribution function g ⃗ and velocity u/c will be so small that these can still be represented without risking an overflow, as discussed next.

The choice of suitable values for N _E and exponent bias for the scaled and normalized D1Q3 model can be made based on the flow physics that is modelled. For the D1Q3 model, the lattice speed of sound is defined as c s = c / 3 . The iso-thermal model used here was derived for flows with weak or negligible compressibility effects. For a compressible fluid, a local Mach number can be defined as: M = | u | / c s = 3 | u | / c . For the weakly compressible-flow conditions it is required that M < 0.3 or smaller. Therefore, u/c should generally be limited to maximum values of 0.1–0.15. Clearly, in the modified and re-scaled D1Q3 model employing a symmetric bias will introduce a range of numbers far from optimal for the D1Q3 model. As shown in Table 2, for N _E = 3, a bias of 8 gives ample margin at the upper end of the floating-point range when representing u/c. However, since the (scaled) and re-normalized equilibrium distribution functions combine terms O(u) and O(u ²) it was decided that a bias of 9 could potentially leave too little margin in the floating-point representation of g _i components.

In Lattice Boltzmann models, the non-linear convection term of the Navier-Stokes equation is effectively (partly) replaced by a linear streaming of the distribution functions from one lattice point to a neighboring lattice point during a time step. For the modified D1Q3 model, only two directions specified by d v 1 | d v 0 = 00 (“−1”) and d v 1 | d v 0 = 11 (“+1”) need to be considered since the remaining two are “at rest.” Assuming a uniformly-spaced one-dimensional domain with 64 lattice points and periodic boundary conditions, a possible quantum-circuit implementation for the streaming operations is illustrated in Figure 3. This type of linear operators was previously detailed by Todorova and Steijl (2020) and will therefore not be analyzed further here.

The quantum-circuit was designed with input data encoded in qubits at the top of the circuit (most significant qubits), followed by qubits representing the output of the computation performed. The remaining qubits further “down” (i.e., less significant in the employed memory indexing) generally act as ancillae qubits or as workspace. These ancillae and workspace qubits are all initialized in state 0 and will be returned to 0 at the time of completion of the quantum circuit. In later sections, a modified and transformed design will be considered were the ordering is altered to facilitate efficient simulation. For the original circuit, the qubit register for N _M = 4 and N _E = 3 can be summarized as follows: d v 1 | d v 0 indices for  4  discrete velocities e u 2 | e u 1 | e u 0 3 - bit representation of exponent of  u s u sign bit of input velocity  u m u 2 | m u 1 | m u 0 3 - bit representation of mantissa of  u e g 2 | e g 1 | e g 0 3 - bit representation of exponent of  g e q s g sign bit of output  g e q m g 2 | m g 1 | m g 0 3 - bit representation of mantissa of  g e q e s q 2 | e s q 1 | e s q 0 3 - bit representation of exponent of  u 2 m s q 2 | m s q 1 | m s q 0 3 - bit representation of mantissa of  u 2 c u t state  1  defines that  u 2 is truncated to 0 r 4 | q u 3 | r 3 | q u 2 | r 2 | q u 1 | r 1 | q u 0 workspace qubits ordered for  4 - qubit addition c workspace qubit named to represent ’carry’ bit in  4 - qubit addition r 7 | r 6 | r 5 workspace qubits ordered for  4 - qubit squaring a n c workspace qubit  -  mostly used as control qubit

This shows that 37 qubits are required in this original design. Data related to g ₀ and g ₃ are stored in the vector for d v 1 | d v 0 = 00 and d v 1 | d v 0 = 11 . Similarly, states with d v 1 | d v 0 = 01 and d v 1 | d v 0 = 10 represent data for the (identical) distribution functions g ₁ and g ₂. For the floating-point representations of input u and each component of output g ⃗ e q , 7 qubits are needed using the hidden-qubit approach. For the temporary storage of u ², the sign bit can be omitted. In the following, for N _M = 4 and N _E = 3 an exponent “bias = 8” is used in the floating-point representations.

Figure 4 shows the first part of the quantum circuit designed to compute the equilibrium distribution function g ⃗ e q for N _M = 4 and N _E = 3. The first step involves setting i c u t = 1 for the cases with a guaranteed truncation of u ²—for N _M = 4 and N _E = 3 this truncation always occurs for e u 2 | e u 1 | e u 0 = 000 , 001 and 010 . In the next step, the mantissa of u is set into the required positions in the workspace, defined by q u 3 | q u 2 | q u 1 | q u 0 , for all cases without truncation of u ² to 0. The operation SQ4 then computes the square of the mantissa and stores the results in r 7 | r 6 | … | r 0 . Operation CC _u2 then creates a (temporary) “copy” of the u ² defined in the quantum-floating format in the qubits e s q 2 | e s q 1 | e s q 0 (defining exponent) and m s q 2 | m s q 1 | m s q 0 (defining mantissa using hidden-qubit approach). For u ² no sign qubit is needed. To clear the workspace for further use, operation ISQ4 un-computes the mantissa squaring, followed by the removal of the copy of the mantissa of u from qubits q u 3 | q u 2 | q u 1 | q u 0 .

The quantum-circuit implementation for SQ4 and ISQ4 are shown in Figures 5, 6, respectively. Here, FAdd represents a 4-qubit Cuccaro full adder, and Rmv the un-computation of this adder. The required shift in the used shift-and-add approach are performed by Sh (in SQ4) and Sh’ (shift in reversed direction in ISQ4). Following the definition of u ² in quantum floating-point format, the equilibrium distribution for directions e ₁ and e ₂ (defined by d v 1 | d v 0 = 01 and d v 1 | d v 0 = 10 ) is defined using the operator 0110 in Figure 4.

To define the equilibrium distribution functions for directions e ₀ and e ₃ (defined by d v 1 | d v 0 = 00 and d v 1 | d v 0 = 11 ), the addition of ± u and u ² is required. This addition operator for the signed values defined in the N _M = 4 and N _E = 3 format is performed using a modulo-5 Cuccaro adder, denoted by MA5. Operation 0011a initializes this addition step, followed by the operation 0011b that uses the created result to define the equilibrium distribution functions for directions e ₀ and e ₃ in qubits e g 2 | e g 1 | e g 0 (exponent), s g (sign qubit) and m g 2 | m g 1 | m g 0 (mantissa qubits using hidden-qubit approach).

Figure 7 shows the second part of the quantum circuit designed to compute the equilibrium distribution function g ⃗ e q for N _M = 4 and N _E = 3. The first step involves the un-computation of the modulo-5 addition (denoted by UMA5), followed by operation 0011c used to clear the inputs to this addition. Upon completion of operation 0011c, the 14 workspace qubits are all in state 0 . However, at this stage, the temporary copy of u ² still resides in qubits e s q 2 | e s q 1 | e s q 0 and m s q 2 | m s q 1 | m s q 0 . To clear these qubits to state 0 , the square of the mantissa of u ² needs to be re-computed using SQ4. Then, CC _u2 is used to set the 6 qubits defining u ² in quantum floating point format to 0 . Then, ISQ4 un-computes the mantissa squaring step. Finally the remaining workspace qubits can be cleared along with the qubit c u t , which for cases with truncation of u ² to 0 is re-set to 0 .

At this stage, the output of the quantum circuit has the required format: the qubits e g 2 | e g 1 | e g 0 (exponent), s g (sign qubit) and m g 2 | m g 1 | m g 0 define the equilibrium distribution function for all four directions in the modified D1Q3 model, while the rest of the qubits is left unchanged.

In operator 0110, the previously computed term u ² is used to set −u ²/2 for the two “rest” directions ( d v 1 | d v 0 = 01 and d v 1 | d v 0 = 10 ). For these directions the “−” sign can be trivially introduced by setting the sign qubit of g e q , s g = 1 . For the direction 0 defined by d v 1 | d v 0 = 00 , the equilibrium distribution function is of the form −u/2 + u ²/2. For this direction (and for direction 3 with distribution function of the form u/2 + u ²/2), first the terms −u + u ² and +u + u ² are computed using a 5-qubit modulo adder MA5 (and its reverse UMA5), while the division by 2 is introduced when setting the result in quantum-floating point format. The sign change to ′ − u′ for direction 0 is created by switching to 2’s complement notation or the 4 mantissa qubits. Based on the assumption of positive input velocity u, the 5 input qubits (as “a” input to Cuccaro modulo adder) representing mantissa of u are initially set as 0 | m u 3 | m u 2 | m u 1 | m u 0 (with m u 3 = 0 for sub-normal numbers). The quantum-circuit implementation of operator SgnA5 performing this sign change is shown in Figure 8. With u ≪ 1 in the D1Q3 model, the term −u + u ² will always be a negative number. Since the 4 mantissa qubits of g ⃗ e q are not stored in 2’s complement formulation for negative values (i.e., the sign is defined by s g ), a similar sign change to SgnA5 needs to be applied on the output of the 5-qubit modulo adder. With the guaranteed negative result for direction 0, it is sufficient to perform the 2’s complement conversion on only to the 4 qubits (for N _M = 4 considered here) used to define mantissa qubits of g ⃗ e q . The operator SgnB4 performs this change. Its quantum-circuit implementation is also shown in Figure 8.

The transformation steps detailed in the previous section enabled a reduction to computational kernels with 27 qubits as compared to the 37-qubit full circuit. For a further reduction, the 14-qubit workspace needs to be transformed. Specifically, the arithmetic operations in SQ4, ISQ4, MA5 and UMA5 need to be transformed and specialized for specific inputs. Since the Quantum Computer simulator used here employs a memory allocation with the top qubits in the circuits shown acting as the most significant qubits, the overall circuit design shown in Figures 4, 7 needs to be changed: the qubits defining workspace need to be moved towards the top of the quantum circuit, and therefore, the qubits storing u ² and g ⃗ e q in quantum floating-point format need to be moved down toward less significant bit locations.

Following this re-ordering, a further reduction requires specializing (and factoring out) one or more mantissa qubits, starting from the most significant mantissa qubit of u, i.e., q u 3 , followed by m u 2 , etc. Using this approach, requires transformation to the u ² computations as well as modulo-additions, as discussed in following sections.

To illustrate the further reduction of the number of qubits, the computation of u ² is analyzed first. In the interest of clarity, this section will consider the reduced circuits for the evaluation of g 1 e q ( d v 1 | d v 0 = 01 ) and g 2 e q ( d v 1 | d v 0 = 10 ) , since these identical terms only involve the term u ²/2. Figure 9 shows the reduced quantum circuit with 26 qubit resulting from eliminating the two most-significant mantissa qubits in squaring operations for N _M = 4 and N _E = 3. As can be seen, the circuit only involves the two mantissa qubits q u 1 | q u 0 , and SQ4 and ISQ4 represent reduced squaring (and un-computation) for specific choices of q u 3 | q u 2 . For the further reduction by 2 qubits to 26-qubits, the operations SQ4 and ISQ4 involve 12 qubits in the workspace. The operation CC _u2 sets the squared-velocity values in terms of quantum-floating point format in e s q 2 | e s q 1 | e s q 0 (exponent) and m s q 2 | m s q 1 | m s q 0 (mantissa), and for the considered reduction to 26 qubits, the quantum circuit implementations for CC _u2 are detailed for e u 2 | e u 1 | e u 0 = 011 , 100 or 101 in Figure 10. For e u 2 | e u 1 | e u 0 = 110 , the operator is identical to that for 101 , apart from the NOT operation on e s q 1 that for e u 2 | e u 1 | e u 0 = 110 is performed on e s q 2 instead. Figure 11 shows the quantum-circuit definitions of the reduced mantissa-squaring and adders (as well as the reverse circuits), for the 26-qubit reduced quantum-circuit implementation.

Figure 12 shows the reduced quantum circuit with 25 qubit resulting from eliminating the three most-significant mantissa qubits in squaring operations for N _M = 4 and N _E = 3. As can be seen, the circuit only involves the mantissa qubit q u 0 , and SQ4 and ISQ4 represent reduced squaring (and un-computation) for specific choices of q u 3 | q u 2 | q u 1 . For the further reduction by 3 qubits to 25-qubits, the operations SQ4 and ISQ4 involve 11 qubits in the workspace as can be seen in the quantum-circuit implementations shown in Figure 13. For the reduced quantum circuit with 25 qubits, Figure 14 shows the quantum-circuit implementation of the “shift” operations used in the “shift-and-add” approach. For N _M = 4, a further reduction to 24 qubis can be made, by reducing out qubit q u 0 from the SQ4 and ISQ4 operations. For this further reduction to 24 qubits, quantum-circuit implementation of operations SQ4 and ISQ4 are shown in Figure 15. The addition and remove operators are specialized for the 4 mantissa qubits q u 3 | q u 2 | q u 1 | q u 0 = 1001 .

The operation CC _u2 defines squared-velocity values in terms of quantum-floating point format, and for the considered reduction to 25 qubits, the quantum circuit implementations follows directly from those discussed for the reduction to 26 qubits, since q u 1 | q u 0 are not involved in this step. The operation 0110 finally sets the equilibrium distributions functions g 1 e q (for d v 1 | d v 0 = 01 ) and g 2 e q (for d v 1 | d v 0 = 10 ) in quantum floating point format, based on the definition of u ² in floating-point format defined previously. The quantum circuit implementation of 0110 for the reduced circuits with 25 and 26 qubits (identical in both cases) is shown in Figure 16.

Following the analysis of the evaluation of the u ² terms in quantum circuits with partial reduction of the workspace qubits, the focus now moves to the more complex case of also reducing the modulo adder (and its reverse) used in computing the summations u + u ² and −u + u ². The use of 2’s complement method in the current implementation of g ⃗ e q evaluation poses particular challenges in the further reduction process, as detailed in this section. In the interest of clarity, this section will consider the reduced circuits for the evaluation of g 0 e q ( d v 1 | d v 0 = 00 ) and g 3 e q ( d v 1 | d v 0 = 00 ) , since these represent the only directions for which the modulo-addition steps is required. Figure 17 shows the reduced quantum circuit with 25 qubit resulting from eliminating the three most-significant mantissa qubits in the arithmetic operations for N _M = 4 and N _E = 3. As can be seen, the circuit only involves the mantissa qubit q u 0 . The operation 0011a prepares the modulo-adder step and depends on the mantissa qubits of u as well as on the exponent of u (represented by e u 2 | e u 1 | e u 0 ). Similarly 0011c un-computes these steps to reset workspace to 0 after completion of the addition step and setting g 0 e q and g 3 e q in quantum floating point format. The (reduced) modulo-5 adder MA5 (and its reverse UMA5) only depend on the specific choice of mantissa qubits that were reduced, i.e., in this step there is no dependency on exponent e u 2 | e u 1 | e u 0 . Based on the outcome of modulo-5 adder, the operation 0011b sets g 0 e q and g 3 e q in m g 2 | m g 1 | m g 0 (mantissa), s g (sign) and e g 2 | e g 1 | e g 0 (exponents) for d v 1 | d v 0 = 00 and d v 1 | d v 0 = 11 , respectively.

Following the process detailed in the previous section, for a further reduction by 2 qubits to 26 qubits would result in a similar quantum circuit, then involving q u 1 and q u 0 . For brevity, this quantum circuit is not shown here.

The operation 0011a comprises two parts. The first part uses quantum-floating point representation of u ² to set the “b” register of the modulo-5 Cuccaro adder (i.e., the register that gets overwritten with addition result). Here, a shift is used to account for the difference in exponent of u and u ². This part of the 0011a is not affected by the partial reduction of workspace qubits. The second part of 0011a set −u (for d v 1 | d v 0 = 00 ) or u (for d v 1 | d v 0 = 11 ) into “a” register of modulo-5 adder. Figure 18 shows the quantum-circuit implementation of this second step before a reduction of workspace qubits is performed. It can be seen that in a partial reduction of the workspace qubits, one or more of the qubits in the “a” register need to be removed. For the example of the reduction to 26 qubits, the qubits a 4 | a 3 | a 2 will be eliminated. For the reduction to 25 qubits, qubits a 4 | a 3 | a 2 | a 1 will be involved as indicated with the red box in Figure 18.

The quantum circuit illustrated in Figure 18 first “copies” the mantissa qubits (including the “hidden” qubit) into the “a” register for d v 1 | d v 0 = 00 and for d v 1 | d v 0 . To perform −u + u ², the second part of the circuit transforms the qubits to 2’s complement for d v 1 | d v 0 = 00 . Without this change to 2’s complement formulation, i.e., for d v 1 | d v 0 = 11 the reduction of the workspace for the modulo-5 adder can be performed using the approach previously shown for the 4-qubit adders in the evaluation of u ².

For d v 1 | d v 0 = 00 , in most cases, the need arises to use different computational kernels (derived for different choice of the reduced qubits), depending of the state of the remaining mantissa qubits of u. The following three examples illustrate this:

(1) Reduction by 2 qubits to 26-qubit circuit, m u 2 | m u 1 = 0 . For normalized inputs, we then have q u 3 | q u 2 = 10 in the representation without hidden qubit. For d v 1 | d v 0 = 00 : “-u” into modulo-5 adder, now 4 cases need to be considered:

For the first two examples, Figure 19 shows the quantum-circuit implementation of the reduced MA5 operation for d v 1 | d v 0 = 00 . The circuits were derived by first creating two separate kernels, e.g. for a 4 | a 3 | a 2 = 110 and for a 4 | a 3 | a 2 = 101 in the reduction-by-2 qubits example. Then, the differences between the kernels are performed conditional and ancilla qubit a n c in the combined circuits shown.

Note: The quantum circuit reductions shown here were performed manually and certainly pose major challenges for automation by circuit compilation methods.

For the D1Q3 model 2 qubits d v 1 | d v 0 are required independent of N _M and N _E . Using the hidden-qubit approach, the u-velocity is defined in terms of exponent, sign and mantissa using N _E + 1 + N _M − 1 = N _E + N _M qubits. For the output g ⃗ e q , similarly N _E + N _M are required since the same floating-point format is used. A single qubit c u t is used to identify cases with truncation to 0 of u ². The temporary storage of u ² in floating-point representation requires N _E + N _M − 1 qubits since a “sign” qubit is not required. The remaining qubits represent the “workspace” of the algorithm. Two computational steps are performed (as well as their un-computation) within this space: the shift-and-add based evaluation of u ² and the “signed” addition ±u + u ² for directions d v 1 | d v 0 = 00 and d v 1 | d v 0 = 00 . For the example N _M = 4 this addition is performed using a 5-qubit modulo adder. The modulo adder requires 2N _M + 1 qubits. The u ² evaluation a 2N _M results register, N _M qubits for unsigned velocity input, 1 “carry” qubit as well as 1 ancilla qubit. In total, the u ² evaluation therefore requires 3N _M + 2 qubits. This exceeds the requirement of the modulo adder and therefore the workspace needs a total of 3N _M + 2 qubits. The total space complexity of the original quantum circuit therefore is: 2 + 2 N E + N M + 1 + N E + N M − 1 + 3 N M + 2 = 4 + 3 N E + 6 N M

After the first reduction step introduced in Section 9, the following qubits were removed from the original circuit: 2 qubits d v 1 | d v 0 , N _E + N _M (input “u” in signed floating-point format) as well as the single c u t qubit. Therefore, the total space complexity for the quantum circuit following this first reduction step is therefore: 4 + 3 N E + 6 N M − 2 + N E + N M + 1 = 1 + 2 N E + 5 N M

The further reduction detailed in Section 9.1, one or more of the N _M qubits defining the unsigned velocity input into the u ² evaluation were eliminated. In its most aggressive form, this reduction step can eliminate all N _M qubits defining the unsigned velocity input, so that the space complexity reduces further to: 1 + 2 N E + 5 N M − N M = 1 + 2 N E + 4 N M

To illustrate the complexity for different choices of N _M , Table 3 summarized the required number of qubits for N _E = 3 and increasing N _M for the original quantum circuits as well as the reduction steps 1 and 2.

Figure 20 shows the memory required to store the qubit information as a pair of 32-bit floating point numbers. For reference, the red lines show the FPGA board memory and 1 TB, 1 PB, 1 EB and 1 ZB. To put this into context: the UK’s largest supercomputer, Archer, comprises 4,920 nodes with each 64 GB of memory, so the total memory is still less than 1 PB.

We evaluated our approach on a single-node FPGA system. Our host system has a dual Intel Xeon E5-2609 V2 2.5 GHz processor and 64 GB RAM (DDR3, 1.6 GHz). This system hosts a Nallatech PCIe-385N A7 FPGA board with 8 GB RAM (DDR3) connected through a PCIe 2 connection. The system runs Scientific Linux 6.8 and we used the Intel SDK for OpenCL Version 17.1 to communicate with and program the FPGA device. The CPU used for evaluation is the same as the FPGA host. The host code is compiled using G++ (GCC version 4.7.2). To evaluate the CPU approach, we run the same OpenCL kernel with the same host code. Our GPU system has access to an NVIDIA GK110B (GeForce GTX TITAN Black), with access to 6 GB of VRAM, which is used for the GPU evaluation target. The GPU system compiles the C++ host code with G++ (GCC version 5.4.0).

The above results show the potential of our approach of circuit transformations to reduce the number of qubits in exchange for a higher number of control bits. Table 6 shows a summary of the runtimes used to infer the following performance and power consumption results. With the current status of our work, in absolute terms, the performance of the simulation on the FPGA is about half as fast as on the CPU host. The GPU benefits from very high memory bandwidth and thousands of compute threads, and thus it shines for this particular approach of simulation (QWM with a full state vector), beating the FPGA by a factor of 20 for high-qubit circuits.

However, as explained before, to scale to larger numbers of qubits, it is necessary to move to clusters of FPGAs. For supercomputer clusters, the dominant factor is the energy consumption and therefore we need to consider the performance per Watt. From that perspective, the picture is quite different: the measured power consumption of the FPGA board is 25 W (Segal et al., 2014); the host CPU consumes 160 W (not including RAM power consumption). So already the FPGA simulation has more than 3× better performance per Watt than the CPU. Our evaluation GPU consumes 250 W of power, meaning it is about 2× better than the FPGA in terms of performance per Watt.

Furthermore, there is a lot of scope for improvement of the FPGA performance. We used a baseline over-PCIe architecture without pipeline control to obtain the presented results. For parallelisation on the device, the architecture currently relies on the dynamic scheduler of the OpenCL NDRange construct. It is unclear how well this pipelines our design and we estimate the architecture would likely benefit from tighter control over pipelining. This would involve using OpenCL channels to explicitly control the flow of the amplitudes through the design. In addition, we are currently not applying any caching techniques in our design, several of which can be applied. Each of our kernels currently reads two values at a time from the distributed RAM of the FPGA, performs some computation on them, and writes them back. We can make better utilization of the DRAM’s bandwidth by reading the required amplitudes in consecutive blocks. Due to the exponential strides between the amplitudes, it may not always be possible to fill the cache in one memory read with matching amplitudes. However, we can show that this is possible with two contiguous memory reads. Furthermore, another optimization technique which would be particularly useful for circuits like the ones presented here is state vector compression. For circuits whose intermediary state vectors present repetition in their amplitude values, it is not always necessary to store the full state vector in memory. Instead, memory access is replaced by computation over a compressed state vector space; a pattern very suitable for application on an FPGA.

With these improvements, which are the focus of our research at the moment, it is expected that the FPGA will outperform the CPU even in absolute terms, so that we are confident that the FPGA simulation can achieve an order of magnitude better performance per Watt. We are also aiming to show that an optimized FPGA architecture can outperform a GPU in performance per Watt and total energy consumption.