Ilia Iliashenko

Vincent Zucca
email: vincent.zucca@univ-perp.fr

Faster homomorphic comparison operations for BGV and BFV

Fully homomorphic encryption (FHE) allows to compute any function on encrypted values. However, in practice, there is no universal FHE scheme that is efficient in all possible use cases. In this work, we show that FHE schemes suitable for arithmetic circuits (e.g. BGV or BFV) have a similar performance as FHE schemes for non-arithmetic circuits (TFHE) in basic comparison tasks such as less-than, maximum and minimum operations. Our implementation of the less-than function in the HElib library is up to 3 times faster than the prior work based on BGV/BFV. It allows to compare a pair of 64-bit integers in 11 milliseconds, sort 64 32-bit integers in 19 seconds and find the minimum of 64 32-bit integers in 9.5 seconds on an average laptop without multithreading.

Introduction

Fully Homomorphic Encryption (FHE) can perform any kind of computations directly on encrypted data. It is therefore a natural candidate for privacypreserving outsourced storage and computation techniques. Since Gentry's breakthrough in 2009 [START_REF] Gentry | Fully homomorphic encryption using ideal lattices[END_REF], FHE has received a worldwide attention which has resulted in numerous improvements. As a result, FHE can now be used in practice in many practical scenarios, e.g. genome analysis [START_REF] Kim | Private Genome Analysis Through Homomorphic Encryption[END_REF], energy forecasting [START_REF] Bos | Privacy-friendly forecasting for the smart grid using homomorphic encryption and the group method of data handling[END_REF], image recognition [START_REF] Bourse | Fast Homomorphic Evaluation of Deep Discretized Neural Networks[END_REF] and secure messaging [START_REF] Angel | PIR with compressed queries and amortized query processing[END_REF]. In addition, FHE is currently going through a standardization process [START_REF] Albrecht | Homomorphic encryption security standard[END_REF].

In practice, homomorphic encryption (HE) schemes can be classified into three main categories:

-The schemes encrypting their input bit-wise meaning that each bit of the input is encrypted into a different ciphertext. From there, the operations are carried over each bit separately. Examples of such schemes include FHEW [START_REF] Ducas | FHEW: Bootstrapping Homomorphic Encryption in Less Than a Second[END_REF] and TFHE [START_REF] Chillotti | Faster Fully Homomorphic Encryption: Bootstrapping in Less Than 0.1 Seconds[END_REF]. These schemes are believed to be the most efficient in practice with relation to the total running time.

-The second category corresponds to word-wise encryption schemes that allow to pack multiple data values into one ciphertext and perform computations on these values in a Single Instruction Multiple Data (SIMD) fashion [START_REF] Smart | Fully Homomorphic SIMD Operations[END_REF].

In particular, encrypted values are packed in different slots such that the operations carried over a single ciphertext are automatically carried over each slot independently. Schemes with these features include BGV [START_REF] Brakerski | Leveled) Fully Homomorphic Encryption without Bootstrapping[END_REF] and BFV [START_REF] Brakerski | Fully homomorphic encryption without modulus switching from classical GapSVP[END_REF][START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF]. Although homomorphic operations in these schemes are less efficient than for bit-wise encryption schemes, their running time per SIMD slot can be better than of the binary-friendly schemes above. We refer to this performance metric as the amortized running time. -The CKKS scheme [START_REF] Hee Cheon | Homomorphic encryption for arithmetic of approximate numbers[END_REF], which allows to perform computations over approximated numbers, forms the third category. It is similar to the second category in the sense that one can pack several numbers and compute on them in a SIMD manner. The CKKS scheme does not have the algebraic constraints that lower the packing capacity of BGV and BFV. Hence, it is usually possible to pack more elements in CKKS ciphertext, thus resulting with the best amortized cost. Unlike previous schemes, CKKS encodes complex, and thus real, numbers natively. However, homomorphic computations are not exact, which means that decrypted results are only valid up to a certain precision.

Each category of schemes is more efficient for a certain application. Thus, when comparing the efficiency of different homomorphic schemes, one must take into account the given use case.

It is commonly admitted that schemes of the first category are the most efficient ones for generic applications. Since they operate at the bit level, they can compute every logical gate very efficiently. The total running time being in this case the sum of the times needed to evaluate each gate of the circuit. As a result, to optimize the computations for a given application, the only possibility is to reduce the length of the critical computational path and parallelize the related circuit as much as possible. However, as this becomes more and more difficult as the size of the circuit grows, it is possible to optimize only some parts of the circuit by identifying some patterns [START_REF] Aubry | Faster Homomorphic Encryption is not Enough: Improved Heuristic for Multiplicative Depth Minimization of Boolean Circuits[END_REF]. Another advantage of these schemes is that they have very fast so-called 'bootstrapping' algorithms that 'refresh' ciphertexts for further computation. This is very convenient in practice as one can set a standard set of encryption parameters without knowing what function should be computed.

Schemes of the second category operate naturally on p-ary arithmetic circuits, i.e. they are very efficient to evaluate polynomial functions over F p , for a prime p. However, these schemes become much less efficient when considering other kinds of computations, e.g. comparison operations, step functions. To alleviate this problem, one can use tools from number theory to evaluate specific functions with relatively efficient p-ary circuits. Nonetheless, in general this techniques are too weak to outperform schemes of the first category. Moreover, bootstrapping algorithms of these schemes are quite heavy and usually avoided in practice.

CKKS, similarly to second category schemes, is very efficient when operating on arithmetic circuits. However, unlike other schemes which perform modular arithmetic, it allows to perform computations on complex (and thus real) numbers. Although this is an important advantage for many use cases, CKKS lacks simplification tools for evaluation of certain functions due to number-theoretic phenomena as for the second category. However, since CKKS usually supports huge packing capacity, it usually presents the best amortized cost. The bootstrapping algorithm of CKKS is fundamentally different from the above schemes as it refreshes ciphertexts only partially and introduces additional loss of output precision. Therefore, the CKKS bootstrapping is usually avoided in practice.

Although FHE now offers a relatively efficient alternative for secure computation, some functions remain difficult to evaluate efficiently regardless the considered scheme.

Step functions, which are required in many practical applications, form a good example of such functions because of their discontinuous nature. The difficulty to evaluate discontinuous functions comes from the hardness to evaluate a quite basic and relatively simple function: the comparison function. Although comparison is an elementary operation required in many applications including the famous Millionaires problem of Yao [START_REF] Andrew | Protocols for Secure Computations[END_REF] or advance machine learning tasks of the iDASH competition † , it remains difficult to evaluate homomorphically.

By now, schemes of the first category look much more suitable for such nonarithmetic tasks, but they are hopelessly inefficient for evaluating arithmetic functions. Hence, one should resort to heavy conversion algorithms [START_REF] Boura | Chimera: Combining ring-lwe-based fully homomorphic encryption schemes[END_REF] to leverage the properties of different schemes.

Contributions

In this work, we describe the structure of the circuits corresponding to comparison functions for the BGV and BFV schemes. For theses schemes, there exists two approaches: either compare two numbers x and y directly by evaluating a bivariate polynomial in x and y, or study the sign of the difference z = x -y by evaluating a univariate polynomial in z.

By exploiting the structure of these two polynomials, we show that it is possible to evaluate them more efficiently than what was proposed in the state of the art. In particular, we prove that these polynomials have multiple zero coefficients that can be ignored during polynomial evaluation.

The benefit of our approach results in significant performance enhancement for both methods. On the one hand, our bivariate circuit can compare two 64-bit integers with an amortized cost of 21ms, which is a gain of 40% with relation to the best previously reported results of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] (See Table 1). On the other hand, our univariate circuit shows even better results with an amortized cost of 11ms for 64-bit numbers -which is, to the best of our knowledge, more than 3 times faster than previously reported results for this kind of scheme [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]. Note that we can compare two 20-bit numbers with an amortized cost of 3ms, which is better by a factor 1.9 than what can be achieved with CKKS-based algorithms and is comparable to TFHE-based implementations (see Table 5).

