Integer polynomial recovery from outputs and its application to cryptanalysis of a protocol for secure sorting

3.1 k -Nearest Neighbor protocol

k -NN algorithm is a basic method used in data mining, machine learning, and pattern recognition used for the classification [10,11]. k -NN is used to search, as per a given metric, k neighbors closest to a given Δ -tuple query point in a database containing n many Δ -tuples. In the privacy preserving version of the k -NN algorithm, the query point is an encrypted Δ -tuple and the database contains n encrypted Δ -tuples. Many efficient solutions have been proposed for determining k -NN on encrypted data [5,12,13, 14], to name a few. Among these, [5] is a solution in a noncolluding federated two-cloud setting more described in Section 3.3.

3.2 Computation on encrypted data

Homomorphic encryption schemes enable computation on encrypted data such that the result of computation is available only after decryption. Current fully/somewhat homomorphic encryption (F/SHE) [15,16,17] schemes come with a price in terms of the large size of ciphertexts and are inefficient when high multiplicative depth circuits like search and sort operations are involved, and hence are impractical when handling a few hundred data elements [18,19,20]. One way to overcome this is a solution that employs two servers, Server A where encrypted data are stored and Server B that has keys needed to decrypt and perform efficient plaintext computations. Server A does some basic operations on encrypted data followed by some transformation so as to “obscure” it before sharing it with Server B , Server B then decrypts the received information, performs operations on plaintext data efficiently, re-encrypts the results and shares it with Server A to give the processed data back to the end-user client. In these models, it needs to be ensured that Servers A and B do not collude with one another. Moreover, protocols between them should ensure that Server A obscures the data in such a way that Server B is provided only as much data as required for the required computation and in a form from which Server B cannot glean anything more. There are many examples of noncolluding multiserver solutions in the literature deployed in disparate settings, and we present a couple of example scenarios. Aono et al. [21] considered a two-cloud model to compute logistic regression securely, and a solution by Hardy et al. [22] used a three-cloud federated setting to do learning on vertically partitioned data, which are examples in a machine learning environment; while pRide [23] and PSRide [24] are two solutions in a ride-hailing environment that use two-cloud models.

3.3 k -NN protocol over encrypted data [5]

At EDBT 2018, Kesarwani et al. [5] presented a new secure k -NN protocol in the two-party honest-but-curious federated cloud model for computing k -NN on encrypted data. Their solution uses a semantically secure (levelled) FHE (LFHE) [25] to compute squared euclidean distances (EDs) directly on encrypted data. To compute the ranking among the distances, the distances are suitably transformed to preserve the order. A federated noncolluding public cloud having the secret key performs the comparison using the transformed data. Since this cloud has access only to transformed results of computations performed on plaintext data, this article claims that public cloud does not learn anything useful about the original database, the results, or the query. The federated cloud setting consisting of two Servers A and B is depicted in Figure 1.

Figure 1

k -NN setting from [5].

This article claims to provide an asymptotically faster solution than the current state-of-the-art protocol. They implement their protocol and run experiments on two real-world datasets from UCI machine learning repository [26], one from healthcare domain and another from financial domain having a relatively large number of dimensions (32 and 23, respectively) and containing 858 and 30,000 data points, respectively; both of which contain personally identifiable sensitive information like gender, age, marital status, etc. This article states that the data were preprocessed to remove negative integer values. They run their experiments for k = 2 , 8, and 16. While the datatype in each dimension is not explicitly mentioned in this article, by examining the datasets, we find that they are nonnegative integers of size ≤ 16 bits.

This article makes the following security guarantee regarding the leakage profile for the Server B as given in ref. [5, Theorem 4.2], reproduced here for reference.

Theorem 4.2. Secure k -NN guarantee: Server B : Our secure k -NN protocol leaks no information to Server B except the number of nearest neighbors k to be returned by the protocol and the number of equidistant points in the database with respect to a given query point q .

