Fast Arithmetic Hardware Library For RLWE-Based Homomorphic Encryption

In this work, we propose an open-source, first-of-its-kind, arithmetic hardware library with a focus on accelerating the arithmetic operations involved in Ring Learning with Error (RLWE)-based somewhat homomorphic encryption (SHE). We design and implement a hardware accelerator consisting of submodules like Residue Number System (RNS), Chinese Remainder Theorem (CRT), NTT-based polynomial multiplication, modulo inverse, modulo reduction, and all the other polynomial and scalar operations involved in SHE. For all of these operations, wherever possible, we include a hardware-cost efficient serial and a fast parallel implementation in the library. A modular and parameterized design approach helps in easy customization and also provides flexibility to extend these operations for use in most homomorphic encryption applications that fit well into emerging FPGA-equipped cloud architectures. Using the submodules from the library, we prototype a hardware accelerator on FPGA. The evaluation of this hardware accelerator shows a speed up of approximately 4200x and 2950x to evaluate a homomorphic multiplication and addition respectively when compared to an existing software implementation.


Introduction
As the internet becomes easily accessible, almost all electronic devices collect enormous amounts of private and sensitive data from routine activities. These electronic devices may be as small as wearable electronics [PJ03], like a smart watch collecting personal health information or a cell phone collecting location information, or may be as large as an IoT-based smart home [BFZY11] [ KJBB16] collecting routine information like room temperature, door status (open or closed), smart meter reading, and other such details. These electronic devices have limited power, storage, and compute capabilities and often need external support to process the collected information. Cloud computing [MG + 11] [Hay08] provides a convenient means not only to store the collected information but also 2 Fast Arithmetic Hardware Library For RLWE-Based Homomorphic Encryption to apply various compute functions to this stored data. The processed information can be used easily for various purposes, like machine learning predictions.
As cloud computing services become readily available and affordable, many industrial sectors have begun to use cloud services instead of setting up their own infrastructure. Sectors like automotive production, education, finance, banking, health care, manufacturing, and many more leverage cloud services for some or all of their storage and computing needs. This, in turn, leads to the storage of a lot of sensitive data in the cloud. Hence, while cloud computing provides the convenience of sharing resources and compute capabilities for individuals and business owners equally, it brings its own challenges in maintaining data privacy [PH10] [KIA11]. A cloud owner has access to all of the private data pertaining to the clients and can also observe the computations being carried out on this private data. Moreover, cloud services are shared among many clients; therefore, even if a client may assume that the cloud service provider is honest and will ensure that there is no data breach in their environment, the chances of data leakage remain high, due to the shared storage space or compute node on the cloud.
Homomorphic encryption [G + 09] is a ground-breaking technique to enable secure private cloud storage and computation services. Homomorphic encryption allows evaluating functions on encrypted data to generate an encrypted result. This result, when decrypted, matches the result of the same operations performed on the unencrypted data. Thus, a data owner can encrypt the data and then send it to cloud for processing. The cloud running the homomorphic encryption based services will perform computations on the encrypted data and send the results back to the data owner. The data owner, having access to the private key, performs the decryption and obtains the result. The cloud does not have access to the private key or the plain data, and, hence, the security concerns related to private data processing on the cloud can be mitigated. An illustrative scenario is shown in Figure 1. 09] provided a framework to make fully homomorphic encryption feasible, and almost a decade's work has now made it practical [NLV11]. While homomorphic encryption has become realistic, it still remains several magnitudes too slow, making it expensive and resource intensive. There are no existing homomorphic encryption schemes with performance levels that would allow large-scale practical usage. Substantial efforts have been put forward to develop full-fledged software libraries for homomorphic encryption. Such libraries include SEAL [CLP17], Palisade [Tec19], cuHE [DS15], HElib [HS14], NFLLib [AMBG + 16], Lattigo [LT219], and HEAAN [Kim18]. All of these libraries are based on the RLWE-based encryption scheme, and they generally implement Brakerski-Gentry-Vaikuntanathan (BGV) [BGV14], Fan-Vercauteren (FV) [FV12], and Cheon-Kim-Kim-Song (CKKS) [CKKS17] homomorphic encryption schemes with very similar parameters. Rashmi Agrawal, Lake Bu, Alan Ehret and Michel A. Kinsy 3 Although the software implementations are impressive, they are still incapable of gaining the required performance, as they are limited by the underlying hardware. For example, Gentry et al. [GHS12], in their homomorphic evaluation of an AES circuit, reported approximately 48 hours of execution time on an Intel Xeon CPU running at 2.0GHz. Even their parallel SIMD style implementation took around 40 minutes per block to evaluate 54 AES blocks. A modern Intel Xeon CPU takes about 20ns to perform a regular AES encryption block, hence it is evident that homomorphic evaluation of an AES block is about 1.2 × 10 11 times slower than a regular evaluation. Similarly, logistic regression, a popular machine learning tool, is often used to make predictions using client's private data in the cloud. A software-based homomorphic logistic regression prediction takes about 1.6 hour while a regular logistic prediction takes about 95ns.
If homomorphic encryption's full potential and power can be unleashed by realizing the required performance levels, it will make cloud computing more reliable via enhanced trust of service providers and their mechanisms for protecting users' data. Hence, there is a need to accelerate the homomorphic encryption operation directly on the hardware to achieve maximum throughput with a low latency. With this in mind, we propose an arithmetic hardware library that includes the major arithmetic operations involved in homomorphic encryption. A hardware accelerator designed using the modules from this library can reduce the computational time for HE operations. To lower the power usage and improve performance, new cloud architectures integrate FPGAs to offload and accelerate compute tasks such as deep learning, encryption, and video conversion. The FPGA-based design and optimization approach introduced in this work fits into this class of FPGA-equipped cloud architectures.
The key contributions of the work are as follows: • A fast and hardware-cost efficient hardware arithmetic library to individually accelerate all operations within homomorphic encryption. A speedup of 4200× and 2950× is observed to evaluate homomorphic multiplication and addition respectively.
• An open-source, FPGA-board agnostic, parameterized design implementation of the modules to provide flexibility to adjust parameters so as to meet the desired security levels, hardware cost and multiplication depth.
• A modular and hierarchical implementation of a hardware accelerator using the modules of the proposed arithmetic library to demonstrate the speedup achievable in hardware.
The rest of the paper is organized as follows. In Section 2, we briefly present the underlying scheme and discuss the required arithmetic operations. Section 3 introduces these arithmetic operations and their efficient implementation. In Section 4, we evaluate the associated hardware cost and latency and then conclude the paper in Section 5 along with future work.

