Game on: a performance comparison of interpolation techniques applied to Shamir's secret sharing

Received: May 2, 2024. Revised: August 9, 2024. Accepted: October 4, 2024 © The Author(s) 2024. Published by Oxford University Press on behalf of The British Computer Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. The Computer Journal, 2024, 1–12 https://doi.org/10.1093/comjnl/bxae109 Original Article Game on: a performance comparison of interpolation techniques applied to Shamir’s secret sharing Anastassis Voudouris1,*,Aristomenis Tressos1,Apostolis Zarras2,Christos Xenakis1 1University of Piraeus, Greece 2Foundation for Research and Technology—Hellas, Greece *Corresponding author: University of Piraeus, Greece. E-mail: [email protected] Abstract Public-key encryption is typically managed through a public key infrastructure. However, it relies on a central control point, the certification authority, which acts as a single point of failure. Recent technological advancements have led to the need for decentralized cryptographic protocols. This paper presents a comprehensive study on enhancing public-key encryption via threshold cryptography and multiparty computation to ensure robust security in decentralized systems. The focus lies in exploring various polynomial interpolation techniques within Shamir’s secret sharing scheme, particularly addressing the efficiency and practicality of Newton interpolation, fast Fourier transformation (FFT), and advanced versions of Lagrange’s method. Utilizing SageMath for a dedicated testing environment, the research investigates the swiftest interpolation methods for secret recovery, introducing new shares into the system, and evaluating the impact of optimizations on performance. The findings highlight FFT as the most effective interpolation method in speed and efficiency, albeit with limitations on the number of shares that can be processed. This paper critically evaluates these interpolation techniques against practical constraints and aims to answer pivotal research questions regarding the optimal approach for large-scale scenarios, challenging existing notions on the efficiency of Newton’s method and providing experimental evidence to support the superiority of FFT in specific contexts. 1. INTRODUCTION Public-key encryption in cryptographic schemes is commonly managed through a public key infrastructure (PKI). PKI ensures secure communication across insecure public networks by employing digital signatures to authenticate entities. PKI relies on the security of a central control point, named the certification authority (CA), and is considered trusted by all; however, it also acts as a single point of failure. If this point gets compromised, the entire network’s security will fall. Recent technological advancements showed that the research community needs to strengthen current cryptographic protocols in a decentralized manner [1]. New ideas are driven by a well-studied area of cryptography called threshold cryptography [2]. A (n,k)threshold cryptography scheme allows nentities to share the ability to perform a cryptographic operation so that any kentities suffice to perform this operation jointly. In contrast, it is infeasible for at most k−1entities to do so, even by collusion. A cryptographic key, which will be used to encrypt the message transmitted from one person to another, is the fundamental tool for both symmetric and asymmetric cryptography and must be protected, from the generation part to the secret sharing among the group’s entities. Secret sharing and threshold cryptography are both branches of the same tree, called, in general terms, multiparty computation (MPC). One security application of MPC is to protect cryptographic keys from theft and misuse by never having them exposed at any single point at any single time and by enforcing usage policies at multiple entities. By embracing MPC’s ability to protect cryptographic keys in software, an alternative to legacy hardware solutions, such as Trusted Platform Module’s root of trust, is provided, mitigating many operational challenges in today’s hybrid and multi-cloud environments. Shamir’s secret sharing scheme (SSS) is the main sub-module of threshold cryptographic protocols [3]. Assuming that trust is established, a dealer owning a secret desires to distribute this secret among a group of nshareholders. A shared secret is generated by a trusted dealer, who distributes a share (of a cryptographic key, for example) to each member of the group, such that any subset of size greater than a threshold (t)can reveal or use the shared secret, while smaller subsets do not have any knowledge about it. The core subroutine of SSS is Lagrange’s interpolation, through which the shared secret is eventually reconstructed. It is noteworthy that numerous optimizations for Lagrange’s interpolation have been documented in the literature, aimed at accelerating the computation of Lagrange polynomials and obtaining the constant term of the interpolating polynomial, which serves as the secret in the standard SSS scheme [4]. Additionally, there are various other methods available for interpolation. Here are a few notable cases: multi-secret schemes involve using the secret sas the highest degree coefficient, thereby increasing the efficiency of another interpolation method called Newton’s interpolation [5]. The adoption of Newton’s interpolation has also been proposed in a recent academic publication where the authors explore the application of this Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 2|Voudouris et al. approach in the context of Multi-Factor authentication during the introduction of a new node in a communication protocol [6]. Furthermore, Tomescu et al.[ 7] aggregate a Boneh–Lynn–Shacham signature [8] 3000 times faster in 46 s compared with 1.5 days if done naively. This improvement is achieved by utilizing a Fast Lagrange interpolation algorithm, which has a time complexity of (tlog 2t), where in ntotal users, the t>n/2must be honest,and incorporates a Fast Fourier interpolation method.This approach outperforms the naive Lagrange algorithm, which has a time complexity of (t2). This centers on exploring alternative interpolation techniques instead of the traditional Lagrange method, focusing on enhancing performance. In response to concerns regarding interpolation methods within the context of SSS, this article presents a comprehensive exploration of various polynomial interpolation techniques. The methodologies scrutinized include Newton interpolation, the fast Fourier transformation (FFT), and two refinements of Lagrange’s method. The research utilizes a dedicated testing environment, employing SageMath [9], a freely available opensource mathematics software system licensed under the GPL, and utilizing Python-based language. The testing environment accommodates fields with diverse user counts and specific mathematical properties. Additionally, the authors delve into the implications of introducing new shares into the system, akin to a simple communication protocol. Experimental findings demonstrate that, in terms of performance, the FFT emerges as the most efficient interpolation method in both scenarios. However, it is noteworthy that limitations exist regarding the number of shares that can be effectively utilized. Overall, this article poses three main research questions: (i) What is the swiftest method and most suitable in large-scale scenarios to interpolate the polynomial in the context of SSS to recover the secret from the constant term? (ii) What is the quickest approach to interpolate the polynomial for recovering the secret when adding additional members/shares? (iii) Do optimizations (e.g. precomputing values) affect the overall performance of interpolation methods? These research questions form the basis of this work. In summary, this article makes the following main contributions: •A comprehensive examination of different interpolation techniques and their implementation in the SSS scheme. • A specific use case addressing the addition of new shares and the subsequent reconstruction of the secret, challenging the assertion by Bezzateev et al.[ 6] regarding Newton’s method as the swiftest. • Experimental evidence demonstrating the superior efficiency of the FFT compared with other methods. The remainder of this paper is organized as follows. In Section 2, we elaborate on several fundamental concepts in cryptography, encompassing finite fields and roots of unity, culminating in an exhaustive exploration of SSS. Moving forward, Section 3 conducts an in-depth analysis of various interpolation methods, including Lagrange’s method, its optimizations, Newton’s method, and a specialized case involving FFT. Section 4encapsulates the culmination of all the experiments. In Section 5, a concise discussion addresses the original research questions. Additionally, Section 6includes a number of research papers that employ various interpolation methods within the context of SSS, following the notable contributions of Shamir and Blakley. Finally, Section 7wraps up the paper by outlining avenues for future research and posing pertinent questions. 2. BACKGROUND Before delving into the concepts of secret sharing and polynomial interpolation, it is imperative to establish a foundational understanding. To enhance the clarity of subsequent sections in our work, we aim to lucidly present essential mathematical concepts from various fields, articulating definitions and theorems without proof. Readers are guided to the relevant references for comprehensive insight when further elucidation is required. 2.1. Mathematical foundations As a prerequisite for the fundamentals of cryptographic theory, a grasp of finite fields and groups is indispensable. The ensuing paragraphs furnish a concise introduction to vital mathematical concepts, laying the groundwork for a comprehensive understanding of our implementations and approaches to polynomial interpolation. Finite fields. The set of real numbers, R, has infinite order, i.e. contains an infinite number of elements, whereas a finite field, or Galois Field (GF), has a finite order. GF(p)(or Fp),where pis prime, is a finite field with finite order pand is mostly used in cryptographic applications due to its excellent properties. Note that the order of a finite field is always a prime por a prime power, pn[10]. Finite fields have two operations (i.e. addition and multiplication), which adopt them from the corresponding groups. A fundamental intuition of groups, especially the cyclic ones, is the generator of a group. The order of an element x∈Fp, denoted as |x|,isthe number of unique elements that xcan generate raised to the elements of the group 1, 2, 3, ..., p−1. The elements that generate all p−1elements of the group are called generators or primitive roots, and the groups generated by a generator are called cyclic groups. The order of generators is always equal to p−1.Ifanelement g∈Fpis a generator, then gk= 1∀k= p−1. 2.2. Roots of unity The following points are crucial for grasping the motivation behind the FFT and the discrete Fourier transformation (DFT), which we will demonstrate to be the fastest method in terms of time for computing the secret of SSS, under specific conditions. A complex root of unity is a convoluted number ωsuch that ωn=1. There are exactly ncomplex nth roots of unity: ei2kπ n· k=0.1. ···n−1. According to the cancelation lemma [11]we take that for any integers n≥0, k≥0,andd≥0, ωdk dn =ωk n, and that for any even integer n≥0, ω n 2 n=ω2=−1. Lastly, according to the halving lemma [11], if n≥0is even, the squares of the ncomplex nth roots of unity are the n/2complex n/2th roots of unity. The halving lemma is crucial to the divide-andconquer technique for converting between coefficient and pointvalue representations of polynomials, which we will analyze in the polynomial interpolation section, since it ensures that the recursive subproblems are half as large as the original problem. However, considering that we are working on a finite field and not on complex numbers, the definition of the roots of unity has a small twist.Let pbe a prime and Fpbe the finite field of pelements. Let ωbe an element of Fp.Wesaythatωis an nth root of unity if ωn≡1(mod p)but not ωm=1, for all 0<m<n. Readers are referred to Stewart [12] for an in-depth subject analysis. Lastly, a criterion to determine a root of unity is that the finite field Fphas a primitive nth root of unity if and only if ndivides the order of Fp, which is p−1. Primitive nth roots of unity can be a convenient way to specify elements of a group that have order n[13]. Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 Interpolation techniques applied to SSS |3 2.3. Secret sharing Secret sharing is now one of the fundamental elements of cryptographic protocols. Consider a group S=s1,s2,...,snof parties called shareholders or nodes that are mutually distrusting of one another and a centralized adversary with the ability to control up to tparties, often referred to as a dealer or data owner. Usually, the designated dealer D∈Pmay share a secret sby providing each si ashareσiof the secret svia a secret sharing scheme. On a high level, a secret sharing scheme is secure if any unauthorized subset of participants learns nothing about the message, while any authorized subset of participants can fully reconstruct the secret, also known as accessibility. SSS is perfectly secure, i.e. secure in an information-theoretic sense, meaning that unauthorized coalitions cannot deduce any information about the secret. Shamir’s secret sharing scheme. We decided to present SSS under the preliminaries section firstly because it is one of the most well-known and researched secret sharing schemes, and secondly for the readers to have an intuition of which is the end goal of this paper. Our following experiments are entirely based on this scheme. The SSS scheme is an algorithm first proposed in 1979 by the renowned Israeli cryptographer Adi Shamir [3]. Shamir’s (t,n)- threshold scheme, where tis the reconstructing threshold and n is the size of the set S=s1,...,snof shareholders, uses Lagrange interpolation. The core components of SSS scheme consist of two key algorithms: the Share algorithm and the Reconstruction algorithm. The Share algorithm takes as input a message m∈Fq,whereq is a prime number. It selects a polynomial p(x)=a0+a1x+···+ at−1xt−1of degree deg(p(x)) =t−1,wherea0:=mand coefficients a1,...,at−1∈Fqare chosen uniformly at random. It computes share σi∈Fqfor shareholder si∈Sas a point on polynomial p(x), i.e. σi:=(i,p(i)),whereiis the ID of shareholder si. It distributes share σito shareholder sithrough an information-theoretically secure channel, for i=1, ...,n. The Reconstruct algorithm takes as input a set of shares σ1,...,σr held by a subset R∈Sof shareholders, and with sufficient shares the interpolating polynomial p(x)is reconstructed using the Lagrange’s interpolation formula and outputs message m∈Fq as p(0)=m. Other schemes for sharing secrets, such as Blakley’s [14]and Mignotte’s [15], may offer alternative security assumptions, complexity, and performance trade-offs. SSS was chosen for this analysis because it strikes a balance between simplicity, efficiency, security, and versatility, which has led to its widespread adoption and preference across a wide range of applications. Due to its sole reliance on polynomial interpolation, SSS is relatively straightforward to comprehend and implement. In addition, it offers information-theoretic security, also known as perfect security, which means that the scheme’s security is based on mathematical principles rather than computational assumptions and cannot be improved. This makes it exceptionally resistant to various attacks, including those utilizing quantum computers. It must be noted, however, that the fundamental security relies on one’s trust in the dealer’s honesty. If the dealer is dishonest, a share verification mechanism [16] must be used, which is beyond the scope of this paper and will be the subject of future research. The ability to customize the threshold is a crucial component of SSS. As we have seen, the threshold determines the number of shares necessary to reconstruct the secret, allowing for flexibility in establishing the desired level of security and availability. SSS can accommodate any number of participants. In addition, SSS is resistant to share loss and compromise. If the threshold is not exceeded, the secret cannot be reconstructed even if some shares are lost or compromised.Lastly, the protocol’s robustness makes it suitable for situations where the integrity of shares cannot always be guaranteed. 3. POLYNOMIAL INTERPOLATION The primary focus of this article is on the application of polynomial interpolation in SSS over finite fields, aiming to demonstrate that, in specific scenarios, alternative interpolation methods might offer better performance than the conventional Lagrange method. Following a concise overview of polynomial interpolation and its connection to polynomial evaluation, we analyze four interpolation techniques. Our exploration begins with the standard Lagrange interpolation, which Shamir also advocated in his scheme. Subsequently, we introduce two optimizations of the standard method: the optimized Lagrange for the constant term and the well-established Barycentric Lagrange (BL) formula. While elucidating the Lagrange forms, we scrutinize the Newton interpolation and delineate its trade-offs compared with the Lagrange method. The concluding segment of this section will center on the Fast Fourier method, culminating in a succinct assessment of the advantages and limitations of each technique. 3.1. Two ways of polynomial representation Polynomial interpolation finds wide application in approximating complex curves, such as evaluating trigonometric functions or the natural logarithm. Its efficient computation makes it a preferred choice for sub-quadratic multiplication and squaring operations, exemplified by techniques like the Toom-Cook multiplication [17] and the Karatsuba multiplication [18]. In Euclidean geometry, two points determine a line, three points a parabola, and so forth. However, given the variety of approaches available, reconstructing a smooth function f(x)from a set of collected sample values lacks a straightforward answer. A forthcoming analysis aims to establish a mathematical language that bridges the understanding between the reader and the authors. Polynomials are commonly represented in two forms: (i) coefficient-form and (ii) point-value form. Nevertheless, every polynomial in coefficient form has a unique counterpart in point-value form and vice versa. The coefficient representation of the n-degree polynomial p(x)is the vector of coefficients p0,p1,...,pn. Conversely, the point-value representation of the n-degree polynomial p(x)consists of a set of npoint-value pairs, denoted as {(x0,P(x0)),(x1,P(x1)),...,(xn,P(xn))},whereallxiare distinct. For conciseness, let yi=P(xi). Consequently, the shares will be denoted as {(x0,y0),(x1,y1),...,(xn,yn)}. The process of computing the point-value representation of a polynomial P(x), bounded by degree n, given its coefficient representation p∈Rn,isknownasevaluation. By selecting n distinct points x0,x1,...,xn−1, for each xi,wherei=0, 1, ...,n−1, P(xi)can be computed.The evaluation of the value of any arbitrary polynomial P(x)of degree nat a given point xcan be achieved in (n)time, utilizing Horner’s method [19]. This method employs the following formula: P(x)=p0+x(p1+x(p2+ ··· + x(pn−2+ x(pn−1)) ···)). Conversely, if one were to compute the point-value representation of a polynomial P(x)of degree n, it would require (n2)time. On the other hand, the process of computing the coefficient representation of a polynomial P(x),bounded by degree n, given its point-value representation {(x0,y0),(x1,y1),...,(xn,yn)}, is referred to as interpolation. Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 4|Voudouris et al. 3.2. Lagrange interpolation Our examination of interpolation methods commences with the classical Lagrange interpolation technique [20], widely regarded as the most prevalent and extensively utilized method for data interpolation. This method, initially proposed by Shamir in his seminal work, forms the cornerstone of many interpolation procedures. In the following text, we provide a concise overview of its functioning and subsequently introduce two optimizations. At the outset, given a set of n+1points (x0,y0),(x1,y1),..., (xn,yn), the interpolation problem, as described in Section 3.1, entails defining a unique polynomial P(x)of degree nsuch that yi=P(xi),fori=0, 1, ...,n. For the set of points (x0,y0),(x1,y1),...,(xk,yk),wherexj= xmfor all j= m,the Lagrange basis for polynomials of degree ≤kfor those points consists of the basis polynomials 0(x),1(x),...,k(x), each of degree k. These polynomials satisfy j(xm)=0if m= j,orj(xj)=1 if m=j. Each basis polynomial j(x)can be expressed as lj(x)= 0≤m≤k x−xm xj−xm (1) The Lagrange interpolating polynomial for those points through the corresponding values y0,y1,...,ykis the linear combination: L(x)= k  j=0 yjj(x)(2) In this context, P(xi)=yifor all 0 ≤i≤t. Each basis polynomial possesses a degree of k, hence the summation L(x)has a degree of ≤k, ensuring interpolation of the data as L(xm)=k j=0yjj(xm)= k j=0yjδmj =ym. The polynomial P(x)satisfying P(xi)=yifor all 0≤i≤t, as derived from the aforementioned expression L(x)= k j=0yjj(x), is indeed unique. 3.2.1. Optimizations of standard lagrange method The primary objective of interpolation within the context of SSS is to ascertain the value of P(0), representing the secret value. In our comprehension of SSS, it is unnecessary to interpolate the entire polynomial up to the maximum degree, solely the constant term necessitates extraction. Hence, by setting the variable xto zero, we introduce our initial optimization, ptimized Lagrange (Opt.L). Leveraging Lagrange polynomials, a simple optimization technique (achieved through substitutions and basic arithmetic operations), can reduce the number of multiplications by nearly 50%: P(0)=n−1  i=0 −xin−1  i=0 yi −xi(n−1 j=0,j=i(xi−xj)) .(3) The previously mentioned optimization lacks clarity in the existing literature. While numerous schemes employing SSS reference the standard Lagrange interpolation method, practical implementations frequently diverge, balancing security with speed. Hence, we contend that providing developers of SSSbased secret-sharing applications with access to and comprehension of this optimization technique would be invaluable, thus warranting its inclusion in our research. Nevertheless, subsequent sections will demonstrate the existence of superior optimizations. 3.2.2. Barycentric Lagrange formula Continuing with the Lagrange optimizations, the subsequent one is termed the BL formula. Here, we will provide a concise overview, adhering explicitly to the notation and definitions outlined in Berrut et al.[21]. Commencing with the standard Lagrange formula, the numerator of ljin Equation 1can be expressed as (x)=(x−x0)x−x1···x−xn(4) divided by x−xj. By defining the Barycentric weights: wj=1 k=i(xj−xk),j=0, ...,n,(5) we can thus write ljas j(x)=(x)wj x−xj .(6) So, in the end, we define the polynomial P(x)as P(x)=(x) n  j=0 wj x−xj fj.(7) The fundamental contrast between the current optimization and the standard Lagrange interpolation formula lies in the methodology employed for constructing the interpolation polynomial. The BL formula entails expressing the numerator of each Lagrange basis function as the quotient of the product of differences between the interpolation point and all other points, and the difference between the interpolation point and the current point. The determination of Barycentric weights is based on these differences, as shown in Equation 5,which are subsequently utilized to compute the Lagrange basis functions. Eventually, the polynomial is obtained by multiplying the Lagrange basis functions with the function values at the interpolation points while accounting for the Barycentric weights. Section 4will show that the BL method significantly reduces interpolation times. 3.3. Newton interpolation Newton’s interpolation formula provides an alternative approach to compute P(x), given a set of n+1 data points (x0,y0),(x1,y1),..., (xn,yn),wherethexiare distinct. The initial step of this method involves computing Newton’s table of divided differences. This table constitutes a triangular array of numbers defined by the recursive formula: f[xi,xi+1,...,xi+k]=f[xi+1,xi+2, ..., xi+k]−f[xi,xi+1, ..., xi+k−1] xi+k−xi (8) Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 Interpolation techniques applied to SSS |5 with initial conditions f[xi]=f(xi). The table has the following form: x0f[x0] x1f[x1]f[x0,x1] x2f[x2]f[x1,x2]f[x0,x1,x2] . . .. . .. . .. . . xnf[xn]f[xn−1,xn]f[xn−2,xn−1,xn]... f[x0,x1,...,xn] Once the table is computed, the subsequent step involves evaluating P(x)using the following Newton interpolation formula: n k=1f[x0,x1,...,xk](x−x0)···(x−xk−1), where the coefficients f[x0,x1,...,xk] are the entries in the table. This formula can be interpreted as a weighted sum of terms (x−x0)···(x−xk−1), where each term is weighted by the corresponding coefficient f[x0,x1,...,xk]. It should be emphasized that when the interpolation points x0,...,xnare distinct, finding a polynomial passing through the points (xi,yi)is equivalent to solving a system of linear equations Ax=bthat has a unique solution. The matrix Ais determined by the choice of basis for the space of polynomials of degree n or less. In Newton interpolation, the basis functions are the set Nj(x)n j=0defined recursively as N0(x)=1, and for j=1, ...,n: Nj(x)=j−1 k=0(x−xk). These basis functions are structured such that Nj(x)is a polynomial of degree j,andNj(xk)=0for k= 0, ...,j−1. Consequently, any polynomial of degree at most ncan be expressed as a linear combination of these basis functions, i.e. Pn(x)=a0N0(x)+a1N1(x)+···+anNn(x). Updating the interpolating polynomial becomes straightforward as interpolation points are added because the basis functions {Nj(x)},j=0, ...,n, remain unchanged. New interpolation points can be added to update the interpolating polynomial using Newton’s method, and then the divided differences table is updated accordingly. For instance, suppose we have a set of n+1distinct points xi,yi, and we have computed the divided differences table for these points. If we wish to add a new point (xn+1,yn+1), we update the divided differences table by adding a new column to the right. The first entry in this new column will be the divided difference of the new point, which is simply the function value yn+1. Subsequent entries in the column can be recursively computed using divided differences from the previous column. Specifically, the (k+1)th entry in the new column is determined as follows: f[xn+1,x0,···,xk]=f[xn+1,x1,...,xk]−f[xn,x0,...,xk−1] xn+1−xn. Here, f[xn+1,x1,...,xk]represents the divided difference computed in the previous column, while f[xn,x0,...,xk−1] denotes the divided difference from the same row but in the preceding column. Once the divided differences table is updated, the updated interpolating polynomial can be computed using the same formula as before Pn(x)=n k=0f[x0,···,xk]Nk(x),whereNk(x)signifies the kth Newton basis function. 3.4. Fourier interpolation As demonstrated earlier, polynomials can be defined over a field denoted as F. However, in the context of cryptography, we have chosen to define them specifically on GF(p),wherepis a prime. The Fourier transform is typically defined on the infinite set C, which is the field of complex numbers. For the intended scope of our research, it is necessary to establish its definition on GF(p) as well. Subsequently, within this section, we will examine the methodology for achieving this objective through the utilization of the roots of unity on finite fields. This research involves the utilization of the FFT to alter the representation of a polynomial P(x), which has a degree ¡nand consists of coefficients belonging to the field GF(p), from a coefficients-based representation to an evaluation-based representation (and vice versa). The evaluation is conducted exclusively over ndistinct field elements, which are the nth roots of unity. As a result, sufficient knowledge of the preliminaries that were described at the beginning is required. We analyze the DFT, the FFT, and lastly the inverses of both, which will lead us to the interpolation. DFT over finite fields. Consider a polynomial P(x)of degree n expressed in coefficient form, denoted as p∈Fn,wherepis a column vector given by p=[p0p1...pn−1]T. We aim to represent P(x) in its point-value form by computing its values at the ndistinct nth roots of unity, namely ω0 n,ω1 n,...,ωn−1 n.Letyjdenote P(ωj n)for jranging from 0to n−1.Thenth order DFT is a mathematical operation that takes the coefficients of a polynomial Pof degree n−1as input. Its output is an n-dimensional vector, where each component represents the evaluation of Pat one of the nth roots of unity. Conventionally, we denote this operation as y=DFTn(p),wherepsignifies the polynomial of degree n− 1. A visual depiction of this problem can be seen as a matrix multiplication: ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ y0 y1 . . . yn−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ =⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1x0x2 0··· xn−1 0 1x1x2 1··· xn−1 1 . . .. . .. . ..... . . 1xn−1x2 n−1··· xn−1 n−1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ·⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a0 a1 . . . an−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (9) What the DFT essentially provides is a transformation: ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ F0 F1 . . . Fn−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ =⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 11 1 ··· 1 1ωω 2··· ωn−1 . . .. . .. . ..... . . 1ωn−1ω2(n−1)··· ω(n−1)2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ·⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a0 a1 . . . an−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (10) In the transformed relation, observe that xicorresponds to ωi. We designate the Vandermonde matrix of 1, ω,ω2,...,ωn−1as the DFT matrix. Alternatively, we can express this algebraically as Fk= n−1  i=0 aiωki (11) After rearranging, we derive the inverse discrete Fourier transformation (IDFT) formula, designed to recompute the coefficients ai. ai=1 n n−1  k=0 Fkω−ik (12) Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 6|Voudouris et al. This formula arises from the symmetric properties inherent in the primitive nth roots of unity. Similarly, the IDFT can also be expressed through matrix multiplication. ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a0 a1 . . . an−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ =1 n ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 11 1 ··· 1 1ω−1ω−2··· ω−(n−1) . . .. . .. . ..... . . 1ω−(n−1)ω−2(n−1)··· ω−(n−1)2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ·⎛ ⎜ ⎜ ⎜ ⎜ ⎝ F0 F1 . . . Fn−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (13) In line with the DFT matrix, we denote the matrix containing the inverse powers of ωas the IDFT matrix. 3.5. Fast Fourier transform over finite fields We know that Horner’s algorithm [22] facilitates the computation of DFTnin (n2)steps. However, when dealing with significantly large values of n(e.g. hundreds of thousands or millions), executing DFTnmay demand excessive time. An alternative algorithm for computing the DFT of size nin (nlog n)operations, known as the FFT, addresses this challenge. This approach leverages the symmetric properties of the nth roots of unity and is rooted in the technique introduced by Cooley and Tukey [23]. We adhere explicitly to the mathematical notations and definitions outlined by Fateman [24], as they offer a more algebraic perspective and are utilized in our implementations. Fast Fourier transform. The FFT algorithm leverages the symmetry and periodicity properties inherent in roots of unity to compute the Fourier transform more efficiently. By employing the Cooley and Tukey algorithm [23], the complexity is reduced from O(n2)to O(nlog n)operations. The fundamental concept involves breaking down the initial DFT problem into smaller DFT problems of half the size, employing a recursive divideand-conquer approach. This characteristic makes the FFT particularly valuable for handling large datasets and real-time applications. For demonstration purposes and to showcase the algorithm’s optimal performance, we opt for the number of shares nto be a power of 2, i.e. n=2k. While the same methodology can be extended to other values of n,suchas n=pkor n=p1p2p3...pk, these choices may not underscore the significant advantage of the Cooley and Tukey algorithm as effectively. Consider the polynomial Pdefined as P(x)=n−1 i=0aixi.We can partition Pinto two polynomials: one comprising the evenindexed coefficients of P, denoted as P0, and another consisting of the odd-indexed coefficients, denoted as P1. P(x)=(a0+a2x2+... an−2xn−2)   P0 +(a1x+a3x3+... +an−1xn−1)   P1 (14) P(x)= n 2−1  i=0 a2ix2i    P0 + n 2−1  i=0 a2i+1x2i+1    P1 (15) P(x)= n 2−1  i=0 a2i(x2)i    P0 +x n 2−1  i=0 a2i+1(x2)i    P1 (16) Notice that the first term in the summation represents the evaluation of the (n/2)-degree polynomial P0at x2,whilethe second term represents the evaluation of the (n/2)-degree polynomial P1at x2. Consequently, we can express Pas P(x)=P0(x2)+xP1(x2)(17) We iterate this procedure for each of P0and P1until the size of the FFT reduces to 1 and can no longer be split further. By recursively evaluating the smaller sized FFTs and substituting these evaluations into the polynomial of the higher layer, we eventually derive the shares of the original polynomial. The inverse Fourier transform is the problem of determining the left-hand side coefficient vector of the following relation: ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a0 a1 . . . an−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ =1 n ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 11 1 ··· 1 1ω−1ω−2··· ω−(n−1) . . .. . .. . ..... . . 1ω−(n−1)ω−2(n−1)··· ω−(n−1)2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ·⎛ ⎜ ⎜ ⎜ ⎜ ⎝ F0 F1 . . . Fn−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ (18) The algorithm resembles the standard Fourier Transform, but in this case, the required matrix comprises negative powers of ω. Thus, the objective is to apply the FFT as described above, but this time evaluating at negative powers. As a final step, the inverse Fourier transform divides the inner product of each multiplied vector by nto obtain the coefficients. Section 4presents our testbed environment and the rationale behind selecting the respective parameters. 4. EVALUATION After establishing the essential mathematical groundwork and conducting a thorough analysis of each interpolation method, this section encapsulates the culmination of our experiments. The experimental results demonstrate how various parameters influence the efficiency and performance of each interpolation method. This contribution is motivated by the lack of related work on this topic, particularly when combined with secret sharing. To underscore the feasibility and practicality of our approach, we implemented a prototype using Python 3.11.18 and SageMath 10.3 [9] on a MacBook running macOS Sonoma 14.4. 4.1. Domain parameters configuration Considering that our interpolation techniques will be employed for SSS across finite fields, defining the finite field within the specified domain we will be working with is crucial. The domain parameters encompass a set of predefined characteristics that hold universal applicability and will be consistently utilized across all experimental procedures. These parameters are influenced by the constraints of the FFT in the following manner. The initial step involves constructing the finite field Fp. Since the FFT does not function with a randomly chosen p,wemust manually construct it. To compute the primitive nth roots of unity, the division p−1 nmust yield an integer, indicating that ndivides p−1.Weoptforpto be a safe prime following the format p= nmaxq+1,whereqrepresents a random 28-bit prime number. Here, we set nmax to 2100 as an upper bound since, in our specific use Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 Interpolation techniques applied to SSS |7 Table 1. Performance of unoptimized Newton, Barycentric, and FFT. Methods 128 256 512 1024 2048 4096 8192 Newton 0.05 0.12 0.6445 1.9967 8.8905 31.8654 320.3828 Barycentric 0.1483 0.3906 0.9381 2.998 12.5994 38.0322 114.7491 FFT 0.0025 0.0081 0.0112 0.0344 0.1144 0.1393 0.2699 cases, nwill never reach such magnitudes. Hence, the value of p is determined as p=2100 ∗135783853 +1 =172126482756751667503106328023403593729 (19) Thesizeofpis arbitrary as it does not affect the performance of the interpolation methods. Another important domain parameter crucial mostly for FFT is the group generator g. We can find gusing any classical method for determining generators, but for simplicity, as stated by Lynn et al.[ 25], we let g=3. For each specific scenario, a 128-bit AES key will be consistently designated as the constant term of the polynomial and treated as the secret value s. The selection of a 128-bit key aligns with the length of the prime number p(Equation (19)). While the algorithm could support longer keys, such as those with 256 or 512 bits, maintaining equivalent security necessitates proportionally increasing the length of the prime number. Deviating from this balance by selecting a shorter prime number or extending the key length beyond what is proportional compromises security. Shorter primes risk insecurity due to modular arithmetic properties and loss of information, while longer keys necessitate padding. The objective is to distribute this confidential information among participants of the secret-sharing protocol to securely reach a consensus on the same key. 4.2. A comparative analysis of different interpolation methods Regarding our initial experiment, certain adjustments are necessary to accommodate the constraints of the FFT algorithm. 1. The number of shares nshould be a power of 2. Specifically, we set nmin =27and nmax =213 as these selections effectively highlight performance differences. 2. Maintaining consistency in the value of the primitive root, ω, across various nnecessitates minor adjustments, which are not aligned with our objectives. While precomputing the roots of unity could potentially resolve this issue, our current approach utilizes unoptimized versions of each algorithm. 3. Due to the considerations above, it is imperative that the degree of the interpolating polynomial remains n−1 for all values of n. Consequently, it is not possible to determine a fixed interpolating polynomial, necessitating the generation of an arbitrary polynomial of degree (n−1)for each value of n. We will now analyze the results obtained from the unoptimized implementations of the Newton, Barycentric, and FFT algorithms. Each cell within the Table 1 denotes the duration, measured in seconds, required to complete the interpolation process. It is worth noting that, based on the FFT algorithm’s computational complexity of O(nlog2n), its performance outperforms by far all the other interpolation methods. Figure 1. Interpolation times for unoptimized Newton, barycentric, and FFT diagram. Table 2. Performance of unoptimized naive Lagrange. Number of users Interpolation time (s) 100 0.5433 200 3.1385 300 6.5299 400 16.1861 500 25.2504 600 49.3655 700 64.2438 800 80.3796 900 107.0505 1000 152.7803 We notice that for values of nup to 4096, both Newton and Barycentric methods demonstrate nearly comparable performance. However, Newton’s time complexity undergoes a notable and rapid increase beyond this threshold. The disparity between the computational complexities of O(nlog2(n)) and O(n2)becomes evident in Fig. 1, as demonstrated by the logarithmic nature of the FFT graph, which converges toward a constant O(1)with an infinitesimally small rate of change. In the second segment of our initial experiment (i.e. the standard Lagrange method), we set the lower bound nmin to 100 and the upper bound nmax to 1000, incrementing the value of nby 100 in each successive iteration. Due to the absence of constraints associated with the FFT, it becomes feasible to establish a deterministic solution for a randomly selected polynomial of degree five. The outcomes are depicted in Table 2. It is evident that for n=1000, the Lagrange algorithm necessitates ∼2 min to complete, as depicted in Fig. 2. Naive Lagrange interpolation time increases quickly even with a small number of shares (n=1000), as anticipated due to its quadratic time complexity O(n2). Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 8|Voudouris et al. Table 3. Performance of Newton and Barycentric without precomputed values. Methods 5000 5001 5002 5003 5004 5005 5006 5007 5008 5009 5010 Newton 8.5294 9.376 8.7036 9.3542 9.3054 9.5331 8.9018 8.8316 10.1682 9.0897 8.8205 Barycentric 5.6065 5.5197 5.7283 5.6771 5.6474 5.6487 5.7032 5.5013 5.6746 5.7882 5.524 Figure 2. Standard Lagrange method for increasing number of members. 4.3. Enhancing Newton and barycentric methods with new shares This experiment aims to explore an application employing secret sharing for key exchange. Both the enhanced Newton and Barycentric methods utilize precomputed divided differences and weights. These precomputed values are reused for each new user added to the system, eliminating the need for recalculations from scratch. To illustrate this concept, we set the variable representing the number of shares to n=5000. The precomputed values are stored and gradually expanded with each new share until the value of n reaches 5010. Let us first examine the scenario where no values have been precomputed. Table 3 exhibits the results obtained from the Newton and Barycentric methods when the precomputation of divided differences and weights is not conducted. It is crucial to note that the utilized versions of these approaches are optimized, focusing solely on recovering the constant term a0 rather than the entire polynomial. Each cell denotes the duration in seconds required to complete the interpolation process. The advantage of Barycentric over Newton is evident from the very beginning with the n=5000 shares. The number of shares is selected arbitrarily and is intentionally set to a large value to represent a real-life scenario. This is in contrast to the n=1000 used in the naive Lagrange method, which is not suitable for high values of nand has no impact on the result. Moreover, the time advantage persists even when incorporating additional shares. Bezzateev et al.[ 6] proposed that using Newton interpolation for adding a single new share in the multifactor authentication setting is beneficial, compared with naive Lagrange. While this is correct, we bolster this assumption by demonstrating that the Barycentric technique is more advantageous, even for a larger number of new shares. In the subsequent experiments, we showcase that the FFT outperforms all methods with the correct setup. Now, let us examine the outcomes of the Newton and Barycentric methods when utilizing precomputed values with the same number of shares. Table 4 highlights the significant difference in time complexity when divided differences and weights are precomputed. For the initial value n=5000, the time in the first column is nearly identical to Table 3, as the experiment is repeated to establish the initial divided differences and weights. However, subsequent interpolations with n+1shares show notable changes in interpolation time. This is because interpolations with n+1 shares utilize the precomputed ndivided differences and weight values, while interpolations with n+2shares utilize the precomputed n+1 values, and so forth. Both methods demonstrate substantial improvements in interpolation times compared with Table 3, with the Barycentric method still maintaining a computational advantage over Newton. It can be concluded that irrespective of the application, employing Barycentric interpolation is more advantageous for a secret-sharing application experiencing a progressive increase in the number of users. 4.4. Comparison of all optimized interpolation methods In our final experiment, we compare all the aforementioned interpolation methods. It should be noted that the restrictions of the FFT method are considered, such as the requirement for the number of shares nto be a power of 2. All approaches are optimized by precomputing values for divided differences in Newton’s method, weights in BL’s method, and roots of unity in the FFT method. This experiment aims to determine the most effective interpolation method for standard SSS based only on performance metrics such as interpolation time and shared secret retrieval a0 among nusers. Table 5 presents interpolation times in seconds for varying numbers of shares. It is worth noting that the FFT method does not even surpass 1 s for n=16 384 shares. Furthermore, Fig. 3 illustrates all optimal variations of the methodologies. Computational bounds are set at nmin =27and nmax =213. The FFT approach demonstrates almost a straight line on the x−yplane due to its time complexity of O(nlog2(n)), which is significantly lower than the O(n2)complexity of other methods, except Barycentric, which has an O(n)retrieval time with precomputed weights. 5. DISCUSSION Through an examination of the research questions, we conducted an in-depth investigation into the practicalities and relative efficiencies of different polynomial interpolation techniques as they apply to secret-sharing schemes and more extensive realworld implementations. This research has been thorough, encompassing both theoretical frameworks and empirical validations through carefully designed experiments. Our journey through this analytical process illuminated the distinct capabilities, limitations, and situational advantages of the FFT, the BL formula, and Newton’s interpolation method. Every technique was evaluated in terms of its computational complexity, practicability in finite fields, and ease of implementation in dynamic environments. Our main objective while analyzing these methodologies was Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024 Interpolation techniques applied to SSS |9 Table 4. Performance of Newton and Barycentric with precomputed values. Methods 5000 5001 5002 5003 5004 5005 5006 5007 5008 5009 5010 Newton 8.557 0.0085 0.0086 0.0087 0.3562 0.0082 0.008 0.0083 0.008 0.0083 0.0083 Barycentric 5.7163 0.0091 0.0092 0.0088 0.0086 0.0092 0.009 0.009 0.0095 0.0093 0.0098 Table 5. Performance of FFT and optimized standard Lagrange, Newton, and Barycentric. Methods 128 256 512 1024 2048 4096 8192 Standard Lagrange 0.0112 0.0465 0.1719 0.6817 2.7803 10.9333 43.7564 Newton 0.0322 0.0194 0.1042 0.4909 1.6527 6.1578 27.8844 Barycentric 0.004 0.0157 0.0615 0.246 0.9626 3.7881 15.1438 FFT 0.001 0.0019 0.0044 0.0092 0.019 0.0405 0.0864 Figure 3. Optimized methods performance. to unequivocally respond to the research questions that were introduced at the beginning of this paper. Having reached the final phase of our investigation, equipped with a comprehensive body of evidence and insights derived from our rigorous analysis, we confidently conclude and respond to the initial research questions. This serves as concrete proof of our dedication to expanding the knowledge of the research community into applying polynomial interpolation methods into practical, largescale applications, in addition to being an endorsement of the depth of our investigation. RQ1:What is the swiftest method and most suitable in large-scale scenarios to interpolate the polynomial in the context of SSS to recover the secret from the constant term? The FFT emerges as the premier choice for polynomial evaluation and interpolation in SSS scheme, due to its unparalleled computational efficiency and rapid execution times, particularly when handling large datasets. Shamir’s scheme relies on constructing a polynomial of degree k−1for a k-threshold scheme, where the secret is embedded in the constant term. For reconstruction, traditionally, Lagrange interpolation is used, which has a computational complexity of O(n2)for nshares. In contrast, FFT can perform the necessary polynomial evaluations or interpolations in O(nlog n)time. This efficiency is largely due to the FFT’s ability to perform polynomial operations in O(nlogn)time, leveraging divide-and-conquer strategies [23] and the discrete Fourier transform’s symmetry properties. Such a significant reduction in computational complexity positions the FFT as exceptionally advantageous for secret reconstruction within large-scale operations. Nonetheless, the FFT requires the number of participants and the polynomial’s degree to adhere to powers of 2, a constraint that can be circumvented through algorithmic enhancements designed to broaden its applicability [26]. On the other hand, the Lagrange interpolation method and its evolution through the BL formula offer a different set of advantages and limitations. While the Lagrange method is renowned for its versatility, it faces practical difficulties in dynamically changing scenarios, such as the addition of new participants to a sharing scheme. This limitation arises from the need to recalculate Lagrange polynomials from scratch with each new addition, rendering the process cumbersome in rapidly evolving settings. The BL formula, in contrast, addresses some of the practical challenges inherent in the standard Lagrange approach by eliminating the need for complete polynomial recalculations for every new interpolation point. This refinement significantly enhances efficiency, particularly in large datasets, by reducing the overall computational demand to O(n)for evaluating interpolating polynomials once Barycentric weights are computed. The primary challenge of BL method lies in managing small datasets where the overhead of computing Barycentric weights may outweigh its computational benefits. Furthermore, Newton’s interpolation method, while sharing foundational similarities with the Lagrange approach in terms of employing a specific set of basis functions, distinguishes itself through the utilization of divided difference polynomials. This method offers advantages in computational efficiency for iterative polynomial evaluations, albeit at the cost of increased memory complexity due to the storage requirements of the divided difference table. Despite these advantages, Newton’s method is not the primary choice compared with the FFT and BL methods, especially when quick interpolation of polynomials at new points is required. In delving deeper into the comparative analysis of these methods, it becomes evident that the FFT stands out for its robust applicability and superior efficiency in large-scale secret-sharing schemes, assuming its specific constraints are effectively managed.The BL formula, with its focus on overcoming the limitations of the Lagrange method in dynamic and large datasets, provides a complementary approach that excels under certain conditions, particularly in environments not restricted by the size of the dataset. Consequently, when assessing the most suitable interpolation method for large-scale scenarios, the FFT unequivocally offers the most promising solution, subject to further refinements and adaptations to extend its utility beyond the confines of Downloaded from https://academic.oup.com/comjnl/advance-article/doi/10.1093/comjnl/bxae109/7900561 by guest on 02 December 2024