It may be noted that though the protocol of [5] has been described specifically in the context of securely evaluating k -NN, their technique of transforming through a random monotonic polynomial has applications in many settings where sorting of SHE or LFHE encrypted data is needed. We remark here that if sorting is the only functionality required, then order-preserving or order-revealing encryption schemes would suffice for the purpose [27,28,29].

3.4 Problem setting

We refer to Figure 1 where the data owner outsources his/her database in an encrypted form using a semantically secure homomorphic encryption scheme for storage in Server A . Each data item in the database is of dimension Δ . Server A does not have access to any secret keys and operates only on encrypted data. Server A provides storage for the database and provides services on the encrypted database homomorphically. One of these services is the computation of the k -NN of a given query point q , of dimension Δ , received from an end-user or client. Server A homomorphically computes squared EDs between the query point and each of the n data points in the database. Because a single ED 2 computation is of multiplicative depth 1, it can be efficiently evaluated using an LFHE scheme. Since Server A does not possess the decryption key, it will not be able to efficiently uncover the underlying plaintext information of either the query point or the data points.

Computing ED is of low multiplicative depth, but sorting is not. Hence, Server A transforms the distances while preserving its order. It picks a monotonic polynomial p ( x ) of degree d of the form a 0 + a 1 ⋅ x + a 2 ⋅ x 2 + … + a d ⋅ x d for some chosen d ∈ N , where each of the integer coefficients a i are picked uniform randomly and independently in the range [ 1 , 2 α − 1 ] , for example, in the range [ 1 , 2 32 − 1 ] as presented in ref. [5, Section 3.4]. This polynomial is then evaluated homomorphically for each of the EDs and the output ciphertexts are re-ordered using a permutation σ picked uniformly at random, before sending them to the Server B for sorting.

Server B possesses the decryption key using which it will decrypt the values received from Server A and sorts them. As the decrypted values are outputs of a random polynomial, the original distances as computed by Server A are “hidden” from Server B . Server B then sends the indices of k -NNs to Server A , which then applies σ − 1 to the received ordering of the indices and forwards the same to the client. We refer to the security guarantee with respect to Server B given in [5, Theorem 4.2] in Section 3.3 and ask the question whether Server B can glean more information about the polynomial chosen by Server A and, eventually, the plaintext integer inputs used by Server A for its homomorphic polynomial evaluation.

Remark: The plaintext data points and the query point are encoded as tuples of integers, since in the context of F/SHE schemes, fixed-point values too are (exactly) encoded using essentially the scaled-integer representation [30]. The only exception is the CKKS FHE scheme [31] that natively supports floating-point arithmetic but was not considered in ref. [5]. However, it may be noted that a successful key recovery attack on CKKS was published in EUROCRYPT 2021 [32]. We leave it for future work to extend our attack to also include any implementation based on the CKKS scheme.

3.5 Adversarial model

Servers A and B are assumed to be honest but curious. Each will perform the computations correctly but is keen to learn more about the distances between the query point and the data set points. Server A has access only to encrypted data. The semantic security of the underlying encryption scheme will guarantee the protocol is secure against Server A .

Using the decryption key, Server B decrypts the permuted results of computation performed by Server A . After decryption, the Server B would observe only the outputs of the polynomial evaluation (and not the input squared distances). That is, it only sees the values on the L.H.S. of the following set of equations:

It is assumed that the Server B , seen here as the adversary, knows the degree d , which is usually small since the homomorphic evaluation of polynomials in encrypted form are efficient only for small degrees. It also knows the range [ 1 , 2 α − 1 ] for the unknown coefficients a i , and the range [ 0 , 2 β − 1 ] for the unknown inputs x i . For our attack, we need not know the exact values for the three parameters, just an upper bound on them would suffice. Also, note that all the parameters above take nonnegative integer values.