Homomorphic Encryption
To present the hardware library, we first start by introducing the underlying RLWE-based homomorphic encryption scheme. For this purpose, we chose the Fan-Vercauteren (FV) [FV12] scheme, as it has more controlled noise growth while performing homomorphic operations when compared to approaches like the BGV (Brakerski-Gentry-Vaikuntanathan) scheme [BGV14]. Moreover, Costache et al. [CLP] presented results showing FV scheme outperforming BGV for large plaintext moduli.

FV Scheme
The FV scheme operates in the ring The plaintext m is chosen in the ring R t for some small t and a ciphertext consists of only one element in the ring R Q for a large integer Q. The security of the scheme is governed by the degree of this polynomial f (x) and the size of Q.
The secret key, s k , is sampled from the ring R 2 or a Gaussian distribution, χ k . The public key, a, is sampled from the ring R Q and the error vector, e, is sampled from a second Gaussian distribution, χ err . The other public key, b, is computed as follows: RQ (1) Encryption of the plaintext message yields a pair of ciphertexts as follows: The homomorphic addition operation adds two such pairs of ciphertexts: After the addition operation, decryption is done as: The homomorphic multiplication operation multiplies two such pairs of ciphertexts using the following equations: Decryption, after the multiplication operation, using the secret key is computed as: Since after the multiplication operation a degree 2 ciphertext is obtained, to continue further multiplication operations this degree 2 ciphertext needs to be reduced to a degree 1 ciphertext. In the FV scheme, this is achieved by performing a relinearisation operation. The scheme facilitates two different approaches for performing relinearisation. Relinearisation version 1 operation consists of generating the relinearisation key, decomposing c 2 to limit the noise explosion, and then conversion to a degree 1 ciphertext by using generated relinearisation keys and the decomposed c 2 . The relinearisation keys are generated as follows: Here, T is independent of t and ℓ = ⌊log T (Q)⌋. Decomposition of c 2 involves rewriting c 2 in base T and can be computed using the following equation: Next the relinearisation operation can be performed as follows: Rashmi Agrawal, Lake Bu, Alan Ehret and Michel A. Kinsy 5 Since, we have obtained a degree 1 ciphertext after the relinearisation operation, the decryption can be performed without using the s 2 k term, as in equation 6. Therefore, the decryption operation simplifies to the equation: Note that the choice of T will determine the size of relinearisation keys and the noise growth during the relinearisation operation. The larger the value of T , the smaller the relinearisation keys will be, but the noise introduced by relinearisation will be higher. And the smaller the value of T , the larger the relinearisation keys will be, with smaller noise introduction. So the value of T must be picked in a balanced way.
Relinearisation version 2 is a modified form of modulus switching and hence requires choosing a second modulus p such that p ≥ Q 3 for small enough error samples. Now, the relinearisation keys can be generated as follows: Here, a ∈ R p·Q and e ← − χ ′ err . We can perform the relinearisation using the following computation: Once the c 2 component is removed, we can perform decryption using equation 10.

Required Operations
Our proposed arithmetic library includes highly optimized hardware-based implementations of Residue Number System (RNS), Chinese Remainder Theorem (CRT), modulo inverse, fast polynomial multiplication using Number Theoretic Transform(NTT), polynomial addition, modulo reduction, Gaussian noise sampler and relinearisation operations. While implementing these operations for our arithmetic library, the design choices are highly motivated by the parameter selection. This is because an RLWE-based encryption scheme requires adding a small noise vector to obfuscate the plaintext message as shown previously in equations 1 and 2. While performing homomorphic addition and multiplication, the noise present in the ciphertexts gets doubled and squared respectively. Due to this noise growth, the ring R Q , along with the degree of the polynomial needs to be large, so as to compute a circuit of certain depth and still enable successful decryption of the result. Hence, the parameter Q with the degree of the polynomial f (x) needs to be large. The operations on this large parameter set not only increase the hardware cost but also slow down the homomorphic encryption.