We also apply our comparison methods to speed up popular computational tasks such as sorting and computing minimum/maximum of an array with N elements. For example, for N = 64, we obtain an amortized cost of 6.5 seconds to sort 8-bit integers and 19.2 seconds for 32-bit integers, which is faster than the prior work by a factor 9 and 2.5 respectively (see Table 3). For N = 64, we can find the minimum of 8-bit integers with an amortized running time of 404 ms and of 32-bit integers with an amortized time of 9.57 seconds (see Table 4).

Related Art

Comparison is a common function required in many applications; as a consequence, its homomorphic evaluation has been the object of several works. Since inputs are encrypted, a comparison algorithm cannot terminate whenever it finds the first difference between most significant bits. As a result, homomorphic comparison has a complexity corresponding to the worst-case complexity in the plain domain. The practical efficiency of homomorphic comparison depends on the type of HE schemes considered.

For bit-wise HE schemes (FHEW, TFHE), Chillotti et al. [START_REF] Chillotti | Faster packed homomorphic operations and efficient circuit bootstrapping for TFHE[END_REF][START_REF] Chillotti | TFHE: Fast fully homomorphic encryption over the torus[END_REF] showed that one could compare two n-bit integers by evaluating a deterministic weighted automata made of 5n CMux gates. Using the TFHE scheme, evaluating a CMux gate takes around 34 microseconds on a commodity laptop, meaning that one can homomorphically compare two n-bit numbers in around 170n microseconds. Note that these estimations correspond to the fastest (leveled) version of TFHE which avoids bootstrapping. If one wants to use the bootstrapped version then the best method requires to evaluate 7n Mux gates, where each gate takes around 26 millisecond to be evaluated, which makes a total of 182n millisecond.

Schemes from the second category (BGV and BFV) can use SIMD techniques to batch several plaintexts into a single ciphertext [START_REF] Smart | Fully Homomorphic SIMD Operations[END_REF]. Therefore, a natural idea would be to pack the input bits into a single ciphertext. Cheon et al. [START_REF] Hee Cheon | Search-and-compute on encrypted data[END_REF][START_REF] Cheon | Optimized Search-and-Compute Circuits and Their Application to Query Evaluation on Encrypted Data[END_REF] studied comparison functions in this context using the bivariate polynomial interpolation. Some of the algorithmic tools they have used -e.g. computation of running sums and products -are optimal in the homomorphic setting, i.e. regarding the multiplicative depth and the number of multiplications, and have laid the ground for future works in this direction. Some works have tried to exploit other features of these schemes by encoding integers modulo an odd prime p instead of bits. In [START_REF] Narumanchi | Performance analysis of sorting of fhe data: Integer-wise comparison vs bit-wise comparison[END_REF], Narumanchi et al. compare integer-wise comparison algorithms based on the univariate interpolation with bit-wise counterparts. The SIMD packing was ignored in this study. They concluded that bit-wise methods are more efficient because they have a smaller multiplicative depth. In particular, n-bit numbers can be compared with a circuit of depth O(log n) instead of O(n) in the case of integer-wise algorithms. This comes from the fact that integer-wise comparison circuits require to evaluate a Lagrange interpolation polynomial of degree p -1 ≥ 2 n .

In [START_REF] Kim | On the Efficiency of FHE-Based Private Queries[END_REF], Kim et al. noticed that SIMD packing techniques reduce the multiplicative complexity of homomorphic comparison circuits. In addition, they took advantage of the nature of the finite field F p d , which corresponds to the plaintext space of a SIMD slot. Namely, any power x p i can be evaluated with the homomorphic Frobenius automorphism x → x p , which does not consume any homomorphic multiplicative depth level. This allowed to reduce the depth of the equality circuit EQ(x, y) = 1 -(x -y) p d -1 from d log 2 (p) to log 2 (d) + log 2 (p -1) .

Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] proposed a method to perform digit-wise comparison using SIMD and the bivariate polynomial interpolation. Their idea consists in decomposing input integers into digits of size p r encoded into a subfield of F p d , with r|d, in order to reduce the degree of the Lagrange interpolation polynomial of a comparison function. To compare input integers, one should extract digits, compare them and combine the results of digit comparison using the lexicographical order. Note that their evaluation of the lexicographical order intensively uses the efficient equality circuit of Kim et al. [START_REF] Kim | On the Efficiency of FHE-Based Private Queries[END_REF]. Overall, they have used their method to compare integers up to 64-bit while reporting, to the best of our knowledge, the current best timings for performing homomorphic comparison with BGV scheme.

Finally, in [START_REF] Shaul | Secure k-ish nearest neighbors classifier[END_REF], Shaul et al. used the univariate approach to evaluate comparison functions in the context of top-k selection with integer-wise circuits. However, they did not use the decomposition method of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF], thus obtaining relatively poor performance of comparison.

Note that all these works did not exploit the structure of comparison interpolation polynomials neither in the bivariate nor in the univariate case. Kaji et al. [START_REF] Kaji | Polynomial expressions of p-ary auction functions[END_REF] described basic properties of the polynomial expressions of max, argmax and other non-arithmetic functions over non-binary prime fields F p . However, their results do not allow to evaluate these functions very efficiently, as an example their homomorphic circuit to evaluate max has a quadratic complexity in p.

The situation for the CKKS scheme is quite different since its plaintext space natively supports complex/real numbers. Therefore, circuit optimizations related to data encoded into finite fields are not applicable for CKKS. Nonetheless, the approximated nature of computations in CKKS makes it suitable to use iterative methods from real analysis to compute close approximations of non-arithmetic functions. Bajard et al. [START_REF] Bajard | Improving the Efficiency of SVM Classification With FHE[END_REF] used Newton iteration to evaluate the sign function, while independently Cheon et al. [START_REF] Hee Cheon | Numerical method for comparison on homomorphically encrypted numbers[END_REF][START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF] generalized this approach and studied its efficiency in more details. Using the methods of [START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF], one can compare 20-bit numbers with an amortized cost comparable, although slower, to TFHE. However, to obtain these timings one has to use quite large encryption parameters (ring dimension 2 17 and ciphertext modulus up to 2200 bits). Since the running time of homomorphic operations increases quasi-linearly with the ring dimension, while the number of slots only increases linearly, increasing the dimension would affect significantly the timings. Therefore, it would be interesting to know whether these methods can be used in practice to compare larger inputs -e.g. 64-bit integers -without degrading the performance.

Background

Notations

Vectors will be written in column form and denoted by boldface lower-case letters. The set of integers { , . . . , k} is denoted by [ , k]. For a non-negative integer a, let wt(a) be the Hamming weight of its binary expansion. We denote the set of residue classes modulo p by Z p and the class representatives of Z p are taken from the half-open interval [-p/2, p/2).

Comparison of integers with finite fields operations

Let S be a totally ordered set with a binary relation <. For any x, y ∈ S, we can define the less-than and the equality functions as follows.

LT S (x, y) = 1, if x < y; 0, if x ≥ y, EQ S (x, y) = 1, if x = y; 0, if x = y.
Functions over finite fields The map defined by χ : x → x p d -1 from F p d to the binary set {0, 1} is called the principal character. According to Euler's theorem, it returns 1 if x is non-zero and 0 otherwise. Using the principal character, every function from F l p d to F p d can be interpolated by a unique polynomial according to the following well-know lemma.

Lemma 1. Every function f : F l p d → F p d is a polynomial function represented by a unique polynomial P f (X 1 , . . . , X l ) of degree at most p d -1 in each variable. In particular,

P f (X 1 , . . . , X l ) = a∈F l p d f (a) l i=1 (1 -χ(X i -a i )) .
where a i is the ith coordinate of vector a.

Comparison of integers Let