As noted earlier, Server A homomorphically evaluates the polynomial p ( x ) at n (squared ED) integer values x 1 , … , x n , and we can assume without loss of generality that 0 ≤ x 1 ≤ x 2 ≤ … ≤ x n . Since p ( x ) is monotonic, the ordering of p ( x i ) is identical to the ordering of x i . If a 0 ≤ x 1 , then for any given tuple of coefficients ( a 0 , a 1 , … , a d ) , there will be a set of positive real roots ( χ 1 , χ 2 , … , χ d ) to (1). Hence, the authors of ref. [5] seem to argue that if the range for a i is large enough, then it will be infeasible to search for all possible polynomials satisfying (1). The authors claim that the probability that Server B successfully recovers the coefficients a i , followed by x i , is approximately 1 / 2 α ⋅ ( d + 1 ) , which is negligible when α is large. Referring to the example given in ref. [5, Section 4.2], for α = 16 and d = 9, this probability is approximately 2 − 160 , which is negligible. Hence, the protocol leaks only a negligible amount of information, as claimed in ref. [5], about either a i or x i to the Server B and nothing else other than the ordering of the x i ’s. We note here that the x i values may never be uniquely recovered in (1) with probability = 1 since p ( x + c ) is also an equivalent polynomial satisfying the equation for c ∈ Z , and there may be many values of c such that 0 ≤ x i − c < 2 β . Hence, the inputs and the polynomial may only be recovered up to an affine scaling. Other nonlinear transformations may also result in an equivalent solution. For instance, p ( x ) can be a potential solution when all the x i are perfect squares and p ( x ) contains only even powers. But these possibilities will likely be significantly small when n ≫ d , and the coefficients a i and/or the inputs x i come from a random source with sufficiently high entropy, which indeed is the scenario in ref. [5]. Please refer to Section 4.10 where we describe the uniqueness of the obtained solution.

3.6 Our contribution

We give an algorithm (see Algorithm 1 on p. 36) to the earlier defined polynomial recovery problem where the goal is to recover the integer inputs (up to an affine scaling) by observing only the outputs of the monotonic polynomial with positive integer coefficients, assuming the number of evaluation points is much greater in number compared to the degree of the polynomial. Sections 4.1–4.7 explain the core idea of our proposed algorithm. Our algorithm has a heuristic expected time complexity that is exponential in the size of the inputs ( β ) and the degree ( d ) of the chosen polynomial, but polynomially dependent on the size of the coefficients ( α ). As derived in equation (11) on p. 11, the heuristic expected time complexity of our algorithm is O ˜ ( d 2 β 2 β ( α + d β ) + n d 3 ( α + d β ) d ) . We would like to note that in many real-world scenarios the inputs are or can be encoded as integers of 16- or 32-bit length, and our method runs efficiently for inputs of such size. Also note that in SHE applications d is required to be not large as well. Since there can be many solutions to the aforementioned polynomial reconstruction problem, hence we will output one solution that satisfies all the output points (there is also a possibility to enumerate all the solutions). But as discussed earlier, when the number of output values is far bigger than the input degree, and the coefficients and/or inputs are sufficiently random, then the number of equivalent solutions will likely be small. Indeed, we show in Section 4.10 that when the polynomial coefficients are chosen uniform randomly, the original polynomial along with the inputs are recovered by our algorithm with a high probability. Our algorithm extends to recovering any integer polynomial (not necessarily a monotonic integer polynomial) and any input range (not necessarily [ 0 , 2 β − 1 ] ), as long as the range of the coefficients and inputs are known apriori.

The aforementioned results invalidate the security claim in ref. [5, Theorem 4.2] regarding the leakage profile for the Server B . In particular, Server B will be able to learn the square of the EDs between the query point and the data set points. It may not be able tell the exact distance to a given point due to random re-ordering but will be able to know all such values. Such an information can potentially help the adversary to narrow down further if it has access to extra information about the underlying data set or query point. For the concrete parameters suggested in ref. [5, Section 4.2], i.e., β = 16 and d = 9 , we can recover the inputs for α = 16 in a few seconds (see Table 1 on p. 29). We then extend it further up to d = 96 , α = 128 , and β = 48 and show that we can recover the polynomial and the inputs efficiently (see Table 2 on p. 28). Finally, we investigate in Section 6 another variant of the protocol of [5] where the (homomorphically) transformed polynomial outputs are perturbed by a noise yet maintaining the monotonicity. In this case, our previously mentioned attack will not work. But we show, in a straightforward way, that it is still possible to recover the ratios of the inputs.