Residue Number System
Chinese Remainder Theorem
We use the concepts of modular arithmetic to speed up HE computations. The underlying FV scheme does not restrict Q to a prime number; instead, Q can be the product of small primes. When working modulo a product of numbers, say Q = q 1 × q 2 × · · · × q k , Residue Number System (RNS) helps reduce the coefficients in each of the R qi and Chinese Remainder Theorem (CRT) lets us work in each modulus q i separately. Since the computational cost is proportional to the size of operands, this is faster than working in the full modulus Q. Moreover, breaking down the coefficients into smaller integers using Rashmi Agrawal, Lake Bu, Alan Ehret and Michel A. Kinsy 7 RNS also limits the noise expansion. Figure 2 illustrates the sequence of operations that will be required in HE using the FV scheme, for Q = q 1 × q 2 × q 3 × q 4 .
In our implementation of the arithmetic library, we consider Q to be 1200 bits and degree of the polynomial, n as 1024 with a 128-bit security level. With these values of Q and n, we can evaluate a binary circuit of depth 56 using somewhat homomorphic encryption. The RNS module will take a 1200-bit wide integer coefficient as an input, perform x i = x mod q i , and thus break x into 40 small integers of 30 bits each. This enables us to set up 40 pipelines to perform 40 operations in parallel, providing the required performance boost. Once the homomorphic add or multiplication operation is done on these small integers, CRT can combine them to map back to the original 1200-bit width. Note that the bit width selection for the small integers is a design decision one can make based on available resources. Our implementations take q i as a parameter, which facilitates different bit width selections. Figure 3 shows all the core building blocks required to perform somewhat homomorphic encryption (SHE) addition and multiplication operations using FV scheme. The client-side building blocks include key generation, relinearisation (versions 1 and 2) key generation, encryption, decryption (for both degree-1 and degree-2 ciphertexts), residue number system (RNS) along with modular reduction, and Chinese remainder theorem (CRT) along with modulo inversion. The cloud provider has blocks to perform homomorphic addition, homomorphic multiplication, and relinearisation (versions 1 and 2). While RNS and CRT are standalone modules, the rest of the main blocks share the following submodules: polynomial multiplication, polynomial addition, modular reduction, scalar multiplication, scalar division to nearest binary integer, noise sampler, and true random number generator. Certain operations like decompositions, powers-of-2 computations, divide and round operations are specifically required for the purpose of relinearisation, which we will discuss in detail in Section 3.10.

Arithmetic Hardware Library For HE
We will start with discussing the standalone modules first (i.e., RNS and CRT) along with the supplemental operations to these modules. Then, we will describe other arithmetic operations shared between all the core building blocks, along with their design implementations. All the operations are customized for hardware-based implementation, and we include both hardware-cost efficient serial and fast parallel implementations in the library. It is worth noting that, in all the algorithms, Q will denote the large integer and q i or q will denote a prime factor of Q.

Residue Number System
is a mathematical way of representing an integer by its value modulo a set of k integers q 1 , q 2 , q 3 , . . . , q k , called the moduli, which generally should be pairwise coprime. An integer, x, can be represented in the residue number system by a set of its remainders x 1 , x 2 , x 3 , . . . , x k under Euclidean division by the respective moduli. That is, The serial implementation of RNS is shown in Figure 4. Each 1200-bit coefficient is modulo reduced by a q i and stored in the respective BRAM. For k moduli, it takes k cycles to perform all the computations. When modulo reductions are performed in parallel, all the computations can be performed in a single clock cycle instead. The parallel implementation is shown in Figure 4. Since the mod operation is the key operation in RNS, we next optimize the modulo reduction operation. This will allow us to reduce the hardware cost substantially.

Modular reduction
Modular reduction is not only at the core of many asymmetric cryptosystems, it is the most performed operation in encryption schemes based on R-LWE. This is because, in RLWE, all the operations are required to be performed over large finite rings. The function of the modular reduction operation is to compute the remainder of an integer division.
Mathematically, it is written as r = a (mod q).
While the modular reduction operation sounds relatively simple, the division of two large integers is very costly. Moreover, the moduli in RLWE-based schemes are prime numbers and not power-of-2 numbers, which makes the operation non-trivial. Therefore, the hardware implementation of the modulo operation is quite expensive. For example, the use of the inbuilt modulo operator, % in Verilog for 30-bit operands utilizes about 800 LUTs, and when there are many such modular operations involved, the hardware cost quickly adds up. Hence, optimization of modulo operation can lead to significant hardware cost reductions.
One well-known modular reduction optimization algorithm is Barrett reduction [Bar87]. It is preferred over Montgomery reduction [Mon85], as it operates on the given integer number directly, while Montgomery reduction requires numbers to be converted into and out of Montgomery form, which is expensive in itself. We will discuss the Barrett reduction next, and then later, we propose some modifications to the existing Barrett reduction algorithm to reduce the hardware cost further.