F p d = F p [X]/ f (X) for some irreducible monic polynomial f (X) of degree d. Let d ≤ d and S ⊆ [0, p d -1],
we can map S into F p d by using the decomposition of integers in base p:

ι p : S → F p d , d -1 i=0 a i p i → d -1 i=0 a i X i-1 .
Note that the map ι p is injective and gives a one-to-one correspondence between S and F p d when d = d and S = [0, p d -1]. Therefore, we identify integers belonging to S with their image by ι p and thus omit ι p when the situation is clear from the context. Let a, b ∈ S be two integers to be compared and

d -1 i=0 a i X i and d -1
i=0 b i X i with a i , b i ∈ F p are their respective encodings into F p d . The order of the set S induces a polynomial function in F p d , which can be interpolated by Lemma 1. However, since LT S (0, y) = 1 for any non-zero y ∈ F p d and zero otherwise, we obtain LT S (0, y) = χ(y). Hence, the total degree of the interpolation polynomial is at least p d -1, which might be prohibitive in practice.

Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] proposed an alternative approach where S = [0, p ], they are encoded into F p d as above, but their comparison is performed via the lexicographic order defined on F d p . This method is based on the extraction of coefficients a i , b i ∈ F p thanks to the following result from the theory of finite fields (see [START_REF] Lidl | Introduction to finite fields and their applications[END_REF]Theorem 2.24] for the proof).

Lemma 2. The linear transformations from F p d to F p are exactly the mappings

L α (x) = Tr F p d /Fp (αx) for some α ∈ F p d . Furthermore, L α = L β if α = β.
This lemma implies that for any i ∈ [0, d -1] there exist α i ∈ F p d such that L αi (a) = a i for any a from the vector subspace F d p of F p d . Such α i 's can be computed by solving the following system of equations over

F p d XA = I d where X =      1 1 . . . 1 x x p . . . x p d-1 . . . . . . . . . . . . x d-1 x (d-1)p . . . x (d-1)p d-1      , A =      α 0 α 1 . . . α d-1 α p 0 α p 1 . . . α p d-1 . . . . . . . . . . . . α p d-1 0 α p d-1 1 . . . α p d-1 d-1      and I d ∈ F d×d p d is the identity matrix. Hence, the ith column of X -1 contain the powers α i , α p i , . . . , α p d-1 i
, which define the linear map L αi . Given the input encodings

d -1 i=0 a i X i and d -1 i=0 b i X i ,
we can extract and then compare their vectors of coefficients a = (a 0 , a 1 , . . . , a d -1 ) and b

= (b 0 , b 1 , . . . , b d -1 ) ∈ F d p using the lexicographical order < on F d p defined by a < b ⇔ ∃i ∈ [0, d -1] such that a i < b i and a j = b j ∀j > i .
The corresponding less-than function is equal to

LT S d (a, b) = d -1 i=0 LT S (a i , b i ) d -1 j=i+1 EQ S (a j , b j ) ,
whereas the equality function is defined by

EQ S d (a, b) = d -1 i=0 EQ S (a i , b i ) ,
Notice that the above construction is generic for any set S embedded into F p . For example, if S = [0, s -1] for some s < p, then one can encode input integers via decomposition in base s and compare them using LT S d .

Comparison of large integers

a < b ⇔ ∃i ∈ [0, l -1] such that a i < b i and a j = b j ∀j > i .
As done in [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF], we can employ EQ S d and LT S d to compute the corresponding less-than function LT(a, b) as follows

LT(a, b) = l-1 i=0 LT S d (a i , b i ) l-1 j=i+1 EQ S d (a j , b j ).
(1)

Homomorphic Encryption

We are interested in homomorphic encryption schemes that support SIMD operations on their plaintexts. This section aims at giving the necessary background regarding these schemes.

Cyclotomic fields and Chinese Remainder Theorem Let m be a positive integer and n = ϕ(m) where ϕ is Euler's totient function. Let K = Q(ζ m ) be the cyclotomic number field constructed by adjoining a primitive m-th root of unity ζ m ∈ C to the field of rational numbers. The ring of integers of K, denoted by R, is isomorphic to Z[X]/ Φ m (X) where Φ m (X) is the m-th cyclotomic polynomial. Let p > 1 be a prime number coprime to m, then Φ m (X) splits modulo p into irreducible factors of same degree

d, i.e. Φ m (X) = F 1 (X) • • • F (X) mod p.
The degree d is actually the order of p modulo m, and = n/d. The Chinese Remainder Theorem (CRT) states that in this case the following ring isomorphism holds:

R p = Z p [X]/ Φ m (X) ∼ = Z p [X]/ F 1 (X) × . . . × Z p [X]/ F (X) For each i ∈ [1, ] the quotient ring Z p [X]/ F i (X) is isomorphic to the finite field F p d .
Hence, the above isomorphism can be rewritten as R p ∼ = F p d . We call every copy of F p d in this direct product a slot. Therefore, every element of R p contains slots, which implies that an array of independent F p d -elements can be encoded as a unique element of R p . The slot isomorphic to

Z p [X]/ F i (X)
is referred to as the ith slot.

Additions and multiplications of R p -elements results in the corresponding coefficient-wise operations of their respective slots. In other words, each ring operation on R p is applied to every slot in parallel, which resembles the Single-Instruction Multiple-Data (SIMD) instructions used in parallel computing. Therefore, the above encoding method from F p d to R p is often called the SIMD packing.

The HE schemes that support SIMD packing and exact computations over encrypted data include BGV [START_REF] Brakerski | Leveled) Fully Homomorphic Encryption without Bootstrapping[END_REF] and FV [START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF]. These schemes have a common framework described below.

Basic setup Let λ be the security level of an HE scheme. Let L be the maximal multiplicative depth of homomorphic circuits we want to evaluate. Let d be the order of the plaintext modulus p modulo the order m of R. Assume that the plaintext space R p has SIMD slots, i.e. R p ∼ = F p d . The basic part of any HE schemes consists of key generation, encryption and decryption algorithms.

KeyGen(1 λ , 1 L ) → (sk, pk). Given λ and L, this function outputs the secret key sk and the public key pk.

Encrypt(pt ∈ R p , pk) → ct. The encryption algorithm takes a plaintext pt and the public key pk and outputs a ciphertext ct.

Decrypt(ct, sk) → pt. The decryption algorithm takes a ciphertext ct and the secret key sk and returns a plaintext pt. For freshly encrypted ciphertexts, the decryption correctness means that Decrypt(Encrypt(pt, pk), sk) = pt.

Homomorphic operations

The homomorphic addition (multiplication) algorithm takes two input ciphertexts ct 1 and ct 2 encrypting plaintexts pt 1 and pt 2 respectively. It outputs a ciphertext ct that encrypts the sum (product) of these plaintexts in the ring R p . It implies that homomorphic addition (multiplication) sums (multiplies) respective SIMD slots of pt 1 and pt 2 . Similar operations between ciphertexts and plaintexts are defined as well.

Every homomorphic ciphertext contains a special component called noise that is removed during decryption. However, the decryption function can deal only with noise of small enough magnitude; otherwise, this function fails. This noise bound is defined by encryption parameters in a way that larger parameters result in a larger bound. The ciphertext noise increases after every homomorphic operation and, therefore, approaches its maximal possible bound. It implies that to reduce encryption parameters one needs to avoid homomorphic operations that significantly increase the noise. Therefore, while designing homomorphic circuits, we need to take into account not only the running time of homomorphic operations but also their effect on the noise.

The most expensive homomorphic operation with relation to both noise and running time is ciphertext-ciphertext multiplication (Mul). This operations takes place when two expressions containing input values are multiplied. Such multiplication is called non-scalar. In contrast, ciphertext-plaintext multiplication (MulPlain) is used when an expression with input values is multiplied by an unencrypted or publicly known value. This is a scalar multiplication. Since Mul is much more expensive than MulPlain, the multiplicative depth and complexity of a homomorphic circuit is usually calculated with relation to the number of Mul's, or non-scalar multiplications. Thus, in the following sections we focus on the non-scalar complexity of comparison circuits.