A preliminary version of our work has appeared as part of 17th IMA International Conference of Cryptography and Coding (IMACC’17), 2019 [33]. We have made a number of improvements to our algorithm (and the material in this article has more than 50% new material compared to the earlier version). We have also run our experiments with much larger bounds. We list hereunder some of the important changes and enhancements that are present as part of this submission:

1. The heuristic expected running time of our earlier polynomial recovery algorithm (found in the conference proceedings) was O ˜ ( n 2 β ( α + d β ) d ) , where n is the number of polynomial integer outputs, the parameters d and α denote the degree and the size of coefficients of the original polynomial, and β denotes the size of the original inputs. This earlier algorithm used a brute force approach to iterate over the input space to try and find a suitable polynomial in the second step.

2. In the current version, we have improved the algorithm significantly and the heuristic expected running time of the polynomial recovery algorithm is now O ˜ ( d 2 β 2 β ( α + d β ) + n d 3 ( α + d β ) d ) . This is made possible by avoiding the brute force search mentioned earlier.

3. In addition to running tests with parameters that were claimed secure in the article by Kesarwani et al. mentioned earlier, we show that our algorithm can recover the polynomial and the inputs for larger bounds on the size of inputs, size of coefficients, and the degree of the polynomial. In the earlier version, input values were in the range [ 0 , 2 24 − 1 ] , polynomials up to degree 9, and coefficient space was [ 0 , 2 32 − 1 ] . In this version, the range of the parameters have been greatly increased; the range of input values is now [ 0 , 2 64 − 1 ] , polynomials up to degree 96, and the polynomial coefficient space is now [ 0 , 2 128 − 1 ] .

4. When the coefficients are chosen uniformly randomly, we analyze and show that the recovered polynomial and the input values are, with a very high probability, the same as what was chosen initially.

The earlier version of our polynomial recovery algorithm has also found applications in the privacy-preserving ride-hailing Service (RHS) scenario [34] to recover driver locations that are homomorphically transformed by a service provider as a means to protect drivers’ privacy.

4.1 Divisors from y -differences

Summary of this step: Given a nondecreasing order of n outputs y i = p ( x i ) with y 1 being the smallest, the compute y -differences y i , j = y i − y j , where i > j , followed by finding divisors of y i , j , of which one divisor is equal to x i − x j . Differences of y values are stored in Matrix Y , and the corresponding divisors are stored in Matrix D at the end of this step.

Refer to the procedure GuessTheDifference from Algorithm 1.

Let y i , j = y i − y j . From Lemmas 1 and 2, ( x i − x j ) < 2 β is a divisor of y i , j , and as part of this step, we find all positive divisors of y i , j < 2 β . For β ≤ 28 , we use the sieve method to obtain factors by dividing y i , j by successive primes and forming positive divisors < 2 β . When β > 28 , and since we know that there will exist a small factor (as in our case, since the range of x i are known, and the factors are differences of x i ), we use elliptic-curve factorization method (ECM) [35] to find the required factors.

Time complexity for finding factors. In the expected case, ECM is a subexponential algorithm with complexity O ( e ( log p log log p ) ( 1 + O ( 1 ) ) ) to find a single factor, where p is the smallest factor. We first use the sieve method to obtain all small factors ≤ 2 28 and then invoke SageMath ECM.find_factor() on Q , where Q is the quotient obtained by dividing y i , j by all the small factors and their multiplicity, to obtain the next possible factor f . If f < 2 β , we continue the process of invoking ECM.find_factor() on Q / f as long as we obtain a factor < 2 β or we hit a timeout. The timeout value is empirically set such that the aforementioned function is able to find “small” factors ( < 2 β ) well within the timeout.

Remark. There is a possibility of losing factors due to our timeout being insufficient. Section 5 presents results on loss of factors for different timeout values. But we would like to stress that we must enumerate all the factors of the differences that are less than 2 β if we want to eliminate any risk of missing any satisfying polynomial.

