Faster Polynomial Evaluations for SIMD FHEs and Application to BGV in HElib (CRYPTO 2026)

Okay, uh good morning, everyone. Today, I will present our Crypto 2026 paper, Faster Polynomial Evaluations for SIMD FHE. The core of this work is to reduce the cost of evaluating high-degree polynomials in SIMD fully homomorphic encryption, with a particular focus on BGV bootstrapping. The talk has four parts. First, I will motivate the problem by reviewing FHE schemes, bootstrapping, and the high cost of homomorphic multiplication. Second, I will discuss the algorithms. I will begin with the general-purpose Paterson-Stockmeyer method, review prior norm- and trace-based methods, and then present our own methods. Third, I will show the experimental results from our ACD implementation. Finally, I will close with two open questions. Let us >> [snorts] >> start with the basic motivation. Why does polynomial evaluation matter so much in FHE? Uh FHE is powerful because it allows us to compute an encrypted data without learning the underlying plaintext. A ciphertext contains not only the encoded plaintext, but also noise. Every homomorphic operation increases the noise and consumes part of the available noise capacity. When the noise goes over the limit, the ciphertext cannot be decrypted correctly. You can see these two diagrams. On the left, the noise after the operation still remains within the available capacity. So, the ciphertext can be decrypted correctly. On the right, the noise exceed the capacity boundary. And correct the decryption is no longer guaranteed. The basic operation of FHE are mainly multiplications and additions. For the second generation FHE schemes discussed here, especially BGV and BFV, ciphertext multiplication is far more expensive than addition. It takes substantially more time, and it also cause a much larger increase in ciphertext noise. As homomorphic computation is built from additions and multiplications, more complicated functions are represented or approximated by polynomials. An especially important example appears inside bootstrapping. The pipeline shows the four main stages of bootstrapping. Among which, digit extraction is the dominant runtime bottleneck. It homomorphically evaluates a high-degree polynomial that approximates the mod P operation required by bootstrapping. Therefore, a faster method for polynomial evaluation uh cannot only save time on computing tasks, but also directly improve the performance of FHE schemes itself. So, now let's see what method do we have to evaluate a polynomial. The Paterson-Stockmeyer algorithm is a general purpose method for polynomial evaluation. And it is widely used in FHE. The basic idea is to split a polynomial into blocks and reuses common set of small powers often called baby steps. The slides shows concrete example with target degree equal to 15 and block size equal to four. We write f of t as f0 of t plus t to the fourth times f1 of t plus t to the eighth times f2 of t plus t to the 12th times f3 of t. Each block polynomial is cubic and it reuses the same baby set basis. Means 1 t t squared and t cubed. For the outer combination we precompute t to the fourth and t to the eighth. We adopt binary strategy instead of owner masters. The advantage of this is that we can achieve multiplication depth of log d. With an appropriate block size PS uses o of square root d scalar multiplications. Its advantage is that it works for essentially any polynomial. But we were wondering if it could be more efficient. Prior work by Okada and co-authors exploited the Galois ring structure of BGV plaintext slot. And introduced the norm and trace-based methods. Uh these methods achieve all of log d several times the multiplications, uh but only when two conditions hold at the same time. Uh first is the slot degree must be a power of two, and second is the target polynomial degree must not exceed slot degree. Uh the table here shows why this regime is too narrow for our practical setting. Uh you can see none of these five slot degrees are power of two, and all the corresponding target degree are larger than the slot degree. So, in all five cases, both assumptions fail. We therefore uh need two extensions. One is an orbit product algorithms for arbitrary slot degree, and the other is method for the case target degree greater than slot degree. Now, let's see how the norm method works, since it is an important component of our algorithm. Uh given an arbitrary monic target polynomial f, where the target degree is equal to the slot degree, we decompose f into the sum of two polynomials. The first one is uh monic basic irreducible polynomial h of degree equals to target degree. And the other one is a correction polynomial g of degree less than target degree. The support of g contains only the constant term and the monomials whose exponents are powers of two. Uh therefore, T of T can be evaluated efficiently using only repeated squarings and additions. Uh for the monic basic irreducible polynomial H, we observe that we can find a root alpha in the slot. Uh then, for point T, the value H of T can be expressed as the slot norm of alpha minus T. So, what we need to compute is the norm of alpha minus T. By the definition, uh the norm is the product of all the conjugates of alpha minus T under the action of the Galois group. Since the Galois group is cyclic and generated by pi, its elements can be rewritten as pi to the zero, pi to the one, and so on up to pi to the D minus one. Therefore, the norm can be expressed as the product of the results obtained by applying pi to the zero through pi to the D minus one to alpha minus T. When the slot degree is a power of two, uh this product is easy to compute. We define A sub zero as alpha minus T and A sub J as the product of the first two to the J conjugates of alpha minus T. Then, we observe that A sub J plus one is equal to A sub J times two pi to the two to the J applied to A sub J. Therefore, starting from A sub zero, we can recursively compute A sub L in log D steps where D equals to two to the L and A sub L is exactly the product of all D conjugates, namely the norm of alpha minus T. Now, let's see how we handle an arbitrary slot degree. First, we use D to denote the product of size D consecutive conjugates of X under pi. Then, we split the size D using Euclidean division. D equals Q times D prime plus R. So, we have D prime full blocks of size Q plus one remainder block of size R. We first compute one full block of size Q and denote its product by B. Similarly, we denote the product of the remainder block of size R by C. If there is no remainder, we simply set C to one. Now comes to the key step. We regroup the full blocks. Each full block is just a shifted copy of B obtained by applying a power of pi to the Q. Therefore, the product of all full blocks becomes a new orbit product of size D prime with generator P to the Q. And then, we multiply the by the shifted remainder block. So, one problem of size D is reduced to three smaller problems of size D prime, Q, and R. And then we apply the same procedure recursively. We use the following simple rule to choose the value of Q. If D is even, set Q equal to two, then D prime is D over two, and the remainder is zero. Otherwise, if three divides D, set Q equal to three, then D prime is D over three, and the remainder is again zero. In every other case, set K equal to two, take D prime as the floor of D over two, and the remainder is one. Under this rule, Q is always two or three, and R is always one or zero. These are constant size subproblems. Only D prime continues recursively, and D prime is at most half of D at each level. The multiplication counts per level is also explicit. If Q is two and R is zero, forming B costs one multiplication, and no final combination is needed. So, the total is one. If Q is three and R is zero, forming B costs two multiplications. Again, with no final combination. If Q is two and R is one, forming B costs one multiplication, and then incorporating the remainder costs one more, for a total of two. The two possible costs do not accumulate to three. So, every level uses at most two homomorphic ciphertext multiplications. The previous slide immediately gives the recurrence for multiplicative complexity. If T of D denotes the number of ciphertext multiplications, then T of D is at most T of D over two plus two. Each recursive step reduces the problem size by at least a factor of two. So, after at most the ceiling of log base two of D levels, we reach the base case. Each level contributes no more than two multiplications. Therefore, T of D is at most two times the ceiling of log base two of D, which is O of log D. This is the first main result. The four orbits product returns logarithmic multiplicative complexity even when the slot degree is not a power of two. We now turn to the second limitation. The target degree may be larger than the slot degree. In the regime considered here, the target degree can extend to lower D log P. We decompose f of x as a of x multiplied by the product of h sub one through h sub k plus a correction g of x. Each h sub i is a monic basically reducible polynomial of degree D. The residual factor a is an arbitrary monic polynomial over integers modulo p to the e, and its degree is upper D minus K lower D, which is strictly smaller than lower D. The most important remaining task is to evaluate the product of K degree lower D factors efficiently. That is where we combine Pattersons-Zuckermeyer with the orbit product schedule. >> Let alpha sub J denote a root associated with the J's degree D factor. The structured product upper H of T is initially written as -1 to the KD times a product over J and inside it a product over the D Frobenius conjugates. The key observation is that we can exchange the order of these two products. We first compute piece of K defined as the product of alpha sub J minus T for J from 0 to K minus 1. This is a degree K product. So, the Pattersons-Zuckermeyer method evaluates it using O of square root K multiplications. Once piece of K is available, we apply the generalized orbit product schedule to it. By the arbitrary D results from the previous slides, this stage uses O of log D multiplications. So, the total multiplicative complexity is therefore O of square root K plus log D. As for the other components, we deal with them using the same method as norm based method. So, now we break the second limitation and we can evaluate a target polynomial with degree more than slot degree with multiplicative complexity of low o log square root k. Now, we turn to the experimental results. We implemented the method in ACD and compared it with the work of Ma and other co-authors published in Eurocrypt 2024. To make the comparison controlled, all cryptographic parameters are kept identical. The only component that changes in the strategy used to evaluate the digit extraction polynomial. The table lists five complete parameter sets. Across these sets, slot degree ranges from 14 to 45 while target degree ranges from 81 to 1,035. The new method accelerates digit extraction in all five parameter sets. The speedups are from 1.22 times to 2.16 times. These values are based on average wall clock time over 100 trials. The best case is set two. The baseline from Ma and co-authors takes 1,005.94 seconds while our method takes 465.39 seconds. This is a 2.16 times speedup. The speedup is not free. Four of the five parameters sets use greater multiplicative depth and consume more capacity. To evaluate the overall trade-off, we use the bootstrapping efficiency metric cap over time. This is the remaining capacity of turbo strapping divided by the bootstrapping wall clock time. Higher is better because it means we return more usable capacity per unit of time. Cap over time improves for all five sets. The improvement factors are 1.20 times, 1.51 times, 1.22 times, Oh, sorry. 1.26 times, 1.22 times, and 1.41 times for set one through set five. Therefore, although the new strategy increases depth in four cases, the reduction in run time is large enough that overall bootstrapping efficiency improves in every cases. I will close with two open questions. The first is theoretical. Our structured decomposition relies on heuristically finding a sparse correction G such that the resulting polynomial H is a basically reducible. Can the existence of such a decomposition be proved formally? A proof would replace empirical evidence with an algebraic guarantee. The second question is algorithmic. Can the search for the sparse correction G scale to larger K without prohibitive offline cost? The current brute force search can become a pre-processing bottleneck. So, a more scalable construction would be important in practice. To summarize, we generalized the orbit product schedule from power of two slot degree to arbitrary slot degree while preserving all of log D multiplicative complexity. For higher degree polynomials, sweeping the products gives all of square root K multiplications. In ASH Lib, the resulting digit extraction strategy delivers speed-ups from 1.22 In ASH Lib, the resulting digit extraction strategy delivers speed-ups from 1.22 to 2.16 times and improves cap over time in all five parameter sets. Closing the gap between the existence of a structured composition and its efficient construction would complete the picture. Thank you for your attention. I would be happy to take questions.