Optimising the comparison circuits over F p

In this section, we study the structure of basic comparison circuits over F p used in Section 2 to compare integers, namely LT S and EQ S for some S ⊆ [0, p -1].

For any choice of S, the corresponding equality function over F p is equal to

EQ S (x, y) = 1 -(x -y) p-1 .
Unfortunately, LT S is not that simple and universal and we have to rely on Lagrange interpolation (Lemma 1) to compute it. Yet, there are two different ways to evaluate it. The first method (as done in [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]) uses S = [0, p -1] and directly interpolates LT Fp (x, y) as a bivariate polynomial over F p . The second approach (as done in [START_REF] Narumanchi | Performance analysis of sorting of fhe data: Integer-wise comparison vs bit-wise comparison[END_REF] and [START_REF] Shaul | Secure k-ish nearest neighbors classifier[END_REF]) has S = [0, (p -1)/2] and interpolates a univariate polynomial of LT S (z, 0) with z = x -y. In this section, we show how to exploit the structure of these polynomials to speed-up their evaluation.

Bivariate interpolation of LT S .

Let S = [0, p -1]. The less-than function can be interpolated using Lemma 1 and the following truth table.

< 0 1 2 • • • p -1 0 0 1 1 • • • 1 1 0 0 1 • • • 1 2 0 0 0 • • • 1 . . . . . . . . . . . . . . . . . . p -1 0 0 0 • • • 0
In particular, the interpolation polynomial of LT S over F p is equal to

P LT S (X, Y ) = p-2 a=0
EQ S (X, a)

p-1 b=a+1 EQ S (Y, b) = p-2 a=0 1 -(X -a) p-1 p-1 b=a+1 1 -(Y -b) p-1 .
Surprisingly, the total degree of P LT S (X, Y ) is only p and its coefficients can be described by the following theorem.

Theorem 1. Let p > 2 be a prime number and S = [0, p -1], then the interpolation polynomial of LT S over F p has the following form

P LT S (X, Y ) = Y p-1 - p -1 2 (XY ) p-1 2 + i,j>0, i =j, i+j≤p a ij X i Y j
where a ij = p-2 a=0

p-1 b=a+1 a p-1-i b p-1-j ∈ F p . The total degree of P LT S (X, Y ) is p.

Proof. See Appendix A.

From the definition of LT S , one can easily prove the following facts about P LT S :

-

P LT S (X, 0) = 0, thus Y divides P LT S (X, Y ); -P LT S (X, X) = P LT S (Y, Y ) = 0, thus (X -Y ) divides P LT S (X, Y ); -P LT S (p -1, Y ) = 0 thus X + 1 divides P LT S (X, Y ).
Hence, there exists a bivariate polynomial f (X, Y ) of total degree p -3 over F p such that:

P LT S (X, Y ) = Y (X -Y )(X + 1)f (X, Y ). ( 2
)
The following theorem describes the structure of f (X, Y ).

Theorem 2. Let p be an odd prime and S = [0, p -1]. Let P LT S (X, Y ) be the interpolation polynomial of LT S over F p and

P LT S (X, Y ) = Y (X -Y )(X + 1)f (X, Y ).
Then, for any z ∈ F p we have

f (z, z) = f (z, 0) = f (p -1, z). (3)
As a consequence, there exists (p -1)/2 polynomials f i (X) over F p , 0 ≤ i ≤ (p -3)/2, such that:

f (X, Y ) = (p-3)/2 i=0 f i (X)Z i , (4)
with

Z = Y (X -Y ) and deg(f i (X)) = p -3 -2i.
Since our proof of Theorem 2 is quite long and with no real interest for the purpose of this work, we defer it to an extended version of this paper. In our experiments, we used the decompositions (4) of f (X, Y ) only for small p (between 3 and 7), which you can find in Appendix B.

Complexity analysis. In [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF], the authors proposed to evaluate P LT S (X, Y ) by evaluating each monomials separately before summing them up. Given x, y ∈ F p , they precompute the powers of x and y up to p -1 for a total of 2p -4 non-scalar multiplications. Then, another p -1 non-scalar multiplications are needed to evaluate each monomial (( i c i X i )Y j ) j where c i 's are scalars in F p , before summing them together to get the final result. Overall their evaluation of P LT S (X, Y ) requires 3p -5 non-scalar multiplications.

Following this idea and using the decomposition of f (X, Y ) given in Eq. ( 4) one needs:

-2 multiplications to compute (X + 1)Z if p ≥ 3 or 1 multiplication when p = 2;

and then for p ≥ 5:

p -4 multiplications to compute the X i 's for 2 ≤ i ≤ p -3 required to compute the terms f i (X); -(p -5)/2 multiplications to compute the Z i for 2 ≤ i ≤ (p -3)/2; -(p -5)/2 multiplications to compute the products f i (X)

• Z i ; -1 final multiplication (X + 1)Z • f (X, Y ).
Overall, at most 2p -6 non-scalar multiplications are needed to homomorphically evaluate P LT S (X, Y ) for p ≥ 5. This number can be slightly reduced by optimizing the way of computing f i 's. For instance, it can be done with only 6 multiplications for p = 7, (see Appendix B), which is smaller than 2p -6 = 8.

Overall, for p ≥ 5, our method saves p + 1 multiplications over the method of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]. However, the complexity of the bivariate circuit remains linear in p which is unpractical for performing homomorphic comparisons using large digits (i.e. a large p).

Univariate interpolation of LT S

Unlike bivariate polynomials, it is possible to evaluate univariate polynomials of degree p -1 in O( √ p) non-scalar multiplications using the Paterson-Stockmeyer algorithm [START_REF] Michael | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF]. The Paterson-Stockmeyer algorithm has been used in various works related to homomorphic encryption in order to speed-up polynomial evaluation. As a recent example, it was applied by Shaul et al [START_REF] Shaul | Secure k-ish nearest neighbors classifier[END_REF] in the context of top-k selection, which uses small-number comparison as a subroutine. However, in our case the study of the structure of LT S as a univariate polynomial will allow us to speed-up its evaluation for large p beyond what could be achieved using only the Paterson-Stockmeyer algorithm.

To evaluate LT S as a univariate polynomial, we compute the difference x -y of the two input values and check its sign. To compute the sign function using finite field arithmetic, we need to split finite field elements into two classes: negative (F - p ) and non-negative (F + p ). In addition, for any x, y ∈ S the following property should hold:

x -y ∈ F + p if LT S (x, y) = 0, F - p if LT S (x, y) = 1.
It is easy to see that these constraints are satisfied by

S = [0, (p -1)/2]. Let us split F p into F + p = [0, (p -1)/2] and F - p = [-(p -1)/2, -1].
Notice that for any x, y ∈ S, their difference x -y belongs to F - p if and only if x < y. Let χ F - p (z) be a function that outputs 1 if z is negative and 0 otherwise.

According to Lemma 1, χ F - p (z) is equal to χ F - p (z) = -1 a=-p-1 2 1 -(z -a) p-1 .
Combining the above facts, the LT S function can be interpolated by the following polynomial over

F p Q LT S (X, Y ) = -1 a=-p-1 2 1 -(X -Y -a) p-1 .
The following theorem describes all the coefficients of this interpolation polynomial.

Theorem 3. For an odd prime p and S = [0, (p -1)/2], the LT S function can be interpolated by the following polynomial over

F p Q LT S (X, Y ) = p + 1 2 (X -Y ) p-1 + p-2 i=1,odd c i (X -Y ) i . ( 5
)
where

c i = p-1 2
a=1 a p-1-i .

Proof. See Appendix C Remark 1. The polynomial Q LT S (X, Y ) yields the interpolation polynomial of the sign function sgn S defined on S = [-(p -1)/2, (p -1)/2] as sgn S (x) = 1 if x < 0 and sgn S (x) = 1 if x ≥ 0. In particular, we have