Barrett reduction
The Barrett reduction algorithm was introduced by P. D. Barrett [Bar87] to optimize the modular reduction operation by replacing divisions with multiplications, so as to avoid the slowness of long divisions. The key idea behind the Barrett reduction is to precompute a factor using division for a given prime modulus, q, and thereafter, the computations only involve multiplications, subtractions and shift operations. These operations are faster than the division operation. The Barrett reduction algorithm steps are shown in Algorithm 1 and works as described below.
Since the modulus q is known in advance and the factor r depends only on this modulus, it can be precomputed and stored. Then, the reduction function requires only computing the remainder value, t. While computing the t value, a division operation by 4 k , being power-of-2, can be performed using the right-shift operation. Hence, the entire compu-hmi Agrawal, Lake Bu, Alan Ehret and Michel A. Kinsy 9 Algorithm 1: Barrett Reduction Algorithm , and q is a prime number.
2 3 Precomputation: /* Choose k ∈ N such that 2 k > q */ 4 assign k = ⌈log 2 q⌉; /* Compute the factor */ 5 assign r = ⌊ 4 k q ⌋; 6 7 Reduction Function: Return t as a mod q . 11 end 12 else 13 Return t − q as a mod q . 14 end entire computation reduces to just two multiplications, one right-shift, and one subtraction operation. Furthermore, the computation is performed in one step and, thus, is performed in constant time. The hardware implementation circuit is as shown in Figure ??
For our hardware-based Barrett reduction implementation, we specifically included some additional optimizations. One such optimization is a careful bit width analysis. Say that the modulus requires exactly bits, then the product ar fits in bits. We also observed that the computed values do not need more than + 1 bits. The advantage of this observation is that we can safely ignore the upper bits of the product. This in turn reduces the size of the registers to + 1 bits while performing computations.

Modified Barrett reduction
enplaugh et al. [ ] introduced an iterative folding method as a modification to the Barrett reduction method. This method not only reduces the number of required multiplications via an increased number of precomputations, but also reduces the bit width of the operations performed. We modify their proposed approach and propose Algorithm ?? which computes modulo reduction in a single fold.
When compared to Barrett reduction, the proposed algorithm precomputes with half the bit width and with one third the bit width. This significantly reduces the multiplication bit-width while performing actual computations after the coefficient integers are known. Moreover, we were able to get rid of an additional check, t < q, which is tation reduces to just two multiplications, one right-shift, and one subtraction operation. Furthermore, the computation is performed in one step and, thus, is performed in constant time. The hardware implementation circuit is as shown in Figure 5.
For our hardware-based Barrett reduction implementation, we specifically included some additional optimizations. One such optimization is a careful bit width analysis. Say that the modulus requires exactly k bits, then the product ⌊ ar 4 k ⌋q fits in 2k bits. We also observed that the computed values t do not need more than m + 1 bits. The advantage of this observation is that we can safely ignore the upper m − 1 bits of the product. This in turn reduces the size of the registers to m + 1 bits while performing computations.

Modified Barrett reduction
Hasenplaugh et al. [HGG07] introduced an iterative folding method as a modification to the Barrett reduction method. This method not only reduces the number of required multiplications via an increased number of precomputations, but also reduces the bit width of the operations performed. We modify their proposed approach and propose Algorithm 2, which computes modulo reduction in a single fold. When compared to Barrett reduction, the proposed algorithm precomputes k with half the bit width and r with one third the bit width. This significantly reduces the multiplication bit-width while performing actual computations after the coefficient integers /* Choose k ∈ N such that 2 3k > q */ 4 assign k = ⌊ log 2 q 2 ⌋; /* Compute the factor */ 5 assign r = ⌈ 2 3k q ⌉; 6 7 Reduction Function: required in case of Barrett reduction. Hence, the result for modulo reduction is available with minimal operations using minimal bit width. The modified Barrett reduction's implementation is shown in Figure ??

Polynomial Multiplication using NTT
Polynomial multiplication is the most performed operation in homomorphic encryption and has the highest implementation complexity. Therefore, the latency of the polynomial multiplication module will govern the efficiency of the entire implementation. Hence, it is critical to design an efficient polynomial multiplication module. A conventional approach to implement a polynomial multiplier is to use convolution method. However, this approach is expensive to implement in hardware, as it requires performing multiplications for a degree polynomial. This complexity can be reduced to log multiplications instead by using NTT combined with negative wrapped convolution to perform polynomial multiplication. We leverage the NTT-based multiplication algorithm proposed by Chen et al. [ ] in our implementation. The steps involved in this algorithm are described in Algorithm ??

Algorithm 3: based Polynomial Multiplier
Setup:Let · · · , a and · · · , b } ∈ be two polynomials of length (with coefficients), where ) = + 1 is an irreducible polynomial with a power of 2, and mod is a large prime number). Let be the -th root of unity, the inverse of mod , and the inverse of are known. Moreover, we were able to get rid of an additional check, t < q, which is required in case of Barrett reduction. Hence, the result for modulo reduction is available with minimal operations using minimal bit width. The modified Barrett reduction's implementation is shown in Figure 5.

Polynomial Multiplication using NTT
Polynomial multiplication is the most performed operation in homomorphic encryption and has the highest implementation complexity. Therefore, the latency of the polynomial multiplication module will govern the efficiency of the entire implementation. Hence, it is critical to design an efficient polynomial multiplication module. A conventional approach to implement a polynomial multiplier is to use convolution method. However, this approach is expensive to implement in hardware, as it requires performing O(n 2 ) multiplications for a degree n polynomial. This complexity can be reduced to O(n log 2 n) multiplications instead by using NTT combined with negative wrapped convolution to perform polynomial multiplication. We leverage the NTT-based multiplication algorithm proposed by Chen et al. [CMV + 15] in our implementation. The steps involved in this algorithm are described in Algorithm 3.

Algorithm 2: Modified Barrett Reduction Algorithm
Let , and is a prime number.

Precomputation:
/* Choose such that > q */ gn /* Compute the factor */ gn uction Function: compute ⌊⌊ rn tq as a mod q required in case of Barrett reduction. Hence, the result for modulo reduction is available with minimal operations using minimal bit width. The modified Barrett reduction's implementation is shown in Figure ??