Let D i , j be the set of all divisors of y i , j . This set constitutes the guesses for the differences of the (to be determined) values x i . It turns out that for many values of y i , j , there may be too many divisors that are < 2 β , so we need to sample larger number of output values (i.e., larger n ) and carefully pick d number of y i , j ’s such that the number of elements in each of D i , j is a small positive number (say, ≤ ψ ), whereby the search space for the guesses becomes feasible to enumerate. We also store the corresponding y i , j values in Matrix Y that are to be used later in the Lagrange interpolation step .

Remark. In our experiments as well as in our analysis, we consider ψ = O ( α + d β ) , since the number of divisors of an integer N is bounded by N O 1 log log N and on the average case, it is log N [36].

The output of this step is Matrix D of divisor sets:

Time complexity. The worst-case time complexity of this step (in terms of bit operations) to find the factors less than 2 β is

since the size of each Y i , j is O ( α + d β ) , and there are O ( d 2 ) elements in D i , j .

4.2 Consistency check on the guessed differences

Summary of this step: The key idea in this step is to eliminate as many divisors as possible before performing Lagrange Interpolation in the following step.

Refer to the procedure CheckConsistencyInColumn from Algorithm 1 in Appendix A.

Remark. If x i − x k is a divisor of y i − y k , and x j − x k is a divisor of y j − y k , then x i − x j must be a divisor of y i − y j .

The output of the previous step is the matrix D . Element D i , j ∈ D , is a list of divisors of y i , j , and we know that only one of them is (but as yet unknown) x i − x j . We create sets of consistent divisors for each D i , j ( 2 ≤ i ≤ d + 1 ) and ( j < i ) . All elements that are not part of the consistent set are set equal to 0.

For each nonzero element in D 2 , 1 , pick its corresponding “consistent” element δ i − 1 from D i , 1 , 2 < i ≤ d + 1 to form a tuple of d + 1 elements ( δ 0 = 0 , δ 1 , … , δ d ) . Store each tuple in an array of tuples C . The output of this step is Array C and ∣ C ∣ .

Remark. In the beginning of Section 4, we use the notation δ i to mean the x i − x 1 differences, but in this section, we use the same notation to mean “consistent” elements. This is because they represent the same element. In other words, a “consistent” element is one that is a divisor of the y output difference as well as (yet to be identified) x input difference.

Time complexity. The worst-case time complexity of this step is

as O ( d 2 ) elements each having O ( ψ ) factors of size O ( β ) bits are compared in O ( d ) columns.

4.3 Recursive consistency check

Summary of this step: In Section 4.2, we obtain the consistency set by examining the divisors of differences of the polynomial outputs and using that we obtain a set of tuples that constitute the sets of guessed differences of inputs. In this step, the differences are obtained in a different manner as described later in this section. These two sets of guessed differences are compared to enable us to converge faster on the set of x differences.

Remark. In the worst case, we need to run this step for ψ d number of times, so we enable this procedure only for experiments with small polynomial degree. In practice, this step quickly reduces the number of inconsistent tuples in Array C , thus reducing the overall running time of the algorithm.

Refer to the procedure RecursiveConsistencyCheck from Algorithm 1 in Appendix A.

This step takes the tuple Array C as input. The d elements of the each tuple of Array C, for example, C [ i ] [ 1 ] to C [ i ] [ d ] of the ith tuple are factors of y i output differences, as obtained the previous step. By using the d elements of the first tuple, we can rewrite the first column of Matrix Y as follows:

We start with Column m = 1 of the Matrix Y , and rename y i outputs as z i , 0 (for ease of notation) for 1 ≤ i ≤ d + 1 . In other words, y 1 is renamed z 1 , 0 , y 2 is renamed z 2 , 0 , and so on. In Step 1 of this procedure, we compute quotients [ ⋅ ] 2 , 1 = z 2 , 0 − z 1 , 0 x 2 − x 1 , [ ⋅ ] 3 , 1 = z 3 , 0 − z 1 , 0 x 3 − x 1 , and so on. Generalizing this, in step m , we compute the quotients [ ⋅ ] i , m , for ( m + 1 ) ≤ i ≤ ( d + 1 ) , as follows:

Next, we find the differences z i , m − z m + 1 , m for ( m + 1 ) < i ≤ ( d + 1 ) . We then find divisors ( < 2 β ) of each of these differences. We call this divisor set as D N i , m + 1 . In Matrix D , we set D i , m + 1 = D N i , m + 1 ∩ D i , m + 1 . We then build the consistent set for Column m similar to the way described in Section 4.2. Finally, we pick the nonzero element in Column m in every row to form a tuple τ of ( d − m ) number of elements, which form the guessed difference divisors for the next recursive step. The recursion terminates in d + 1 steps as given in Lemma 3. In each step, the computed τ is used as input to the subsequent call.

As shown in Lemma 3, for a single tuple of guessed differences, this step terminates in at most d steps. If no a d is returned, then the tuple in Array C is discarded, and the procedure is repeated using the next tuple in Array C .

When recursion completes successfully, we obtain

We obtain z d + 1 , d = a d as given in Lemma 3. We note that a tuple in Array C having values with successive differences matching that of the original x i values will successfully go through all the d + 1 steps until an integer a d is obtained.

Time complexity. In this step, the factors are obtained for O ( d 2 ) elements followed by running consistency check on d − 1 columns. For one tuple, the worst-case time complexity of this step based on (3) and (4) is

4.4 Choosing optimal divisors sets

Summary of this step: In Section 4.2, we enumerate over ψ d many tuples of divisors/guesses. We provide here a way to choose the divisor sets such that the enumeration complexity is as small as possible.

We refer to equation (1) giving the polynomial outputs. As explained in Section 4.1, we compute the differences of polynomial outputs and form the Matrix D of divisors sets (less than 2 β ) of the output differences. We need to find a set of d many D i , j elements of the Matrix D such that ∏ D i , j is minimum. In other words, the product of the number of guesses is minimum. Graph G is the Output Difference Graph of Matrix D . Now finding the minimum ∏ D i , j is akin to finding the d -minimum spanning tree (i.e., a minimum weight tree with d edges only) in G , where the weight of the tree is represented by the product of the weights. The requirement that the subgraph is a tree comes from the linear independence requirement of the corresponding set of equations. Essentially, we are transforming the problem of finding the small search space of divisors to the problem of finding a d -minimum spanning tree having the least cost (in terms of divisor product) across all divisors sets of Matrix D yet satisfying the linear independence condition. It is shown in ref. [37] that the d -MST problem is NP-hard for points in the Euclidean plane. The same article provides an approximation algorithm to find d -MST [37, Theorem 1.2] and is reproduced here for reference.

Theorem 1.2. There is a polynomial-time algorithm that, given an undirected graph G on n nodes with nonnegative weights on its edges, and a positive integer k ≤ n , constructs a tree spanning at least k nodes of weight at most 2 k times that of a minimum-weight tree spanning any k nodes.

Note that this approximation algorithm also works for multiplication of edge weights (weights greater than 1) since by extraction of logarithms this can be trivially turned into addition of edge weights. By using this algorithm, we can carefully select d < n nodes having close to the minimum enumeration complexity to make our search space feasible to guess the differences ( x i − x j ) with x i ≥ x j . From the d -MST so obtained, we can now go on to find the set of divisors ( x i − x j ) such that they are consistent as explained in Section 4.2 and continue with finding the polynomial coefficients using Lagrange interpolation as described in Algorithm 1. However, we did not implement this optimization in our code as the concrete running time was already small enough. But for larger instances, this optimization will be useful.

4.5 Lagrange interpolation to recover all coefficients

Summary of this step: Use Lagrange interpolation to construct a degree- d polynomial using the guessed difference set.

Refer to the procedure ApplyLagrangeInterpolation from Algorithm 1 in Appendix A.

Each tuple ( δ 0 = 0 , δ 1 , … , δ d ) in the Array C obtained as described in Section 4.2 presents possible candidates for the differences ( x i − x 1 ) , 1 ≤ i ≤ d + 1 . Optionally, we could use the method described in Section 4.3 to narrow down the possible candidates to only those tuples for which an integer a d coefficient is obtained.