Q sgn S (X) = Q LT S (X, 0) . (6)
Complexity analysis. The above theorem implies that the less-than function can be expressed by a univariate polynomial of degree p -1. In general, such polynomials are evaluated in p-1 multiplications according to Horner's method.

To reduce the number of non-scalar multiplications, we can resort to the Paterson-Stockmeyer algorithm [START_REF] Michael | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF] that requires 2(p -1) + log 2 (p -1) + O(1) such multiplications. However, we can improve this complexity by exploiting the fact that the polynomial in (5) has only one coefficient with an even index, the leading one. Thus, if Z = X -Y , we can rewrite (5) as follows

α p-1 Z p-1 + Z p-3 i=0,even α i+1 Z i = α p-1 Z p-1 + Zg(Z 2 )
where

α i = p-1 2
a=1 a p-1-i and g(X) is a polynomial of degree (p -3)/2. To evaluate g(X), the Paterson-Stockmeyer algorithm requires

√ p -3+log 2 p-3 2
+O(1) non-scalar multiplications. Furthermore, the preprocessing phase of this algorithm computes the powers Z 2 , Z 4 , . . . , Z 2k and Z 4k , Z 8k , . . . , Z 2 r k with 2k(2 r -1) = p -3. We can use these powers to compute the leading term in r non-scalar multiplications, namely

Z 2 Z 2k Z 4k • • • Z 2 r k = Z 2+2k(2 r -1) = Z 2+p-3 = Z p-1 .
Since the optimal k is about (p -3)/2, we obtain that r must be about log 2 √ p -3. Hence, the total non-scalar complexity of evaluating ( 5) is equal to

p -3 + 3 log 2 (p -3) 2 + O(1).
Remark 2. A careful reader can notice that the leading term of ( 5) is equal to (X -Y ) p-1 , which is the heaviest part of the equality circuit EQ S (X, Y ). Thus, we can get EQ S (X, Y ) almost for free (at the cost of one homomorphic subtraction) after evaluating LT S (X, Y ), which saves O(log(p-1)) non-scalar multiplications. This feature of the univariate circuit allows to compute all the equality operations while comparing large integers using the less-than function LT from (1). This saves O((d -1)(k-1) log(p-1)) homomorphic multiplications, thus leading to a better running time than for the bivariate circuit.

The downside of the univariate circuit is that only (1/2) d of the plaintext space is used to encode input integers in comparison to the bivariate method.

Min/max function

Given the less-than function LT defined on some set, one can compute the minimum of two elements x, y of this set in the following generic way min(x, y) = x • LT(x, y) + y • (1 -LT(x, y))

= y + (x -y) • LT(x, y) . (7)
Notice that the input difference x -y naturally emerges in this expression, thus hinting that the univariate circuit from ( 5) might be useful here. Indeed, by replacing X -Y with a variable Z we obtain the univariate polynomial representation of the minimum function on the set S = [0, (p -1)/2]

Q min S (X, Y ) = Y + Z • Q LT S (X, Y ) = Y + p + 1 2 Z + p-1 2 i=1 Z 2i p-1 2 a=1 a p-2i = p + 1 2 (X + Y ) + p-1 2 i=1 Z 2i p-1 2 a=1 a p-2i = p + 1 2 (X + Y ) + g(Z 2 ) ,
where g(X) is a polynomial of degree (p -1)/2. As a result, min S (x, y) can be computed with O( √ p -1) non-scalar multiplications via the Paterson-Stockmeyer algorithm.

Following the above reasoning, the maximum function can be computed with the following polynomial

Q max S (X, Y ) = p + 1 2 (X + Y ) - p-1 2 i=1 Z 2i p-1 2 a=1 a p-2i .
Remark 3. Maximum and minimum functions are basic building blocks in the design of neural networks. For example, one of the most popular activation functions in neural networks is the rectifier, or ReLU, which is equal to max(x, 0). By analogy with (7), we have max(x, 0) = x • (1 -sgn S (x)) where S = [-(p -1)/2, (p -1)/2] (see Remark 1). Thus, Eq. ( 6) yields the following interpolation polynomial of the ReLU function on S

Q ReLU S (X) = X • (1 -Q sgn S (X)) = p + 1 2 X - p-1 2 i=1 X 2i p-1 2 a=1
a p-2i .

Impact on the overall complexity

In this section, we summarize the complexities of our method as compared to the work of Tan et al. for the evaluation of the less-than function using the lexicographical order method described in Section 2.

For fixed p and d, Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]Section 4.4] determined that the depth of the circuit evaluating the less-than function of two b-bits integers is equal to

log 2 d + log 2 (p -1) + log 2 (log p 2 b )/d) + 4. ( 8
)
Note that our algorithms do not decrease the depth of the circuit. Similarly, Tan et al. showed that the number of homomorphic multiplications required to evaluate the less-than function of two b-bit integers with the bivariate interpolation is

d • (T + log 2 (p -1) + log 2 d ) + log 2 (log p 2 b )/d) + 2.
where T is the number of homomorphic multiplications required to evaluate the comparison circuit over F p . In the work of Tan et al. T = 3p -5, while in our case it is T = 2p -6.

The univariate method saves even more multiplications since one can extract 1 -EQ S (x, y) while computing LT S (x, y). Hence, we obtain EQ S (x, y) almost for free when evaluating the lexicographical order (Remark 2). Thus, in this case the comparison of two b-bits integers requires

d • (T + log 2 d ) + log 2 (log p 2 b )/d) + 2.
with T ≈ √ 2p -4 + 3(log 2 (2p -4))/2 by using the Paterson-Stockmeyer algorithm.

Applications

The previous results can help improving the performance of any task involving comparisons performed homomorphically such as private database queries, knearest neighbour search, top-k selection or step function evaluation in neural nets. In this section we choose to demonstrate the gain brought by our approach for sorting and min/max search which are subroutines needed for the aforementioned tasks.

Sorting

To demonstrate the efficiency of our comparison algorithms, we applied them to a popular computational task that demands multiple comparisons, sorting. The best homomorphic sorting algorithm in terms of running time is the direct sorting algorithm due to C ¸etin et al. [START_REF] Gizem | Depth optimized efficient homomorphic sorting[END_REF]. For a given array A = [a 0 , . . . , a N -1 ], this algorithm computes a comparison matrix L defined by

L ij =      LT(a i , a j ) if i < j, 0 if i = j, 1 -LT(a j , a i ) if i > j.
Example : for A = [5, 1, 7, 2, 3], the matrix L ∈ {0, 1} 5×5 is given by:

L =       0 0 1 0 0 1 0 1 1 1 0 0 0 0 0 1 0 1 0 1 1 0 1 0 0      
It is easy to see that the Hamming weight of the ith row of L is unique and equal to the array index of a i after sorting the array A in the descending order. For example, the zero weight indicates that there are no elements of A bigger than a i . Thus, a i has a zero index in A after sorting; in other words, a i is the maximum element of A.

Let A be a sorted version of A in the descending order. To compute A [i] for any i ∈ [0, N -1], we homomorphically select an element a j such that wt(L[j]) = i. This can be done with the following sum

A [i] = N -1 j=0 EQ [0,N -1] (i, wt (L[j])) • a j . ( 9
)
Note that the equality function should be defined on the set [0, N -1], which implies that N must be smaller than the plaintext modulus p.

Remark 4. Since the matrix L is defined by N (N -1)/2 elements, it can be costly to keep it in memory for large N . Instead, we can compute the Hamming weights of its rows by iteratively computing one comparison LT(a i , a j ) with i < j at a time. To achieve this, we create an array of size N initialized with zeros that eventually will store the Hamming weights. Then, we add the outcome of LT(a i , a j ) to the ith element of this array and the result of 1 -LT(a i , a j ) to the jth element. In this approach, only N elements of the Hamming weight array are being kept in RAM.