Polynomial Multiplication using NTT
Polynomial multiplication is the most performed operation in homomorphic encryption and has the highest implementation complexity. Therefore, the latency of the polynomial multiplication module will govern the efficiency of the entire implementation. Hence, it is critical to design an efficient polynomial multiplication module. A conventional approach to implement a polynomial multiplier is to use convolution method. However, this approach is expensive to implement in hardware, as it requires performing multiplications for a degree polynomial. This complexity can be reduced to log multiplications instead by using NTT combined with negative wrapped convolution to perform polynomial multiplication. We leverage the NTT-based multiplication algorithm proposed by Chen et al. [ ] in our implementation. The steps involved in this algorithm are described in Algorithm ??

Number Theoretic Transform
A generalization of the Fast Fourier Transform (FFT) over a finite ring R q = R/ q = Z q [x]/ f (x) is represented by Number Theoretic Transform (NTT). The equation of the NTT is as follows: where ω is the n th root of unity in the corresponding polynomial field and for a ring R q , where q is a prime number, the n th root of unity ω must satisfy two conditions: 1. ω n = 1 mod q, 2. The period of ω i for i ∈ {0, 1, 2, · · · n − 1} is exactly n.
One of the efficient ways to compute ω is by using the following approach: 1. First compute the primitive root of q, which must satisfy: • α q−1 = 1 mod q • The period of α i for i ∈ {0, 1, 2, · · · q − 1} is exactly q − 1.
2. And since ω n ≡ α q−1 mod q, we can compute: ω = α (q−1)/n mod q 3. As a final step, verify that this ω meets both the conditions mentioned above.
Applying inverse NTT (iNTT) is straight forward and can be performed using the existing NTT module by replacing ω with ω −1 , where ω −1 = ω n−1 mod q. iNTT computation also requires computing the inverse of n, which can be computed as n −1 · n = 1 mod q.
Although there exist hardware implementations of NTT, they are quite expensive because of the way they compute the indices of the points and the corresponding w i . Investing a large number of multiplications and divisions for these computations may not be an issue with software implementation, however they lead to a higher resource consumption in the hardware counterpart. Therefore, in our implementation of the NTT algorithm, Algorithm 3.3.1, we perform the indices computation using only shift and xor operations. The benefit of doing so is that the shift and xor are not only inexpensive to implement but they conveniently replace the large multiplication and division circuits. By leveraging this highly optimized NTT implementation, we implement the fast polynomial multiplication algorithm (Algorithm 4) very efficiently. Figure 6 shows a high level circuit for polynomial multiplication and the operations within the NTT block.

Polynomial Addition
Polynomial addition is the second most frequently used operation after polynomial multiplication. The schematic of the hardware implementation for polynomial addition is shown in Figure 7. The implementation performs a component-wise addition operation on the coefficients of the polynomial. Note that the results are wrapped either within small modulus q or large modulus Q depending on which main module is utilizing this submodule to perform polynomial addition.

Scalar Multiplication
Since the message space is binary, a conditional assignment operator can be used to implement the scalar multiplication operation. As shown in Figure 7, m, the plaintext message, is an n-bit vector. Thus, computing tm essentially requires choosing t or 0 according to each bit of m. Since we avoid performing actual multiplication operations, hardware cost is greatly reduced.

Fast Arithmetic Hardware Library For RLWE-Based Homomorphic Encryption
Algorithm 4: Number-Theoretic Transform 1 Let a = {a 0 , · · · , a n−1 } ∈ Zq[x]/ f (x) where n is a power of 2, and q ≡ 1 mod 2n is a prime number. Let ω be the n-th root of unity for q. Let ω −1 be the inverse of ω. Precompute ω i and ω −i for i ∈ {0, 1, · · · , n/2 − 1}, and store them in the array element ω[i] and iω[i] respectively. 2 Let a be swapped to A so that A[j] = a [jrev], where jrev is a binary vector bit-reversed from j.
Let i, Stage both be initialized to 0.

Scalar Division to the Nearest Binary Integer
The scheme parameter t = ⌊ Q 2 ⌋ value is already published, and hence, it is known. From the equations 4 and 6 in the scheme, if we denote u = (c 0 + s · c 1 ), then to decrypt the message m correctly, all we need to do is to compute m = ⌈ u t ⌋. The nearest binary integer equivalent of ⌈ u t ⌋ can be computed by measuring the distance between u and t as (Absolute(u − t) < t 2 ) ? 1 : 0. If this distance is larger than half of t, it indicates that u and t are far from each other, and thus, the nearest integer of the quotient u t must be 0. But if this distance is less than half of t, then the nearest integer of the quotient u t must be 1. Thus, we implement the scalar division hardware circuit as shown in Figure 7, without using any hardware division circuit.