Since δ i represents differences with respect to some unknown value x 1 , we treat x 1 as an unknown parameter and add it to each of the d differences to obtain ( x 1 , x 2 , … , x d + 1 ) = ( x 1 , x 1 + δ 1 , … , x 1 + δ d ) . The Lagrange polynomial of degree- d is obtained using these values for inputs x , expressed in terms of x 1 , and the corresponding output values come from the y set as given in Section 4.1. We now have a degree- d polynomial L in x with parameter x 1 as shown in equation (6).

From Lemma 4, each L i is a degree- ( d − i ) polynomial in x 1 , 0 ≤ i ≤ d , and L d ( x 1 ) = a d is a constant polynomial. Thus, we have recovered a d , the coefficient of x d in p ( x ) .

We know that the coefficients of p ( x ) are positive integers less than 2 α , and p ( x ) is evaluated at integers in the range [ 0 , 2 β − 1 ] . Our next objective is to output an integer x 1 ∈ [ 0 , 2 β − 1 ] such that ∀ i = 1 , 2 , … , d , 1 ≤ x i = x 1 + g i ≤ 2 β − 1 and ∀ i = 0 , … , d − 1 , L i ( x 1 ) is an integer with 0 ≤ L i ( x 1 ) ≤ 2 α − 1 as described in Section 4.6. If no such x 1 exists, then we repeat this process using the next tuple in C .

Time complexity. The worst-case time complexity of computing the Lagrange interpolation polynomial for one tuple is:

4.6 Obtaining x 1 by finding roots of polynomials L i ( x 1 )

Summary of this step: Given α , β and polynomials L 0 , L 1 , … L d − 1 (each L i ∈ Q [ x 1 ] and is of degree d i , respectively, where d i = d − i ), find x ∈ Z satisfying the condition: 0 ≤ x ≤ 2 β − 1 and 0 ≤ L i ( x ) ≤ 2 α − 1 ∀ i = 0 , … , d − 1 .

Refer to the procedure RecoverPossibleXValuesAndPoly from Algorithm 1 in Appendix A.

It is easy to see that the following algorithm always outputs one such x if it exists (if no integer satisfies the required condition, then this algorithm returns failure).

1. For each i = 0 , … d − 1 , find the ordered set S i = { [ a k , b k ] : a k , b k ∈ Z , 0 ≤ a k ≤ b k ≤ 2 β − 1 } of disjoint intervals on the real line such that every integer m satisfying 0 ≤ m ≤ 2 β − 1 and 0 ≤ L i ( m ) ≤ 2 α − 1 is present in some interval [ a k , b k ] ∈ S i .

2. If there exists an integer x in the range [ 0 , 2 β − 1 ] that is present in some interval from each S i output it, else return failure.

Step 1. Finding the ordered set S i of intervals for each L i :

Compute the set R i , 0 of all the (approximations to the) real roots of L i ( x ) = 0 that lie in the range [ 0 , 2 β − 1 ] . Similarly let R i , α contain all the real roots of L i ( x ) = 2 α − 1 in the range [ 0 , 2 β − 1 ] . Let R i = R i , 0 ∪ R i , α be an ordered set with its elements in the increasing order (as they appear on the real line).

Intuition. We first provide an intuition for our algorithm by looking at one particular case. Let a , b (with a ≤ b ) be two consecutive elements of the ordered set R i = R i , 0 ∪ R i , α . For now, assume that a ∈ R i , 0 and b ∈ R i , 0 . Then, the following lemma easily follows.

To determine which of the aforementioned possibilities occurs, it is enough to compute L i ( c ) at any c between a and b . If 0 < L i ( c ) < 2 α − 1 , the interval [ ⌈ a ⌉ , ⌊ b ⌋ ] contains all integers between a and b satisfying 0 < L i ( ⋅ ) < 2 α − 1 . Add this interval to S i . (In our algorithm, we take the midpoint c = a + b 2 .)