The direct sorting algorithm requires N (N -1)/2 less-than operations to compute the matrix L and N 2 equality operations to compute sorted elements of A . While computing equalities, we can reduce the total number of non-scalar multiplications if N is large enough. Recall EQ S needs M = log 2 (p -1) + wt(p -1) -1 non-scalar multiplications for any S ⊆ [0, p -1]. Hence, to compute EQ [0,N -1] (i, wt (L[j])) for all i ∈ [0, N -1], we should perform N M multiplications. Using Lemma 3 (see Appendix A), we can rewrite:

EQ [0,N -1] (i, wt (L[j])) = 1 - p-1 k=0 i k • wt (L[j]) p-1-k .
If we precompute the powers wt (L[j]) p-1-k , then we need only p -2 non-scalar multiplications to compute all the equalities EQ(i, wt (L[j])) as the index i is not encrypted. Hence, if N > (p -2)/M , this approach results in a smaller number of non-scalar multiplications. Yet, this method introduces p -1 scalar multiplications (by powers i k ) and p -2 additions. However, these operations are much faster in HE schemes than non-scalar multiplication such that the gain from reducing non-scalar multiplications becomes dominant.

The main advantage of direct sorting is that its multiplicative depth is independent of the array length, namely d = d (LT) + log 2 (p -1) + 1 with d (LT) given in Equation (??).

This allows to avoid large encryption parameters and costly bootstrapping operations.

Minimum and maximum of an array

Another application of our comparison algorithms is concerned with finding a minimum (or maximum) element of an array. To find the minimum of an array with N elements, at least N -1 calls of the pairwise minimum function are required [START_REF] Thomas H Cormen | Introduction to algorithms[END_REF]Chapter 9], which can be achieved, for instance, by the tournament method.

The tournament method consists of log N iterations. In each iteration, the input array is divided into pairs. If the array length is odd, one element is stashed for the next iteration. Then, the maximum of each pair is removed from the array. The algorithm stops when only one element is left; this is the minimum of the input array, see Figure 1. Unfortunately, the tournament method has a Fig. 1: The tournament method of finding the minimum of an array. In each stage, the array elements are divided into pairs. Only minimum of a pair go to the next stage.

big multiplicative complexity, namely log N • d(min(x, y)). In the HE world, this enforces us to use either impractical encryption parameters [START_REF] Togan | Comparison-based applications for fully homomorphic encrypted data[END_REF] or a slow bootstrapping function.

To reduce the depth of the array minimum algorithm, we can combine the tournament method and direct sorting. Let A = [a 0 , . . . , a N -1 ] be an input array. First, we perform T iterations of the tournament algorithm, which leaves us with an array A = [a 0 , . . . , a N -1 ] of length N = N/2 T containing minimal elements of A. Then, we can find the minimum by computing the comparison table L as in direct sorting and extracting one of the minimal elements. If M (f ) is the non-scalar multiplicative complexity of a function f , then the total number of non-scalar multiplications to find the minimum of an array is approximately equal to

(N -N ) • M (min(x, y)) + N (N -1) 2 • M (LT) +M (Extraction).
The extraction of a minimum element can be done with two methods. In the first approach, we use the fact that the Hamming weight of the comparison table row corresponding to the minimum is equal to N -1. Hence, we can retrieve the minimum as in ( 8)

min(A ) = N -1 i=0 EQ [0,N -1] (N -1, wt (L[i])) • a i . (10)
Here, the multiplicative depth is equal to

T • d(min(x, y)) + d(LT) + log 2 (p -1) + 1
which is independent of the input array length N . Since EQ [0,N -1] need log 2 (p -1) + wt(p -1) -1 non-scalar multiplications, then the number of non-scalar multiplications needed to extract the array minimum is equal to

M (Extraction) = N (log 2 (p -1) + wt(p -1) -1).
Shaul et al. [START_REF] Shaul | Secure k-ish nearest neighbors classifier[END_REF] proposed another circuit to extract the array minimum that exploits the fact that the comparison table row related to the minimum contains only 1 except for the main diagonal entry. In other words, the product of N j=1,j =i L ij = 1 if and only if a i = min(A ). Hence, the minimal element is equal to

min(A ) = N -1 i=0 a i • N j=1,j =i L ij . (11)
The resulting depth of this circuit amounts to

T • d (min(x, y)) + d(LT) + log(N -1) + 1.
This extraction circuit requires the following number of multiplications

M (Extraction) = N (N -2).
This implies that for small enough N , Shaul's circuit [START_REF] Brakerski | Fully homomorphic encryption without modulus switching from classical GapSVP[END_REF] has a smaller depth or/and a smaller multiplication complexity than the circuit in [START_REF] Bourse | Fast Homomorphic Evaluation of Deep Discretized Neural Networks[END_REF]. Furthermore, Shaul's circuit supports any length N > p, whereas (9) requires N ≤ p such that wt(L[i]) do not overflow modulo p.

In the experiments conducted in Section 5, we use the best of these approaches for given N , T and p.

Implementation Results

We implemented the lexicographic order algorithm (LT, Eq. 1) using the BGV scheme [START_REF] Brakerski | Leveled) Fully Homomorphic Encryption without Bootstrapping[END_REF]. The code is written in the HElib library [START_REF] Helib | An implementation of homomorphic encryption[END_REF]. For a fair comparison with the prior work, we also implemented the algorithm of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]. The code is publicly available via https://github.com/iliailia/ comparison-circuit-over-fq. In all the experiments, we used an average commodity laptop equipped with an Intel Dual-Core i5-7267U CPU (running at 3.1 GHz) and 8 GB of RAM. Multi-threading was turned off.

In the results presented below, the following notation is used:

p : the plaintext modulus of the BGV scheme; q : the initial ciphertext modulus of the BGV scheme; m : the cyclotomic order of the ring R; n : the degree of the ring R (n = φ(m)); : the number of SIMD slots; d : the dimension of a slot subspace used for digit encoding (d in Eq. 1); l : the dimension of digit vectors encoding input integers over F d p ; k : the number of input integers encoded in one ciphertext (k = /l ).

In all the experiments, the encryption parameters of BGV are chosen according to the following strategy. The plaintext modulus p is a prime number such that SIMD slots are isomorphic to a finite field. Next, we choose the order m of R with large enough n to support our homomorphic algorithms (n > 12000). The order of p modulo m should be as small as possible, which maximizes the number of SIMD slots, thus reducing the amortized running time of our algorithms. In addition, m is chosen to be a prime number or a product of a few primes. This constraint makes sure that the slot permutation group is cyclic or a product of a few cyclic groups, which results in a better performance (for more details, see [START_REF] Gentry | Smart. Fully Homomorphic Encryption with Polylog Overhead[END_REF]

, Appendix C.3]).
Every input integer is encoded into exactly one ciphertext. First, it is decomposed into a vector of l digits over F d p (see Section 2.2). Then each digit embedded into one SIMD slot. Thus, each ciphertext encrypts exactly k = /l integers.

Note that although our algorithms save some homomorphic multiplications as compared to [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF], the depth of the circuits remains unchanged. Since the ciphertext size depends mainly on the depth of the circuit one wants to evaluate and the desired level of security, we obtain similar bandwidth requirements than previous works.

To compare our algorithms with the state of the art, we ran the lexicographic order algorithm (LT) to compute the less-than function on 64-bit integers. The results of these experiments are presented in Table 1. In addition, we ran the direct sorting and the array minimum algorithms from Section 4; see Tables 3 and4, respectively.

We ran our algorithms with p ∈ [START_REF] Martin R Albrecht | On the concrete hardness of learning with errors[END_REF]659]. However, Table 1 contains the best results for small p's, where the bivariate and the univariate comparison circuits have a comparable running time. Table 2 contains the results with the best observed amortized running time. For the sorting and the array minimum applications (Tables 3 and4), we showed only the results with p and m supporting reasonable security levels and giving the best running time for arrays of length N = 64. Fig. 2: Less-than function via different methods. The running time speedup factor of our lexicographic order algorithms over the algorithm of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]. These factors are computed using the data from Table 1.