Chinese Remainder Theorem
The Chinese remainder theorem, CRT [KI07], states that if we know the residue of an integer, a modulo two primes q 1 and q 2 , it is possible to reconstruct < a > q1q2 as follows. Let < a > q1 = a 1 and < a > q2 = a 2 , then the value of a (mod Q), where Q = q 1 · q 2 , can be found by where t 1 is the multiplicative inverse of q 1 (mod q) 2 and t 2 is the multiplicative inverse of q 2 (mod q) 1 . This is feasible as the inverses t 1 and t 2 always exist, since q 1 and q 2 are coprime. Mathematically, < a > q1q2 can also be represented by a set of congruent equations as follows: Using the CRT, we combine all the small integers back into one large integer, so as to generate the required final result. A naive approach to implement CRT would compute pairwise multiplicative inverse or modulo inverse for two given moduli and then use equation 14 to merge values of a 1 and a 2 to get a. This process can be carried out recursively until the final coefficient value is obtained.
14 Fast Arithmetic Hardware Library For RLWE-Based Homomorphic Encryption equations as follows: (mod (mod Using the CRT, we combine all the small integers back into one large integer, so as to generate the required final result. A naive approach to implement CRT would compute pairwise multiplicative inverse or modulo inverse for two given moduli and then use equation ?? to merge values of and to get . This process can be carried out recursively until the final coefficient value is obtained. Algorithm 5: Chinese Remainder Theorem LUT-based Algorithm 1 Let a 1 , · · · , a k ∈ Zq i [x]/ f (x) (1 ≤ i ≤ n), q 1 , · · · , q k are pairwise relatively prime and large modulus Q = q 1 × · · · × q k 2 3 Precomputation: The problem with this approach is that computing the multiplicative inverse at runtime increases latency. Therefore, a better approach is to precompute the modulo inverse values since all the moduli are known in advance. And then the actual computation reduces to just a single step, which can be performed in one clock cycle. The precomputation and computation steps involved are shown in Algorithm 5. We call this approach LUT-based, since the precomputated values are stored in LUTs; its hardware implementation circuit is shown in Figure 8. The hardware cost for LUT-based CRT can be further optimized by breaking down the single step multiplication and addition operation into various steps. This will enable the reuse of multipliers and adders during these steps. We leave this as future work for now.

Modulo Inverse
A modulo inverse or multiplicative inverse is the main computation involved in CRT. Moreover, while working with homomorphic encryption, due to large parameters, the amount of storage available can be a concern. Thus, instead of using LUT-based CRT, we may need to use a regular CRT implementation with the modulo inverse computed on the fly. The multiplicative inverse of a (mod q) exists if and only if a and q are relatively prime (i.e., if gcd(a, q) = 1). Given two integers a and q, the modulo inverse is defined by an integer p such that a · p ≡ 1 (mod q) Here, the value of p should be in 0, 1, 2, . . . q − 1, i.e., in the range of integer modulo q. In our case, we need to compute the multiplicative inverse between pairwise moduli, q i . There are two primary algorithms used in computing a modulo inverse. We discuss both algorithms in detail next.

Fermat's Little Theorem
Fermat's little theorem [Vin16] is typically used to simplify the process of modular exponentiation. But since we know q is prime, we can also use Fermats's little theorem to find the modulo inverse. According to this theorem, we can rewrite equation 16 as follows: If we multiply both sides of this equation with a −1 , we get Equation 18 is the only computation carried out by this theorem to get the value of a −1 and this is what is shown in Algorithm 6. The algorithm has a time complexity of O(log 2 q). For our hardware-based implementation, we precompute the power factors from 1 to q − 2 to save the computation cost. This not only speeds up computation but also significantly reduces the hardware cost. The hardware implementation is shown in Figure  9. modulo inverse. According to this theorem, we can rewrite equation ?? as follows: ? is the only computation carried out by this theorem to get the value of and this is what is shown in Algorithm ??. The algorithm has a time complexity of log . For our hardware-based implementation, we precompute the power factors from to to save the computation cost. This not only speeds up computation but also significantly reduces the hardware cost. The hardware implementation is shown in Figure  ??
Algorithm 6: Fermat's Little Theorem 1 Let a ∈ Zq[x]/ f (x) , and q is a prime number.

Extended Euclidean Algorithm
The extended Euclidean algorithm [Vin16] is an extension to the classic Euclidean algorithm that is used for finding the greatest common divisor (GCD). According to this algorithm, if a and q are relatively prime, there exist integers x and y such that ax+qy = 1, and such integers may be found using the Euclidean algorithm. Considering this equation modulo q, it follows that ax = 1; i.e., x = a −1 (mod q). The algorithm used for our implementation is as shown in Algorithm 7. The algorithm works as follows. Given two integers 0 < a < q, using the classic Euclidean algorithm equations, one can compute gcd(a, q) = r j , where r j is the remainder. In the classic Euclidean algorithm, we start by dividing q by a (integer division with remainder), then repeatedly divide the previous divisor by the previous remainder until there is no remainder. The last remainder we divided by is the greatest common divisor. To avoid division operations, the classic Euclidean algorithm equations can be rewritten as follows: . . .
Then, in the last of these equations, r j = r j−2 − r j−1 · x j , replace r j−1 with its expression in terms of r j−3 and r j−2 from the equation immediately above it. Continue this process successively, replacing r j−2 , r j−3 , . . . , until we obtain the final equation r j = ax + qy, with x and y integers. In our special case that gcd(a, b) = 1, the integer equation reads as 1 = ax + qy and therefore we deduce 1 ≡ ax (mod q) so that the residue of Fast Arithmetic Hardware Library For RLWE-Based Homomorphic Encryption x is the multiplicative inverse of a (mod q). The time complexity of this algorithm is O (log(min(a, q))).
Then, in the last of these equations, , replace with its expression in terms of and from the equation immediately above it. Continue this process successively, replacing , r , . . . , until we obtain the final equation ax qy, with and integers. In our special case that gcd a, b) = 1, the integer equation reads as 1 = ax qy and therefore we deduce ax mod so that the residue of is the multiplicative inverse of mod . The time complexity of this algorithm is log a, q Algorithm 7: Extended Euclidean Algorithm , and q is a prime number.
2 3 Initialization: 4 assign tempa = q; 5 assign tempb = a%q; 6 assign m = tempa tempb ; 7 assign x = 0, prevx = 1; 8 assign y = 1, prevy = 0; The hardware implementation is shown in Figure ??. The algorithm is simplified by removing unnecessary variables and computations to make it more suitable for hardware implementation. The implementation is done in an iterative fashion so that the input parameters gradually decrease while keeping the GCD of the parameters unchanged. The hardware implementation is shown in Figure 9. The algorithm is simplified by removing unnecessary variables and computations to make it more suitable for hardware implementation. The implementation is done in an iterative fashion so that the input parameters gradually decrease while keeping the GCD of the parameters unchanged.

Gaussian Noise Sampler
The security of the RLWE-based encryption scheme is governed by small error samples generated from a Gaussian distribution. Hence, a Gaussian noise sampler lies at the core of maintaining the required security level. However, it is critical to select a sampling algorithm with a high sampling efficiency and throughput so that the key generation and the encryption operations, at the client side, still remain efficient. We leverage the implementation of a Ziggurat-based Gaussian noise sampler done by the authors in [ABK]. Due to space constraints we do not present the implementation details here but the interested readers can refer to the actual paper.

Relinearisation
We discuss relinearisation version 1 implementation details first. In version 1, the key generation will reuse most of the existing submodules except for the powers-of-2 computation. This operation is indicated by the PowersOf2 submodule in Figure 10. The values of T i are not precomputed and stored to reduce the memory overhead. Instead, we take the vector s 2 and perform a left shift operation on all of the elements of this vector. The first set of left shift operations should be by 0 bits to indicate 2 0 multiplication, then by 1 bit for 2 1 multiplication, and so on until 2 ℓ multiplications are performed. Note that we generate the relinearisation keys in R q1 , . . . , R q k rather than R Q . This facilitates the routing of the relinearisation keys correctly to the corresponding q i operation pipeline without the need to perform modular reductions. Although we perform k times more operations, these operations are significantly faster. Additionally, the output of this submodule is arranged in such a fashion that the elements having the same index, from key rlk [0], are treated as a single output. A similar output format holds true for rlk[1] as well. This helps in faster indexing of the relinearisation keys while computing the inner product with c 2 . The schematic of PowersOf2 submodule implementation is as shown in Figure 10.
In the relinearisation version 1 module in Figure 10, the decomposition submodule's task represents the ciphertext c 2 at bit level, i.e., converting the coefficients from R q to R T and here T = 2. We know that in hardware bit-level operations can be performed readily, and hence, this operation becomes trivial. Therefore, we do not specifically provide an implementation of this submodule, and it is shown for completeness in Figure 3. Next, we describe the implementation of the inner product submodule. To avoid performing actual multiplication operations, we leverage our scalar multiplication module (Section 3.5) within this submodule, since c 2 is binary. Hence, a conditional operator does the work of multiplying elements of c 2 and relinearisation keys. We just need to use adders to compute the summation to finish the inner product computations. Implementation of this submodule is shown in Figure 10.
We will explore the relinearisation version 2 submodules now. For key generation, we pick the largest q i from the moduli set, compute q 3 i and then the immediate next power of two is set as the value of p. This p is the scaling factor. Since we choose a power of 2 as p, we can simply perform shift left operations to emulate the multiplication of p with s 2 while generating the relinearisation keys. Note that rlk[1] or a is sampled from R p·qi . Additionally, to maintain the required security, the error samples need to be generated from a different noise sampler. Hence, a second instance of the noise sampler is used here with the required parameter settings. The rest of the submodules are as previously discussed.
While performing the relinearisation operation in version 2, the ciphertext needs to be scaled down. This task is accomplished by using the Div&Round submodule shown in Figure 10. As the scale factor p is a power of 2, division operations can be avoided, and shift right operations can be performed instead. Most other existing implementations (both software and hardware), precompute 1 p , round it down, and perform multiplication operations instead of division. There are two disadvantages to this approach. First rounding leads to loss of precision, generating approximate results and magnifying the errors in decryption as the levels of operations increase. Second, even though the expensive division operations are avoided, multiplications are still costly, requiring large multipliers which are not only expensive but also lead to a lower operating frequency. Figure 10 shows the Div&Round submodule circuit.

Performance Evaluation
We evaluate the performance of all the design implementations through synthesis on a Xilinx Zynq-7000 family xc7z020clg400-1 FPGA. The tool used for synthesis is ISE design suite 14.7, with all designs implemented in Verilog 2001. To generate synthesis results, the input coefficient bit width is considered 1200 bits and there are 40 coprime moduli, q, having 30 bits each. The degree of the polynomial is taken as 1024.

Hardware Cost and Latency
We start by discussing the hardware cost and latency of the individual operations listed in Table 1 and 2. The hardware cost depends on the size of each coefficient (either 1200 bits or 30 bits) and the number of portions into which a 1200-bit coefficient is divided. The latency computation factors in the number of portions, k, or the polynomial degree, n as required by the implementation of a module. That is why in Table 2, the latencies are represented as a factor of k or n. The built in Verilog mod operator is very expensive and takes about 800 LUTs to perform a 30-bit modulo reduction. For the same bit width, the classic Barrett reduction reduces the hardware cost by almost 11 times, and our proposed modified Barrett reduction reduces the cost by about 35 times. Since modular reduction is performed very frequently and used by all the modules, using the modified Barrett reduction method substantially reduces the hardware resources required for implementing the entire homomorphic encryption scheme. Note that the latency of the modular reduction module is not shown, as the implementation comprises combinational logic only. We observe that the RNS parallel implementation utilizes about 12 times more LUT slices, when compared to the serial implementation, while the latency of the serial implementation is about 40 times higher than the parallel one. The RNS serial and parallel modified implementations listed in the table use the modified Barrett reduction to perform the modulo reduction operation, instead of the inbuilt Verilog mod operator. The hardware resource utilization is dramatically reduced by this modification, however, the latency remains the same. The regular CRT implementation is about 3 times more expensive than a LUT-based CRT. This difference is because regular CRT spends a lot of hardware resources for computing the multiplicative inverse, while the LUT-based CRT, with all its precomputations, not only requires fewer hardware resources but also performs computations in 9 times fewer clock cycles. While computing multiplicative inverses, Fermat's little theorem facilitates a low hardware cost implementation, with about half the hardware cost as compared to the widely used extended Euclidean method. However, the extended Euclidean method performs computations about 9 times faster. An NTT-based polynomial multiplication cuts down the latency from n 2 to n log n. The polynomial addition, scalar multiplication, and scalar division submodules avoid the usage of modular reduction, multiplication and division operations respectively. Hence, these submodules are implemented using minimal hardware resources and have a low latency. For the rest of the other submodules involved in relinearisation, because of all the optimizations in implementation, we observe a low hardware cost and latency. Table 3 provides the time, in clock cycles, for computing various homomorphic encryption operations. We represent time in clock cycles due to frequency difference between FPGA and general-purpose CPU. As seen in the table, a single HE addition is about 23× faster than a HE multiplication. Moreover, if one has to choose between relinearisation version 1 and 2, then version 1 would be the unanimous choice, as it is about 6× faster than version 2, even when key generation takes almost the same time. It is worth noting that, although the table lists the relinearisation key generation time for both versions, these keys can be precomputed, and hence, the time required for key generation need not be included in the overall time required for HE operations. Based on the parameters that are used in the implementation, we can evaluate a circuit of depth 56, with relinearisation performed after every multiplication operation. Therefore, using the data from Table 3, we can compute the number of clock cycles required for this entire set of operations. We observe that when the circuit evaluation is done with relinearisation version 1, the total cycles required are 5, 031, 984, while evaluation done with relinearisation version 2 takes 10, 194, 651 cycles. Next, we present the speed up obtained by the hardware accelerator designed using the modules in our hardware library in comparison to its software counterpart. For this comparison, we recorded the number of clock cycles required for encryption, HE multiplication, HE addition and decryption in the Palisade software library using same underlying scheme with RNS implementation using same parameters. Table 4 lists the observed speedup. This evaluation assumes we utilize the maximum possible resources available on the FPGA that we used for our evaluation. There is a scope to further enhance the speedup using more hardware resources.

Hardware library vs Software library Speedup
Finally, we evaluate the time required to make a prediction using logistic regression, a common tool used in machine learning for binary classification problems. Making predictions using logistic regression model with x 0 , x 1 , ..., x n features require computing the logistic regression equation, Y = e X /1 + e X , where X = n i=0 b i x i . We first compute X and then use the value in the Remez algorithm (an iterative minimax approximation algorithm [Fra65]) equation, Y (X) = −0.004X 3 +0.197X +0.5 to compute the probability. The use of Remez algorithm equation helps avoid log computation which is required in the logistic regression equation. However, these equations require working over floating point numbers and the FV scheme does not support the encryption of floating point numbers. Hence, we scale them to integers and perform fixed point operations instead. We observe that our hardware accelerator requires 470, 064 clock cycles to perform one logistic regression prediction using homomorphic encryption, providing a speedup of around 2650× over a software-based prediction as mentioned in Table 5.

Conclusion
We presented a fast hardware arithmetic hardware library with a focus on accelerating the key arithmetic operations involved in RLWE-based somewhat homomorphic encryption. For all of these operations, we include a hardware cost efficient serial implementation and a fast parallel implementation in the library. We also presented a modular and hierarchical implementation of a hardware accelerator using the modules of the proposed arithmetic library to demonstrate the speedup achievable in hardware. The parameterized design implementation approach of the modules and the hardware accelerator provides the flexibility to extend use of the modules for other schemes, such as BGV, and the accelerator for many applications, especially in the FPGA-centric cloud computing environment. Evaluation of the implementation shows that a speed up of about 4200× and 2950× for evaluating homomorphic multiplication and addition respectively is achievable in hardware when compared to software implementation.
As future work, we would like to optimize and implement the arithmetic operations involved in bootstrapping as well. The bootstrap operation is one of the key functions in achieving fully homomorphic encryption, but it remains very expensive to perform. Optimizing the bootstrap operation will render it more practical to use. We are also actively working on integrating other RLWE-based homomorphic encryption schemes, like BGV, into our library so as to leverage inherent advantages that these schemes offer. Once we have the required operations and schemes implemented, we will open-source the arithmetic library and FPGA design examples.