As shown in Table 1 and Figure 2, our lexicographic order circuits have a better running time than the prior work of Tan et al. [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] even for small plaintext moduli. In particular, using our bivariate circuit, the less-than function is 1.71 times faster than with the circuit of Tan et al. on their fastest set of parameters ((p, d, l) = (5, 7, 4)). The best running time per integer was achieved by our univariate circuit, which outperforms any bivariate circuit for any p > 5 at the cost of larger encryption parameters. As shown in Figure 2, the speedup [START_REF] Martin R Albrecht | On the concrete hardness of learning with errors[END_REF]. Note that the amortized timing per slot from [START_REF] Gizem | Depth optimized efficient homomorphic sorting[END_REF] is obtained with the LTV scheme [START_REF] López-Alt | On-the-fly multiparty computation on the cloud via multikey fully homomorphic encryption[END_REF], which was attacked by Albrecht et al. [START_REF] Albrecht | A subfield lattice attack on overstretched NTRU assumptions -cryptanalysis of some FHE and graded encoding schemes[END_REF].

(p, m, n) Type (d, l) log 2 q λ k
Table 3 illustrates that our homomorphic sorting implementation achieves the best running time in the existing literature. Note that the best result [START_REF] Gizem | Depth optimized efficient homomorphic sorting[END_REF] in this area is based on an SHE scheme that was successfully attacked by Albrecht et al. [START_REF] Albrecht | A subfield lattice attack on overstretched NTRU assumptions -cryptanalysis of some FHE and graded encoding schemes[END_REF]. Hence, this result is hard to compare directly to our work. The existing literature on homomorphic array minimum/maximum algorithms [START_REF] Togan | Comparison-based applications for fully homomorphic encrypted data[END_REF][START_REF] Shaul | Secure k-ish nearest neighbors classifier[END_REF] is based on the techniques described in Section 4.2. Hence, our improvement of comparison circuits automatically results in a better performance over these works. For example, Togan et al. [START_REF] Togan | Comparison-based applications for fully homomorphic encrypted data[END_REF] needed 346.9 seconds to find maximal elements of 960 arrays of 16 8-bit integers, which is 361 milliseconds per array. Our work can perform this task in 74 milliseconds, see Table 4. To find the minimum of 64 32-bit integers, our array minimum algorithm requires about 9.5 seconds.

Comparison to other HE schemes

As mentioned in the introduction, there are three types of HE schemes suitable in different use cases including TFHE, CKKS and BGV/BFV. Our algorithms are designed for BGV/BFV, which support SIMD packing and are the most efficient FHE schemes for exact computation in arithmetic circuits. However, the amortized running time per data value of our comparison algorithms is comparable to efficient FHE schemes for binary circuits (TFHE) and HE supporting approximate arithmetic over complex numbers (CKKS).

In Table 5, our implementation of the less-than function is compared to the implementations of this function in TFHE [START_REF] Chillotti | Faster packed homomorphic operations and efficient circuit bootstrapping for TFHE[END_REF][START_REF] Chillotti | TFHE: Fast fully homomorphic encryption over the torus[END_REF] and CKKS [START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF].

Chilloti et al. [START_REF] Chillotti | Faster packed homomorphic operations and efficient circuit bootstrapping for TFHE[END_REF][START_REF] Chillotti | TFHE: Fast fully homomorphic encryption over the torus[END_REF] constructed a deterministic weighted automata that can compute the maximum function using the TFHE scheme. The same automata can compute the less-than function without any performance loss. In this case, the running time of the less-than function of two b-bit integers takes 170b microseconds on a hardware similar to ours and with the encryption parameters supporting at least 152 bits of security. This security level might be lowered by a recent attack of Espitau et al. [START_REF] Espitau | On a hybrid approach to solve small secret LWE[END_REF]. Note that this running time is achieved in the leveled mode of TFHE, i.e. without bootstrapping. Unfortunately, we could not manage to run the code of this implementation on our machine.

Cheon et al. [START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF] designed a polynomial approximation of the less-than function over real numbers. Since the precision of this approximation depends on the ciphertext noise which cannot be reduced in CKKS, only a few consecutive comparisons are possible to perform correctly unlike in TFHE and BGV/BFV. Nevertheless, the number of SIMD slots in CKKS is always half of the ring dimension (n/2), which significantly reduces the amortized running time per value.

We ran the comparison method of Cheon et al. on our machine with multithreading turned off. Since the implementation in [START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF] uses parallelization with 8 cores, our running time of this method presented in Table 5 is significantly larger than in [START_REF] Hee Cheon | Efficient homomorphic comparison methods with optimal complexity[END_REF]. The encryption parameters of the CKKS scheme are set to support a security level of at least 128 bits.

Our implementation runs with p = 131, m = 17293 (n = 17292), which corresponds to 5764 integers of 8-16 bits or 2882 20-bit integers encoded into one ciphertext. The encryption parameters corresponds to a security level of at least 126 bits.

As shown in Table 5, our algorithm for the less-than function demonstrates a similar performance as the TFHE-based implementation and is up to 2 times faster than the CKKS-based work. This means that in use cases that involve arithmetic and non-arithmetic functions (e.g. artificial neural networks) one might resort to only an HE scheme supporting exact arithmetic circuits (i.e. BGV/BFV) instead of combining it with an HE scheme efficient for a nonarithmetic part of the computation (i.e. TFHE).

Conclusion

In this work, we constucted more efficient homomorphic circuits of comparison operations for the BGV and BFV FHE schemes. Our results are based on structural properties of comparison functions over finite fields. We proved that less-than functions of two input variables x and y can be represented either by bivariate polynomials or univariate polynomials (in variable z = x-y) with multiple zero coefficients, which simplifies computation. Moreover, our computation of the univariate less-than functions yields the output of the equality function almost for free, which allowed us to speed up the lexicographic order of vectors over finite fields and thus comparison of large integers encoded by these vectors.

The implementation of our circuits in HElib is faster than the state-of-theart work [START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] by more than a factor of 3. Furthermore, the running time of our circuits is comparable to implementations of comparison algorithms in TFHE, which implies that the total and amortized running time of this scheme are the same. The encryption parameters of each scheme are set to support the following security levels: 152 bits for TFHE (might be smaller due to [START_REF] Espitau | On a hybrid approach to solve small secret LWE[END_REF]), 128 bits for CKKS and 126 bits for BGV. *TFHE timings are estimated from [START_REF] Chillotti | TFHE: Fast fully homomorphic encryption over the torus[END_REF].

which is believed to be the most efficient FHE scheme for non-arithmetic homomorphic computations. As a side contribution, we applied our comparison algorithms to the tasks of sorting and array minimum search. In both cases we achieved the best running time present in the literature. For instance, our sorting algorithm can sort 4676 batches of 64 32-bit integers in about 25 hours, which results in an amortized time equal to about 19 seconds per batch. Our minimum circuit can find minimal elements of 2610 batches of 64 32-bit integers in approximately 7 hours; the amortized time is 9.5 seconds per batch. We hope that this work will draw attention to the study of arithmetic circuits over finite fields representing non-arithmetic functions over integers, thus leading to practically efficient homomorphic implementations of useful algorithms.

A Proof of Theorem 1

To prove Theorem 1 we need the following lemmas. = (-1) p-1-i mod p , we prove the lemma.

Lemma 4. Let P (X) be a polynomial of degree d less than p -1. For any prime number p > 2, it holds p-1 a=0 P (a) = 0 mod p.

Proof. Since

p-1 a=0 b = 0 mod p for any b ∈ F p , it is enough to prove that the sum p-1 a=0 a n = 0 mod p for any 0 ≤ n < p -1. Since the case n = 0 is straightforward, let us assume n > 0. Let g be a primitive element of F p . Since p > 2, we have g = 1. Thus, we can rewrite the above sum as follows.

p-1 i=1 g in = g pn -g n g n -1 .

Since g pn ≡ g n mod p, the sum turns into zero modulo p. Now, we have all the ingredients to prove Theorem 1.

Proof (Proof of Theorem 1). Assume that all computations are done modulo p.

Using Lemma 3, we obtain that P LT S (X, Y ) is equal to -a p-1-i + a p-1-i b p-1 = 0 .

Next, we compute coefficients by Y i with i > 0. If i = p -1, this sum is equal to 1. According to Lemma 4, it is 0 if i < p -1.

To compute coefficients by X i Y j with i, j > 0, we will use Faulhaber's formula below

Note that P (0) = 0. The coefficient by X i Y j for some positive i and j is equal to According to [START_REF] Brakerski | Leveled) Fully Homomorphic Encryption without Bootstrapping[END_REF], there exist a polynomial P i (X) of degree p -i such that b-1 a=0 a p-1-i = P i (b). Since Q ij (X) = X p-1-j P i (X) has degree 2p -1 -i -j, Lemma 4 implies that if i + j > p, then

Q ij (b) = p-1 b=0 Q ij (b) -Q ij (0) = 0.
Thus, all the coefficient X i Y j with i + j > p are zero, which means that the total degree of P LT S (X, Y ) is at most p.

In addition, we consider the case when i = j and i, j ≤ (p -1)/2. Let us consider the following sum Note that the inverse of 2 is well defined modulo p since p is an odd prime. If i < (p -1)/2, then Lemma 4 says the sum on the right side is zero. Thus, the coefficient by X i Y i , which is exactly the sum on the left side, is equal to zero. If i = (p -1)/2, the above equality yields that the coefficient by (XY ) (p-1)/2 is equal to -(p -1)/2.

B Decomposition of f (X, Y ) for 3 ≤ p ≤ 7

Let Z = Y (X -Y ). One non-scalar multiplication is needed to compute Z. p=3. f (X, Y ) = 2 .

Since the polynomial f (X, Y ) is constant, it can be computed without any homomorphic multiplication.

p=5. f (X, Y ) = 4X 2 + 4X + Z .

Two non-scalar multiplications are needed to compute X 2 and Z. p=7.

f (X, Y ) = 1 + 4X(X + 1) + 6[X(X + 1)] 2 + (X 2 + 3X)Z + 6Z 2

In this case, four non-scalar multiplications are needed (indicated in bold) when rewritten as follows

f (X, Y ) = 1 + 2(X 2 + X) • [2 + 3(X 2 + X)] + Z • [(X 2 + 3X) + 6Z] .

C Proof of Theorem 3

Let Z = X -Y . Thus we can rewrite Q LT S (X, Y ) as the univariate function χ F - p , namely

Q LT S (X, Y ) = χ F - p (Z) = -1 a=-p-1 2 1 -(Z -a) p-1 .
Thanks to Lemma 3, we can expand (Z -a) p-1 and obtain

-1 a=-p-1 2 1 - p-1 i=0 Z i a p-1-i = p-1 i=1 Z i -1 a=-p-1 2
(-a p-1-i ).

If i is even and i < p -1, then the ith coefficient is equal to

- -1 a=-p-1 2 a p-1-i = - p-1 2 a=1 a p-1-i = - 1 2 p-1 2 a=-p-1 2 a p-1-i .
This coefficient is equal to 0 for any even 0 < i < p -1 thanks to Lemma 4. The (p -1)-th coefficient is equal to -(p -1)/2 = (p + 1)/2 mod p. If i is odd, then we can rewrite the ith coefficient in the following way

- -1 a=-p-1 2 a p-1-i = p-1 2
a=1 a p-1-i , which finishes the proof.

- 1 ]

 1 . If input integers a and b belong to S, the result of LT S (a,b) is computed with its interpolation polynomial over F p . If a, b are larger, e.g. a, b ∈ [0, p d -1

  When the size of input integers exceeds |S| d , we can decompose integers in base |S| d and then compare their vectors of digits using the lexicographical order < on F d p l , for some d ≤ d. In fact, we compute two lexicographical orders on top of each other. Let a, b ∈ [0, |S| ld -1] be input integers. We represent an integer a = l i=0 a i |S| id by the vector a = (a 0 , a 1 , . . . , a l-1 ) ∈ F d p l of its digits of length l. The comparison of two integers a and b is thus equivalent to the comparison of their vector of digits a = (a 0 , a 1 , . . . , a l-1 ) and b = (b 0 , b 1 , . . . , b l-1 ) using the lexicographical order < on F d p l defined as follows

  , a1, a2, a3)

Lemma 3 .

 3 For all (a, b) ∈ Z 2 we have:(a -b) p-1 = p-1 i=0 a i b p-1-i mod p.Proof. Using the binomial theorem we obtain(a -b) p-1 = p -1 -i)! = (p -1)(p -2) . . . (i + 1) 1 • 2 . . . (p -(1 + i))

X 1 - 1 = 1 -

 111 i Y j a p-1-i b p-1-j .Let us compute individual polynomial coefficients. The constant term is equal to ap-1 -b p-1 + a p-1 b p-a p-1 -1 + a p-1 = 0 .Coefficients by X i with i > 0 can be computed as follows

-

(- 1 )

 1 δie e + 1 i B e+1-i • n i ,where δ ie is the Kronecker delta and B i is the ith Bernoulli number. This implies that there exist a polynomial P (X) ∈ F p [X] of degree e + 1 such that n k=1 k e = P (n).

Table 1 :

 1 Less-than function via different methods. The running time of our lexicographic order algorithms and the algorithm of Tan et al.[START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF] to compare 64-bit integers with encryption parameters supporting λ bits of security. The second column (Type) indicates which comparison circuit is used: the univariate (U), bivariate from this work (B) or bivariate one from[START_REF] Tan | Efficient private comparison queries over encrypted databases using fully homomorphic encryption with finite fields[END_REF]. The total time is averaged over 50 trials.

	Total	Amortized time
	time, s	per slot, ms

Table 2 :

 2 Less-than function via the univariate circuit. The best empirical running time of our lexicographic order algorithm to compare 64-bit integers with encryption parameters supporting λ bits of security. The algorithm is based on the univariate circuit. The total time is averaged over 50 trials. computational overhead and ciphertext noise, thus leading to larger encryption parameters.Table2shows that for p = 131, the univariate circuit takes only 11 milliseconds to compare two 64-bit integers, which is more than 3 times faster than the best running time achieved by the circuit of Tan et al.

	N log 2 q #Trials Average total	Amortized time	Amortized time
				time, s	per slot, ms	per slot, ms [12]
	8-bit integers (d = 2, l = 1, k = 9352)
	4	589	32	186.28	20	140
	8	599	20	867.46	93	690
	16 599	20	3652.23	391	3140
	32 599	20	14769.23	1579	13900
	64 604	10	60351.02	6453	60000
	32-bit integers (d = 3, l = 2, k = 4676)
	4	659	20	299.17	64	200
	8	671	20	1356.19	290	944
	16 671	20	5700.12	1219	4280
	32 684	20	23017.03	4922	18600
	64 684	10	89972.27	19241	49700

Table 3 :

 3 Sorting. The running time needed to sort N 8-bit or 32-bit integers with p = 167, m = 28057 and n = 28056. The minimal security level is 92 bits according to the LWE estimator

Table 4 :

 4 Array minimum. The running time needed to find the minimum of N 8-bit or 32-bit integers with p = 17, m = 41761 and n = 41760. The parameter T denotes the number of the tournament method stages. The minimal security level is 121 bits according to the LWE estimator[START_REF] Martin R Albrecht | On the concrete hardness of learning with errors[END_REF].

Table 5 :

 5 Comparison with other HE schemes. Total and amortized running time of the less-than function implemented with different HE schemes including TFHE, CKKS and BGV. Note that TFHE does not support SIMD packing,